MFEM  v4.6.0 Finite element discretization library
ex37p.cpp
Go to the documentation of this file.
1 // MFEM Example 37 - Parallel Version
2 //
3 // Compile with: make ex37p
4 //
5 // Sample runs:
6 // mpirun -np 4 ex37p -alpha 10 -pv
7 // mpirun -np 4 ex37p -lambda 0.1 -mu 0.1
8 // mpirun -np 4 ex37p -o 2 -alpha 5.0 -mi 50 -vf 0.4 -ntol 1e-5
9 // mpirun -np 4 ex37p -r 6 -o 2 -alpha 10.0 -epsilon 0.02 -mi 50 -ntol 1e-5
10 //
11 // Description: This example code demonstrates the use of MFEM to solve a
12 // density-filtered [3] topology optimization problem. The
13 // objective is to minimize the compliance
14 //
15 // minimize ∫_Ω f⋅u dx over u ∈ [H¹(Ω)]² and ρ ∈ L¹(Ω)
16 //
17 // subject to
18 //
19 // -Div(r(ρ̃)Cε(u)) = f in Ω + BCs
20 // -ϵ²Δρ̃ + ρ̃ = ρ in Ω + Neumann BCs
21 // 0 ≤ ρ ≤ 1 in Ω
22 // ∫_Ω ρ dx = θ vol(Ω)
23 //
24 // Here, r(ρ̃) = ρ₀ + ρ̃³ (1-ρ₀) is the solid isotropic material
25 // penalization (SIMP) law, C is the elasticity tensor for an
26 // isotropic linearly elastic material, ϵ > 0 is the design
27 // length scale, and 0 < θ < 1 is the volume fraction.
28 //
29 // The problem is discretized and gradients are computing using
30 // finite elements [1]. The design is optimized using an entropic
31 // mirror descent algorithm introduced by Keith and Surowiec [2]
32 // that is tailored to the bound constraint 0 ≤ ρ ≤ 1.
33 //
34 // This example highlights the ability of MFEM to deliver high-
35 // order solutions to inverse design problems and showcases how
36 // to set up and solve PDE-constrained optimization problems
37 // using the so-called reduced space approach.
38 //
39 // [1] Andreassen, E., Clausen, A., Schevenels, M., Lazarov, B. S., & Sigmund, O.
40 // (2011). Efficient topology optimization in MATLAB using 88 lines of
41 // code. Structural and Multidisciplinary Optimization, 43(1), 1-16.
42 // [2] Keith, B. and Surowiec, T. (2023) Proximal Galerkin: A structure-
43 // preserving finite element method for pointwise bound constraints.
44 // arXiv:2307.12444 [math.NA]
45 // [3] Lazarov, B. S., & Sigmund, O. (2011). Filters in topology optimization
46 // based on Helmholtz‐type differential equations. International Journal
47 // for Numerical Methods in Engineering, 86(6), 765-781.
48
49 #include "mfem.hpp"
50 #include <iostream>
51 #include <fstream>
52 #include "ex37.hpp"
53
54 using namespace std;
55 using namespace mfem;
56
57 /**
58  * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
59  * ∫_Ω ρ dx = θ vol(Ω) as follows:
60  *
61  * 1. Compute the root of the R → R function
62  * f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
63  * 2. Set ψ ← ψ + c.
64  *
65  * @param psi a GridFunction to be updated
66  * @param target_volume θ vol(Ω)
67  * @param tol Newton iteration tolerance
68  * @param max_its Newton maximum iteration number
69  * @return double Final volume, ∫_Ω sigmoid(ψ)
70  */
71 double proj(ParGridFunction &psi, double target_volume, double tol=1e-12,
72  int max_its=10)
73 {
74  MappedGridFunctionCoefficient sigmoid_psi(&psi, sigmoid);
75  MappedGridFunctionCoefficient der_sigmoid_psi(&psi, der_sigmoid);
76
77  ParLinearForm int_sigmoid_psi(psi.ParFESpace());
79  ParLinearForm int_der_sigmoid_psi(psi.ParFESpace());
81  der_sigmoid_psi));
82  bool done = false;
83  for (int k=0; k<max_its; k++) // Newton iteration
84  {
85  int_sigmoid_psi.Assemble(); // Recompute f(c) with updated ψ
86  double f = int_sigmoid_psi.Sum();
87  MPI_Allreduce(MPI_IN_PLACE, &f, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
88  f -= target_volume;
89
90  int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
91  double df = int_der_sigmoid_psi.Sum();
92  MPI_Allreduce(MPI_IN_PLACE, &df, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
93
94  const double dc = -f/df;
95  psi += dc;
96  if (abs(dc) < tol) { done = true; break; }
97  }
98  if (!done)
99  {
100  mfem_warning("Projection reached maximum iteration without converging. "
101  "Result may not be accurate.");
102  }
103  int_sigmoid_psi.Assemble();
104  double material_volume = int_sigmoid_psi.Sum();
105  MPI_Allreduce(MPI_IN_PLACE, &material_volume, 1, MPI_DOUBLE, MPI_SUM,
106  MPI_COMM_WORLD);
107  return material_volume;
108 }
109
110 /**
111  * ---------------------------------------------------------------
112  * ALGORITHM PREAMBLE
113  * ---------------------------------------------------------------
114  *
115  * The Lagrangian for this problem is
116  *
117  * L(u,ρ,ρ̃,w,w̃) = (f,u) - (r(ρ̃) C ε(u),ε(w)) + (f,w)
118  * - (ϵ² ∇ρ̃,∇w̃) - (ρ̃,w̃) + (ρ,w̃)
119  *
120  * where
121  *
122  * r(ρ̃) = ρ₀ + ρ̃³ (1 - ρ₀) (SIMP rule)
123  *
124  * ε(u) = (∇u + ∇uᵀ)/2 (symmetric gradient)
125  *
126  * C e = λtr(e)I + 2μe (isotropic material)
127  *
128  * NOTE: The Lame parameters can be computed from Young's modulus E
129  * and Poisson's ratio ν as follows:
130  *
131  * λ = E ν/((1+ν)(1-2ν)), μ = E/(2(1+ν))
132  *
133  * ---------------------------------------------------------------
134  *
135  * Discretization choices:
136  *
137  * u ∈ V ⊂ (H¹)ᵈ (order p)
138  * ψ ∈ L² (order p - 1), ρ = sigmoid(ψ)
139  * ρ̃ ∈ H¹ (order p)
140  * w ∈ V (order p)
141  * w̃ ∈ H¹ (order p)
142  *
143  * ---------------------------------------------------------------
144  * ALGORITHM
145  * ---------------------------------------------------------------
146  *
147  * Update ρ with projected mirror descent via the following algorithm.
148  *
149  * 1. Initialize ψ = inv_sigmoid(vol_fraction) so that ∫ sigmoid(ψ) = θ vol(Ω)
150  *
151  * While not converged:
152  *
153  * 2. Solve filter equation ∂_w̃ L = 0; i.e.,
154  *
155  * (ϵ² ∇ ρ̃, ∇ v ) + (ρ̃,v) = (ρ,v) ∀ v ∈ H¹.
156  *
157  * 3. Solve primal problem ∂_w L = 0; i.e.,
158  *
159  * (λ r(ρ̃) ∇⋅u, ∇⋅v) + (2 μ r(ρ̃) ε(u), ε(v)) = (f,v) ∀ v ∈ V.
160  *
161  * NB. The dual problem ∂_u L = 0 is the negative of the primal problem due to symmetry.
162  *
163  * 4. Solve for filtered gradient ∂_ρ̃ L = 0; i.e.,
164  *
165  * (ϵ² ∇ w̃ , ∇ v ) + (w̃ ,v) = (-r'(ρ̃) ( λ |∇⋅u|² + 2 μ |ε(u)|²),v) ∀ v ∈ H¹.
166  *
167  * 5. Project the gradient onto the discrete latent space; i.e., solve
168  *
169  * (G,v) = (w̃,v) ∀ v ∈ L².
170  *
171  * 6. Bregman proximal gradient update; i.e.,
172  *
173  * ψ ← ψ - αG + c,
174  *
175  * where α > 0 is a step size parameter and c ∈ R is a constant ensuring
176  *
177  * ∫_Ω sigmoid(ψ - αG + c) dx = θ vol(Ω).
178  *
179  * end
180  */
181
182 int main(int argc, char *argv[])
183 {
184  // 0. Initialize MPI and HYPRE.
185  Mpi::Init();
186  int num_procs = Mpi::WorldSize();
187  int myid = Mpi::WorldRank();
188  Hypre::Init();
189
190  // 1. Parse command-line options.
191  int ref_levels = 5;
192  int order = 2;
193  double alpha = 1.0;
194  double epsilon = 0.01;
195  double vol_fraction = 0.5;
196  int max_it = 1e3;
197  double itol = 1e-1;
198  double ntol = 1e-4;
199  double rho_min = 1e-6;
200  double lambda = 1.0;
201  double mu = 1.0;
202  bool glvis_visualization = true;
203  bool paraview_output = false;
204
205  OptionsParser args(argc, argv);
207  "Number of times to refine the mesh uniformly.");
209  "Order (degree) of the finite elements.");
211  "Step length for gradient descent.");
213  "Length scale for ρ.");
215  "Maximum number of gradient descent iterations.");
217  "Normalized exit tolerance.");
219  "Increment exit tolerance.");
221  "Volume fraction for the material density.");
223  "Lamé constant λ.");
225  "Lamé constant μ.");
227  "Minimum of density coefficient.");
229  "--no-visualization",
230  "Enable or disable GLVis visualization.");
232  "--no-paraview",
233  "Enable or disable ParaView output.");
234  args.Parse();
235  if (!args.Good())
236  {
237  if (myid == 0)
238  {
239  args.PrintUsage(cout);
240  }
241  MPI_Finalize();
242  return 1;
243  }
244  if (myid == 0)
245  {
246  mfem::out << num_procs << " number of process created.\n";
247  args.PrintOptions(cout);
248  }
249
250  Mesh mesh = Mesh::MakeCartesian2D(3, 1, mfem::Element::Type::QUADRILATERAL,
251  true, 3.0, 1.0);
252  int dim = mesh.Dimension();
253
254  // 2. Set BCs.
255  for (int i = 0; i<mesh.GetNBE(); i++)
256  {
257  Element * be = mesh.GetBdrElement(i);
258  Array<int> vertices;
259  be->GetVertices(vertices);
260
261  double * coords1 = mesh.GetVertex(vertices[0]);
262  double * coords2 = mesh.GetVertex(vertices[1]);
263
264  Vector center(2);
265  center(0) = 0.5*(coords1[0] + coords2[0]);
266  center(1) = 0.5*(coords1[1] + coords2[1]);
267
268  if (abs(center(0) - 0.0) < 1e-10)
269  {
270  // the left edge
271  be->SetAttribute(1);
272  }
273  else
274  {
275  // all other boundaries
276  be->SetAttribute(2);
277  }
278  }
279  mesh.SetAttributes();
280
281  // 3. Refine the mesh.
282  for (int lev = 0; lev < ref_levels; lev++)
283  {
284  mesh.UniformRefinement();
285  }
286
287  ParMesh pmesh(MPI_COMM_WORLD, mesh);
288  mesh.Clear();
289
290  // 4. Define the necessary finite element spaces on the mesh.
291  H1_FECollection state_fec(order, dim); // space for u
292  H1_FECollection filter_fec(order, dim); // space for ρ̃
293  L2_FECollection control_fec(order-1, dim,
294  BasisType::GaussLobatto); // space for ψ
295  ParFiniteElementSpace state_fes(&pmesh, &state_fec,dim);
296  ParFiniteElementSpace filter_fes(&pmesh, &filter_fec);
297  ParFiniteElementSpace control_fes(&pmesh, &control_fec);
298
299  HYPRE_BigInt state_size = state_fes.GlobalTrueVSize();
300  HYPRE_BigInt control_size = control_fes.GlobalTrueVSize();
301  HYPRE_BigInt filter_size = filter_fes.GlobalTrueVSize();
302  if (myid==0)
303  {
304  cout << "Number of state unknowns: " << state_size << endl;
305  cout << "Number of filter unknowns: " << filter_size << endl;
306  cout << "Number of control unknowns: " << control_size << endl;
307  }
308
309  // 5. Set the initial guess for ρ.
310  ParGridFunction u(&state_fes);
311  ParGridFunction psi(&control_fes);
312  ParGridFunction psi_old(&control_fes);
313  ParGridFunction rho_filter(&filter_fes);
314  u = 0.0;
315  rho_filter = vol_fraction;
316  psi = inv_sigmoid(vol_fraction);
317  psi_old = inv_sigmoid(vol_fraction);
318
319  // ρ = sigmoid(ψ)
321  // Interpolation of ρ = sigmoid(ψ) in control fes (for ParaView output)
322  ParGridFunction rho_gf(&control_fes);
323  // ρ - ρ_old = sigmoid(ψ) - sigmoid(ψ_old)
324  DiffMappedGridFunctionCoefficient succ_diff_rho(&psi, &psi_old, sigmoid);
325
326  // 6. Set-up the physics solver.
327  int maxat = pmesh.bdr_attributes.Max();
328  Array<int> ess_bdr(maxat);
329  ess_bdr = 0;
330  ess_bdr[0] = 1;
331  ConstantCoefficient one(1.0);
332  ConstantCoefficient lambda_cf(lambda);
333  ConstantCoefficient mu_cf(mu);
334  LinearElasticitySolver * ElasticitySolver = new LinearElasticitySolver();
335  ElasticitySolver->SetMesh(&pmesh);
336  ElasticitySolver->SetOrder(state_fec.GetOrder());
337  ElasticitySolver->SetupFEM();
338  Vector center(2); center(0) = 2.9; center(1) = 0.5;
339  Vector force(2); force(0) = 0.0; force(1) = -1.0;
340  double r = 0.05;
341  VolumeForceCoefficient vforce_cf(r,center,force);
342  ElasticitySolver->SetRHSCoefficient(&vforce_cf);
343  ElasticitySolver->SetEssentialBoundary(ess_bdr);
344
345  // 7. Set-up the filter solver.
347  DiffusionSolver * FilterSolver = new DiffusionSolver();
348  FilterSolver->SetMesh(&pmesh);
349  FilterSolver->SetOrder(filter_fec.GetOrder());
350  FilterSolver->SetDiffusionCoefficient(&eps2_cf);
351  FilterSolver->SetMassCoefficient(&one);
352  Array<int> ess_bdr_filter;
353  if (pmesh.bdr_attributes.Size())
354  {
355  ess_bdr_filter.SetSize(pmesh.bdr_attributes.Max());
356  ess_bdr_filter = 0;
357  }
358  FilterSolver->SetEssentialBoundary(ess_bdr_filter);
359  FilterSolver->SetupFEM();
360
361  ParBilinearForm mass(&control_fes);
363  mass.Assemble();
364  HypreParMatrix M;
365  Array<int> empty;
366  mass.FormSystemMatrix(empty,M);
367
368  // 8. Define the Lagrange multiplier and gradient functions.
370  ParGridFunction w_filter(&filter_fes);
371
372  // 9. Define some tools for later.
373  ConstantCoefficient zero(0.0);
374  ParGridFunction onegf(&control_fes);
375  onegf = 1.0;
376  ParGridFunction zerogf(&control_fes);
377  zerogf = 0.0;
378  ParLinearForm vol_form(&control_fes);
380  vol_form.Assemble();
381  double domain_volume = vol_form(onegf);
382  const double target_volume = domain_volume * vol_fraction;
383
384  // 10. Connect to GLVis. Prepare for VisIt output.
385  char vishost[] = "localhost";
386  int visport = 19916;
387  socketstream sout_r;
388  if (glvis_visualization)
389  {
390  sout_r.open(vishost, visport);
391  sout_r.precision(8);
392  }
393
394  mfem::ParaViewDataCollection paraview_dc("ex37p", &pmesh);
395  if (paraview_output)
396  {
397  rho_gf.ProjectCoefficient(rho);
398  paraview_dc.SetPrefixPath("ParaView");
399  paraview_dc.SetLevelsOfDetail(order);
400  paraview_dc.SetDataFormat(VTKFormat::BINARY);
401  paraview_dc.SetHighOrderOutput(true);
402  paraview_dc.SetCycle(0);
403  paraview_dc.SetTime(0.0);
404  paraview_dc.RegisterField("displacement",&u);
405  paraview_dc.RegisterField("density",&rho_gf);
406  paraview_dc.RegisterField("filtered_density",&rho_filter);
407  paraview_dc.Save();
408  }
409
410  // 11. Iterate:
411  for (int k = 1; k <= max_it; k++)
412  {
413  if (k > 1) { alpha *= ((double) k) / ((double) k-1); }
414
415  if (myid == 0)
416  {
417  cout << "\nStep = " << k << endl;
418  }
419
420  // Step 1 - Filter solve
421  // Solve (ϵ^2 ∇ ρ̃, ∇ v ) + (ρ̃,v) = (ρ,v)
422  FilterSolver->SetRHSCoefficient(&rho);
423  FilterSolver->Solve();
424  rho_filter = *FilterSolver->GetFEMSolution();
425
426  // Step 2 - State solve
427  // Solve (λ r(ρ̃) ∇⋅u, ∇⋅v) + (2 μ r(ρ̃) ε(u), ε(v)) = (f,v)
428  SIMPInterpolationCoefficient SIMP_cf(&rho_filter,rho_min, 1.0);
429  ProductCoefficient lambda_SIMP_cf(lambda_cf,SIMP_cf);
430  ProductCoefficient mu_SIMP_cf(mu_cf,SIMP_cf);
431  ElasticitySolver->SetLameCoefficients(&lambda_SIMP_cf,&mu_SIMP_cf);
432  ElasticitySolver->Solve();
433  u = *ElasticitySolver->GetFEMSolution();
434
435  // Step 3 - Adjoint filter solve
436  // Solve (ϵ² ∇ w̃, ∇ v) + (w̃ ,v) = (-r'(ρ̃) ( λ |∇⋅u|² + 2 μ |ε(u)|²),v)
437  StrainEnergyDensityCoefficient rhs_cf(&lambda_cf,&mu_cf,&u, &rho_filter,
438  rho_min);
439  FilterSolver->SetRHSCoefficient(&rhs_cf);
440  FilterSolver->Solve();
441  w_filter = *FilterSolver->GetFEMSolution();
442
443  // Step 4 - Compute gradient
444  // Solve G = M⁻¹w̃
445  GridFunctionCoefficient w_cf(&w_filter);
446  ParLinearForm w_rhs(&control_fes);
448  w_rhs.Assemble();
450
451  // Step 5 - Update design variable ψ ← proj(ψ - αG)
453  const double material_volume = proj(psi, target_volume);
454
455  // Compute ||ρ - ρ_old|| in control fes.
456  double norm_increment = zerogf.ComputeL1Error(succ_diff_rho);
458  psi_old = psi;
459
460  double compliance = (*(ElasticitySolver->GetLinearForm()))(u);
461  MPI_Allreduce(MPI_IN_PLACE,&compliance,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
462  if (myid == 0)
463  {
464  mfem::out << "norm of the reduced gradient = " << norm_reduced_gradient << endl;
465  mfem::out << "norm of the increment = " << norm_increment << endl;
466  mfem::out << "compliance = " << compliance << endl;
467  mfem::out << "volume fraction = " << material_volume / domain_volume << endl;
468  }
469
470  if (glvis_visualization)
471  {
472  ParGridFunction r_gf(&filter_fes);
473  r_gf.ProjectCoefficient(SIMP_cf);
474  sout_r << "parallel " << num_procs << " " << myid << "\n";
475  sout_r << "solution\n" << pmesh << r_gf
476  << "window_title 'Design density r(ρ̃)'" << flush;
477  }
478
479  if (paraview_output)
480  {
481  rho_gf.ProjectCoefficient(rho);
482  paraview_dc.SetCycle(k);
483  paraview_dc.SetTime((double)k);
484  paraview_dc.Save();
485  }
486
487  if (norm_reduced_gradient < ntol && norm_increment < itol)
488  {
489  break;
490  }
491  }
492
493  delete ElasticitySolver;
494  delete FilterSolver;
495
496  return 0;
497 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
int visport
GridFunction * GetFEMSolution()
Definition: ex37.hpp:716
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetDataFormat(VTKFormat fmt)
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
LinearForm * GetLinearForm()
Definition: ex37.hpp:353
double epsilon
Definition: ex25.cpp:141
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
Helper class for ParaView visualization data.
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void SetEssentialBoundary(const Array< int > &ess_bdr_)
Definition: ex37.hpp:344
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition: fe_coll.hpp:225
Scalar coefficient defined as the product of two scalar coefficients or a scalar and a scalar coeffic...
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:561
STL namespace.
virtual double ComputeL1Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:270
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1089
double inv_sigmoid(double x)
Inverse sigmoid function.
Definition: ex37.hpp:12
double der_sigmoid(double x)
Derivative of sigmoid function.
Definition: ex37.hpp:33
Integrator that inverts the matrix assembled by another integrator.
Definition: bilininteg.hpp:371
void SetEssentialBoundary(const Array< int > &ess_bdr_)
Definition: ex37.hpp:261
Returns f(u(x)) where u is a scalar GridFunction and f:R → R.
Definition: ex37.hpp:40
Class for parallel linear form.
Definition: plinearform.hpp:26
Returns f(u(x)) - f(v(x)) where u, v are scalar GridFunctions and f:R → R.
Definition: ex37.hpp:65
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
virtual void SetAttributes()
Determine the sets of unique attribute values in domain and boundary elements.
Definition: mesh.cpp:1572
Class for solving linear elasticity:
Definition: ex37.hpp:303
double proj(ParGridFunction &psi, double target_volume, double tol=1e-12, int max_its=10)
Bregman projection of ρ = sigmoid(ψ) onto the subspace ∫_Ω ρ dx = θ vol(Ω) as follows: ...
Definition: ex37p.cpp:71
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Definition: plinearform.cpp:46
Strain energy density coefficient.
Definition: ex37.hpp:121
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void SetRHSCoefficient(Coefficient *rhscf_)
Definition: ex37.hpp:260
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:1125
double sigmoid(double x)
Sigmoid function.
Definition: ex37.hpp:20
void SetHighOrderOutput(bool high_order_output_)
void SetTime(double t)
Set physical time (for time-dependent simulations)
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
GridFunction * GetFEMSolution()
Definition: ex37.hpp:534
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
HYPRE_Int HYPRE_BigInt
void SetMesh(Mesh *mesh_)
Definition: ex37.hpp:248
void SetDiffusionCoefficient(Coefficient *diffcf_)
Definition: ex37.hpp:258
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Solid isotropic material penalization (SIMP) coefficient.
Definition: ex37.hpp:97
void SetOrder(int order_)
Definition: ex37.hpp:341
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:248
void mfem_warning(const char *msg)
Function called by the macro MFEM_WARNING.
Definition: error.cpp:187
void SetMassCoefficient(Coefficient *masscf_)
Definition: ex37.hpp:259
double mu
Definition: ex25.cpp:140
int dim
Definition: ex24.cpp:53
void SetOrder(int order_)
Definition: ex37.hpp:257
Class for solving Poisson&#39;s equation:
Definition: ex37.hpp:215
Adds new Domain Integrator. Assumes ownership of bfi.
void SetAttribute(const int attr)
Set element&#39;s attribute.
Definition: element.hpp:58
Class for parallel bilinear form.
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
void SetLevelsOfDetail(int levels_of_detail_)
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:678
const double alpha
Definition: ex15.cpp:369
Vector data type.
Definition: vector.hpp:58
void SetLameCoefficients(Coefficient *lambda_cf_, Coefficient *mu_cf_)
Definition: ex37.hpp:342
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
virtual void Save() override
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
const Element * GetBdrElement(int i) const
Return pointer to the i&#39;th boundary element object.
Definition: mesh.hpp:1158
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void SetMesh(Mesh *mesh_)
Definition: ex37.hpp:332
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Class for parallel meshes.
Definition: pmesh.hpp:32
Abstract data type element.
Definition: element.hpp:28
Volumetric force for linear elasticity.
Definition: ex37.hpp:167
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:1736
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
void SetRHSCoefficient(VectorCoefficient *rhs_cf_)
Definition: ex37.hpp:343
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
double f(const Vector &p)
int main(int argc, char *argv[])
Definition: ex37p.cpp:182