MFEM  v4.6.0 Finite element discretization library
ex37.cpp
Go to the documentation of this file.
1 // MFEM Example 37
2 //
3 // Compile with: make ex37
4 //
5 // Sample runs:
6 // ex37 -alpha 10
7 // ex37 -alpha 10 -pv
8 // ex37 -lambda 0.1 -mu 0.1
9 // ex37 -o 2 -alpha 5.0 -mi 50 -vf 0.4 -ntol 1e-5
10 // ex37 -r 6 -o 1 -alpha 25.0 -epsilon 0.02 -mi 50 -ntol 1e-5
11 //
12 // Description: This example code demonstrates the use of MFEM to solve a
13 // density-filtered [3] topology optimization problem. The
14 // objective is to minimize the compliance
15 //
16 // minimize ∫_Ω f⋅u dx over u ∈ [H¹(Ω)]² and ρ ∈ L¹(Ω)
17 //
18 // subject to
19 //
20 // -Div(r(ρ̃)Cε(u)) = f in Ω + BCs
21 // -ϵ²Δρ̃ + ρ̃ = ρ in Ω + Neumann BCs
22 // 0 ≤ ρ ≤ 1 in Ω
23 // ∫_Ω ρ dx = θ vol(Ω)
24 //
25 // Here, r(ρ̃) = ρ₀ + ρ̃³ (1-ρ₀) is the solid isotropic material
26 // penalization (SIMP) law, C is the elasticity tensor for an
27 // isotropic linearly elastic material, ϵ > 0 is the design
28 // length scale, and 0 < θ < 1 is the volume fraction.
29 //
30 // The problem is discretized and gradients are computing using
31 // finite elements [1]. The design is optimized using an entropic
32 // mirror descent algorithm introduced by Keith and Surowiec [2]
33 // that is tailored to the bound constraint 0 ≤ ρ ≤ 1.
34 //
35 // This example highlights the ability of MFEM to deliver high-
36 // order solutions to inverse design problems and showcases how
37 // to set up and solve PDE-constrained optimization problems
38 // using the so-called reduced space approach.
39 //
40 // [1] Andreassen, E., Clausen, A., Schevenels, M., Lazarov, B. S., & Sigmund, O.
41 // (2011). Efficient topology optimization in MATLAB using 88 lines of
42 // code. Structural and Multidisciplinary Optimization, 43(1), 1-16.
43 // [2] Keith, B. and Surowiec, T. (2023) Proximal Galerkin: A structure-
44 // preserving finite element method for pointwise bound constraints.
45 // arXiv:2307.12444 [math.NA]
46 // [3] Lazarov, B. S., & Sigmund, O. (2011). Filters in topology optimization
47 // based on Helmholtz‐type differential equations. International Journal
48 // for Numerical Methods in Engineering, 86(6), 765-781.
49
50 #include "mfem.hpp"
51 #include <iostream>
52 #include <fstream>
53 #include "ex37.hpp"
54
55 using namespace std;
56 using namespace mfem;
57
58 /**
59  * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
60  * ∫_Ω ρ dx = θ vol(Ω) as follows:
61  *
62  * 1. Compute the root of the R → R function
63  * f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
64  * 2. Set ψ ← ψ + c.
65  *
66  * @param psi a GridFunction to be updated
67  * @param target_volume θ vol(Ω)
68  * @param tol Newton iteration tolerance
69  * @param max_its Newton maximum iteration number
70  * @return double Final volume, ∫_Ω sigmoid(ψ)
71  */
72 double proj(GridFunction &psi, double target_volume, double tol=1e-12,
73  int max_its=10)
74 {
75  MappedGridFunctionCoefficient sigmoid_psi(&psi, sigmoid);
76  MappedGridFunctionCoefficient der_sigmoid_psi(&psi, der_sigmoid);
77
78  LinearForm int_sigmoid_psi(psi.FESpace());
80  LinearForm int_der_sigmoid_psi(psi.FESpace());
82  der_sigmoid_psi));
83  bool done = false;
84  for (int k=0; k<max_its; k++) // Newton iteration
85  {
86  int_sigmoid_psi.Assemble(); // Recompute f(c) with updated ψ
87  const double f = int_sigmoid_psi.Sum() - target_volume;
88
89  int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
90  const double df = int_der_sigmoid_psi.Sum();
91
92  const double dc = -f/df;
93  psi += dc;
94  if (abs(dc) < tol) { done = true; break; }
95  }
96  if (!done)
97  {
98  mfem_warning("Projection reached maximum iteration without converging. "
99  "Result may not be accurate.");
100  }
101  int_sigmoid_psi.Assemble();
102  return int_sigmoid_psi.Sum();
103 }
104
105 /**
106  * ---------------------------------------------------------------
107  * ALGORITHM PREAMBLE
108  * ---------------------------------------------------------------
109  *
110  * The Lagrangian for this problem is
111  *
112  * L(u,ρ,ρ̃,w,w̃) = (f,u) - (r(ρ̃) C ε(u),ε(w)) + (f,w)
113  * - (ϵ² ∇ρ̃,∇w̃) - (ρ̃,w̃) + (ρ,w̃)
114  *
115  * where
116  *
117  * r(ρ̃) = ρ₀ + ρ̃³ (1 - ρ₀) (SIMP rule)
118  *
119  * ε(u) = (∇u + ∇uᵀ)/2 (symmetric gradient)
120  *
121  * C e = λtr(e)I + 2μe (isotropic material)
122  *
123  * NOTE: The Lame parameters can be computed from Young's modulus E
124  * and Poisson's ratio ν as follows:
125  *
126  * λ = E ν/((1+ν)(1-2ν)), μ = E/(2(1+ν))
127  *
128  * ---------------------------------------------------------------
129  *
130  * Discretization choices:
131  *
132  * u ∈ V ⊂ (H¹)ᵈ (order p)
133  * ψ ∈ L² (order p - 1), ρ = sigmoid(ψ)
134  * ρ̃ ∈ H¹ (order p)
135  * w ∈ V (order p)
136  * w̃ ∈ H¹ (order p)
137  *
138  * ---------------------------------------------------------------
139  * ALGORITHM
140  * ---------------------------------------------------------------
141  *
142  * Update ρ with projected mirror descent via the following algorithm.
143  *
144  * 1. Initialize ψ = inv_sigmoid(vol_fraction) so that ∫ sigmoid(ψ) = θ vol(Ω)
145  *
146  * While not converged:
147  *
148  * 2. Solve filter equation ∂_w̃ L = 0; i.e.,
149  *
150  * (ϵ² ∇ ρ̃, ∇ v ) + (ρ̃,v) = (ρ,v) ∀ v ∈ H¹.
151  *
152  * 3. Solve primal problem ∂_w L = 0; i.e.,
153  *
154  * (λ r(ρ̃) ∇⋅u, ∇⋅v) + (2 μ r(ρ̃) ε(u), ε(v)) = (f,v) ∀ v ∈ V.
155  *
156  * NB. The dual problem ∂_u L = 0 is the negative of the primal problem due to symmetry.
157  *
158  * 4. Solve for filtered gradient ∂_ρ̃ L = 0; i.e.,
159  *
160  * (ϵ² ∇ w̃ , ∇ v ) + (w̃ ,v) = (-r'(ρ̃) ( λ |∇⋅u|² + 2 μ |ε(u)|²),v) ∀ v ∈ H¹.
161  *
162  * 5. Project the gradient onto the discrete latent space; i.e., solve
163  *
164  * (G,v) = (w̃,v) ∀ v ∈ L².
165  *
166  * 6. Bregman proximal gradient update; i.e.,
167  *
168  * ψ ← ψ - αG + c,
169  *
170  * where α > 0 is a step size parameter and c ∈ R is a constant ensuring
171  *
172  * ∫_Ω sigmoid(ψ - αG + c) dx = θ vol(Ω).
173  *
174  * end
175  */
176
177 int main(int argc, char *argv[])
178 {
179  // 1. Parse command-line options.
180  int ref_levels = 5;
181  int order = 2;
182  double alpha = 1.0;
183  double epsilon = 0.01;
184  double vol_fraction = 0.5;
185  int max_it = 1e3;
186  double itol = 1e-1;
187  double ntol = 1e-4;
188  double rho_min = 1e-6;
189  double lambda = 1.0;
190  double mu = 1.0;
191  bool glvis_visualization = true;
192  bool paraview_output = false;
193
194  OptionsParser args(argc, argv);
196  "Number of times to refine the mesh uniformly.");
198  "Order (degree) of the finite elements.");
200  "Step length for gradient descent.");
202  "Length scale for ρ.");
204  "Maximum number of gradient descent iterations.");
206  "Normalized exit tolerance.");
208  "Increment exit tolerance.");
210  "Volume fraction for the material density.");
212  "Lamé constant λ.");
214  "Lamé constant μ.");
216  "Minimum of density coefficient.");
218  "--no-visualization",
219  "Enable or disable GLVis visualization.");
221  "--no-paraview",
222  "Enable or disable ParaView output.");
223  args.Parse();
224  if (!args.Good())
225  {
226  args.PrintUsage(mfem::out);
227  return 1;
228  }
229  args.PrintOptions(mfem::out);
230
231  Mesh mesh = Mesh::MakeCartesian2D(3, 1, mfem::Element::Type::QUADRILATERAL,
232  true, 3.0, 1.0);
233  int dim = mesh.Dimension();
234
235  // 2. Set BCs.
236  for (int i = 0; i<mesh.GetNBE(); i++)
237  {
238  Element * be = mesh.GetBdrElement(i);
239  Array<int> vertices;
240  be->GetVertices(vertices);
241
242  double * coords1 = mesh.GetVertex(vertices[0]);
243  double * coords2 = mesh.GetVertex(vertices[1]);
244
245  Vector center(2);
246  center(0) = 0.5*(coords1[0] + coords2[0]);
247  center(1) = 0.5*(coords1[1] + coords2[1]);
248
249  if (abs(center(0) - 0.0) < 1e-10)
250  {
251  // the left edge
252  be->SetAttribute(1);
253  }
254  else
255  {
256  // all other boundaries
257  be->SetAttribute(2);
258  }
259  }
260  mesh.SetAttributes();
261
262  // 3. Refine the mesh.
263  for (int lev = 0; lev < ref_levels; lev++)
264  {
265  mesh.UniformRefinement();
266  }
267
268  // 4. Define the necessary finite element spaces on the mesh.
269  H1_FECollection state_fec(order, dim); // space for u
270  H1_FECollection filter_fec(order, dim); // space for ρ̃
271  L2_FECollection control_fec(order-1, dim,
272  BasisType::GaussLobatto); // space for ψ
273  FiniteElementSpace state_fes(&mesh, &state_fec,dim);
274  FiniteElementSpace filter_fes(&mesh, &filter_fec);
275  FiniteElementSpace control_fes(&mesh, &control_fec);
276
277  int state_size = state_fes.GetTrueVSize();
278  int control_size = control_fes.GetTrueVSize();
279  int filter_size = filter_fes.GetTrueVSize();
280  mfem::out << "Number of state unknowns: " << state_size << std::endl;
281  mfem::out << "Number of filter unknowns: " << filter_size << std::endl;
282  mfem::out << "Number of control unknowns: " << control_size << std::endl;
283
284  // 5. Set the initial guess for ρ.
285  GridFunction u(&state_fes);
286  GridFunction psi(&control_fes);
287  GridFunction psi_old(&control_fes);
288  GridFunction rho_filter(&filter_fes);
289  u = 0.0;
290  rho_filter = vol_fraction;
291  psi = inv_sigmoid(vol_fraction);
292  psi_old = inv_sigmoid(vol_fraction);
293
294  // ρ = sigmoid(ψ)
296  // Interpolation of ρ = sigmoid(ψ) in control fes (for ParaView output)
297  GridFunction rho_gf(&control_fes);
298  // ρ - ρ_old = sigmoid(ψ) - sigmoid(ψ_old)
299  DiffMappedGridFunctionCoefficient succ_diff_rho(&psi, &psi_old, sigmoid);
300
301  // 6. Set-up the physics solver.
302  int maxat = mesh.bdr_attributes.Max();
303  Array<int> ess_bdr(maxat);
304  ess_bdr = 0;
305  ess_bdr[0] = 1;
306  ConstantCoefficient one(1.0);
307  ConstantCoefficient lambda_cf(lambda);
308  ConstantCoefficient mu_cf(mu);
309  LinearElasticitySolver * ElasticitySolver = new LinearElasticitySolver();
310  ElasticitySolver->SetMesh(&mesh);
311  ElasticitySolver->SetOrder(state_fec.GetOrder());
312  ElasticitySolver->SetupFEM();
313  Vector center(2); center(0) = 2.9; center(1) = 0.5;
314  Vector force(2); force(0) = 0.0; force(1) = -1.0;
315  double r = 0.05;
316  VolumeForceCoefficient vforce_cf(r,center,force);
317  ElasticitySolver->SetRHSCoefficient(&vforce_cf);
318  ElasticitySolver->SetEssentialBoundary(ess_bdr);
319
320  // 7. Set-up the filter solver.
322  DiffusionSolver * FilterSolver = new DiffusionSolver();
323  FilterSolver->SetMesh(&mesh);
324  FilterSolver->SetOrder(filter_fec.GetOrder());
325  FilterSolver->SetDiffusionCoefficient(&eps2_cf);
326  FilterSolver->SetMassCoefficient(&one);
327  Array<int> ess_bdr_filter;
328  if (mesh.bdr_attributes.Size())
329  {
330  ess_bdr_filter.SetSize(mesh.bdr_attributes.Max());
331  ess_bdr_filter = 0;
332  }
333  FilterSolver->SetEssentialBoundary(ess_bdr_filter);
334  FilterSolver->SetupFEM();
335
336  BilinearForm mass(&control_fes);
338  mass.Assemble();
339  SparseMatrix M;
340  Array<int> empty;
341  mass.FormSystemMatrix(empty,M);
342
343  // 8. Define the Lagrange multiplier and gradient functions.
345  GridFunction w_filter(&filter_fes);
346
347  // 9. Define some tools for later.
348  ConstantCoefficient zero(0.0);
349  GridFunction onegf(&control_fes);
350  onegf = 1.0;
351  GridFunction zerogf(&control_fes);
352  zerogf = 0.0;
353  LinearForm vol_form(&control_fes);
355  vol_form.Assemble();
356  double domain_volume = vol_form(onegf);
357  const double target_volume = domain_volume * vol_fraction;
358
359  // 10. Connect to GLVis. Prepare for VisIt output.
360  char vishost[] = "localhost";
361  int visport = 19916;
362  socketstream sout_r;
363  if (glvis_visualization)
364  {
365  sout_r.open(vishost, visport);
366  sout_r.precision(8);
367  }
368
369  mfem::ParaViewDataCollection paraview_dc("ex37", &mesh);
370  if (paraview_output)
371  {
372  rho_gf.ProjectCoefficient(rho);
373  paraview_dc.SetPrefixPath("ParaView");
374  paraview_dc.SetLevelsOfDetail(order);
375  paraview_dc.SetDataFormat(VTKFormat::BINARY);
376  paraview_dc.SetHighOrderOutput(true);
377  paraview_dc.SetCycle(0);
378  paraview_dc.SetTime(0.0);
379  paraview_dc.RegisterField("displacement",&u);
380  paraview_dc.RegisterField("density",&rho_gf);
381  paraview_dc.RegisterField("filtered_density",&rho_filter);
382  paraview_dc.Save();
383  }
384
385  // 11. Iterate:
386  for (int k = 1; k <= max_it; k++)
387  {
388  if (k > 1) { alpha *= ((double) k) / ((double) k-1); }
389
390  mfem::out << "\nStep = " << k << std::endl;
391
392  // Step 1 - Filter solve
393  // Solve (ϵ^2 ∇ ρ̃, ∇ v ) + (ρ̃,v) = (ρ,v)
394  FilterSolver->SetRHSCoefficient(&rho);
395  FilterSolver->Solve();
396  rho_filter = *FilterSolver->GetFEMSolution();
397
398  // Step 2 - State solve
399  // Solve (λ r(ρ̃) ∇⋅u, ∇⋅v) + (2 μ r(ρ̃) ε(u), ε(v)) = (f,v)
400  SIMPInterpolationCoefficient SIMP_cf(&rho_filter,rho_min, 1.0);
401  ProductCoefficient lambda_SIMP_cf(lambda_cf,SIMP_cf);
402  ProductCoefficient mu_SIMP_cf(mu_cf,SIMP_cf);
403  ElasticitySolver->SetLameCoefficients(&lambda_SIMP_cf,&mu_SIMP_cf);
404  ElasticitySolver->Solve();
405  u = *ElasticitySolver->GetFEMSolution();
406
407  // Step 3 - Adjoint filter solve
408  // Solve (ϵ² ∇ w̃, ∇ v) + (w̃ ,v) = (-r'(ρ̃) ( λ |∇⋅u|² + 2 μ |ε(u)|²),v)
409  StrainEnergyDensityCoefficient rhs_cf(&lambda_cf,&mu_cf,&u, &rho_filter,
410  rho_min);
411  FilterSolver->SetRHSCoefficient(&rhs_cf);
412  FilterSolver->Solve();
413  w_filter = *FilterSolver->GetFEMSolution();
414
415  // Step 4 - Compute gradient
416  // Solve G = M⁻¹w̃
417  GridFunctionCoefficient w_cf(&w_filter);
418  LinearForm w_rhs(&control_fes);
420  w_rhs.Assemble();
422
423  // Step 5 - Update design variable ψ ← proj(ψ - αG)
425  const double material_volume = proj(psi, target_volume);
426
427  // Compute ||ρ - ρ_old|| in control fes.
428  double norm_increment = zerogf.ComputeL1Error(succ_diff_rho);
430  psi_old = psi;
431
432  double compliance = (*(ElasticitySolver->GetLinearForm()))(u);
433  mfem::out << "norm of the reduced gradient = " << norm_reduced_gradient <<
434  std::endl;
435  mfem::out << "norm of the increment = " << norm_increment << endl;
436  mfem::out << "compliance = " << compliance << std::endl;
437  mfem::out << "volume fraction = " << material_volume / domain_volume <<
438  std::endl;
439
440  if (glvis_visualization)
441  {
442  GridFunction r_gf(&filter_fes);
443  r_gf.ProjectCoefficient(SIMP_cf);
444  sout_r << "solution\n" << mesh << r_gf
445  << "window_title 'Design density r(ρ̃)'" << flush;
446  }
447
448  if (paraview_output)
449  {
450  rho_gf.ProjectCoefficient(rho);
451  paraview_dc.SetCycle(k);
452  paraview_dc.SetTime((double)k);
453  paraview_dc.Save();
454  }
455
456  if (norm_reduced_gradient < ntol && norm_increment < itol)
457  {
458  break;
459  }
460  }
461
462  delete ElasticitySolver;
463  delete FilterSolver;
464
465  return 0;
466 }
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:733
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
int visport
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
GridFunction * GetFEMSolution()
Definition: ex37.hpp:716
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
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
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:168
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
STL namespace.
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
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[]
Data type sparse matrix.
Definition: sparsemat.hpp:50
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
Strain energy density coefficient.
Definition: ex37.hpp:121
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
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
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
int main(int argc, char *argv[])
Definition: ex37.cpp:177
void SetHighOrderOutput(bool high_order_output_)
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
void SetTime(double t)
Set physical time (for time-dependent simulations)
virtual double ComputeL1Error(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:589
GridFunction * GetFEMSolution()
Definition: ex37.hpp:534
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
double proj(GridFunction &psi, double target_volume, double tol=1e-12, int max_its=10)
Bregman projection of ρ = sigmoid(ψ) onto the subspace ∫_Ω ρ dx = θ vol(Ω) as follows: ...
Definition: ex37.cpp:72
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
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
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
void SetMesh(Mesh *mesh_)
Definition: ex37.hpp:248
void SetDiffusionCoefficient(Coefficient *diffcf_)
Definition: ex37.hpp:258
Solid isotropic material penalization (SIMP) coefficient.
Definition: ex37.hpp:97
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
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
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2415
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_)
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
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
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
void SetMesh(Mesh *mesh_)
Definition: ex37.hpp:332
Abstract data type element.
Definition: element.hpp:28
Volumetric force for linear elasticity.
Definition: ex37.hpp:167
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)