MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
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
55using namespace std;
56using 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 real_t Final volume, ∫_Ω sigmoid(ψ)
71 */
72real_t proj(GridFunction &psi, real_t target_volume, real_t 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());
79 int_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
80 LinearForm int_der_sigmoid_psi(psi.FESpace());
81 int_der_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(
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 real_t f = int_sigmoid_psi.Sum() - target_volume;
88
89 int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
90 const real_t df = int_der_sigmoid_psi.Sum();
91
92 const real_t 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
177int main(int argc, char *argv[])
178{
179 // 1. Parse command-line options.
180 int ref_levels = 5;
181 int order = 2;
182 real_t alpha = 1.0;
183 real_t epsilon = 0.01;
184 real_t vol_fraction = 0.5;
185 int max_it = 1e3;
186 real_t itol = 1e-1;
187 real_t ntol = 1e-4;
188 real_t rho_min = 1e-6;
189 real_t lambda = 1.0;
190 real_t mu = 1.0;
191 bool glvis_visualization = true;
192 bool paraview_output = false;
193
194 OptionsParser args(argc, argv);
195 args.AddOption(&ref_levels, "-r", "--refine",
196 "Number of times to refine the mesh uniformly.");
197 args.AddOption(&order, "-o", "--order",
198 "Order (degree) of the finite elements.");
199 args.AddOption(&alpha, "-alpha", "--alpha-step-length",
200 "Step length for gradient descent.");
201 args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
202 "Length scale for ρ.");
203 args.AddOption(&max_it, "-mi", "--max-it",
204 "Maximum number of gradient descent iterations.");
205 args.AddOption(&ntol, "-ntol", "--rel-tol",
206 "Normalized exit tolerance.");
207 args.AddOption(&itol, "-itol", "--abs-tol",
208 "Increment exit tolerance.");
209 args.AddOption(&vol_fraction, "-vf", "--volume-fraction",
210 "Volume fraction for the material density.");
211 args.AddOption(&lambda, "-lambda", "--lambda",
212 "Lamé constant λ.");
213 args.AddOption(&mu, "-mu", "--mu",
214 "Lamé constant μ.");
215 args.AddOption(&rho_min, "-rmin", "--psi-min",
216 "Minimum of density coefficient.");
217 args.AddOption(&glvis_visualization, "-vis", "--visualization", "-no-vis",
218 "--no-visualization",
219 "Enable or disable GLVis visualization.");
220 args.AddOption(&paraview_output, "-pv", "--paraview", "-no-pv",
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 }
230
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 real_t * coords1 = mesh.GetVertex(vertices[0]);
243 real_t * 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 real_t 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.
344 GridFunction grad(&control_fes);
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);
354 vol_form.AddDomainIntegrator(new DomainLFIntegrator(one));
355 vol_form.Assemble();
356 real_t domain_volume = vol_form(onegf);
357 const real_t 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 *= ((real_t) k) / ((real_t) 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();
421 M.Mult(w_rhs,grad);
422
423 // Step 5 - Update design variable ψ ← proj(ψ - αG)
424 psi.Add(-alpha, grad);
425 const real_t material_volume = proj(psi, target_volume);
426
427 // Compute ||ρ - ρ_old|| in control fes.
428 real_t norm_increment = zerogf.ComputeL1Error(succ_diff_rho);
429 real_t norm_reduced_gradient = norm_increment/alpha;
430 psi_old = psi;
431
432 real_t 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((real_t)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}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
@ GaussLobatto
Closed type.
Definition fe_base.hpp:32
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
A coefficient that is constant across space and time.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
Returns f(u(x)) - f(v(x)) where u, v are scalar GridFunctions and f:R → R.
Definition ex37.hpp:66
Class for solving Poisson's equation:
Definition ex37.hpp:216
void SetDiffusionCoefficient(Coefficient *diffcf_)
Definition ex37.hpp:258
void SetRHSCoefficient(Coefficient *rhscf_)
Definition ex37.hpp:260
void SetOrder(int order_)
Definition ex37.hpp:257
void SetMesh(Mesh *mesh_)
Definition ex37.hpp:248
void SetMassCoefficient(Coefficient *masscf_)
Definition ex37.hpp:259
GridFunction * GetFEMSolution()
Definition ex37.hpp:534
void SetEssentialBoundary(const Array< int > &ess_bdr_)
Definition ex37.hpp:261
Class for domain integration .
Definition lininteg.hpp:109
Abstract data type element.
Definition element.hpp:29
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
void SetAttribute(const int attr)
Set element's attribute.
Definition element.hpp:58
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition fe_coll.hpp:225
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual real_t ComputeL1Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition gridfunc.hpp:590
FiniteElementSpace * FESpace()
Definition gridfunc.hpp:696
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Integrator that inverts the matrix assembled by another integrator.
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Class for solving linear elasticity:
Definition ex37.hpp:304
GridFunction * GetFEMSolution()
Definition ex37.hpp:716
void SetEssentialBoundary(const Array< int > &ess_bdr_)
Definition ex37.hpp:344
void SetRHSCoefficient(VectorCoefficient *rhs_cf_)
Definition ex37.hpp:343
void SetMesh(Mesh *mesh_)
Definition ex37.hpp:332
LinearForm * GetLinearForm()
Definition ex37.hpp:353
void SetOrder(int order_)
Definition ex37.hpp:341
void SetLameCoefficients(Coefficient *lambda_cf_, Coefficient *mu_cf_)
Definition ex37.hpp:342
Vector with associated FE space and LinearFormIntegrators.
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Returns f(u(x)) where u is a scalar GridFunction and f:R → R.
Definition ex37.hpp:41
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition mesh.hpp:1298
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1229
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, real_t sx=1.0, real_t sy=1.0, bool sfc_ordering=true)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny quadrilaterals if type = QUADRILATER...
Definition mesh.cpp:4243
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
virtual void SetAttributes()
Determine the sets of unique attribute values in domain and boundary elements.
Definition mesh.cpp:1604
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Definition mesh.hpp:1265
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Helper class for ParaView visualization data.
void SetHighOrderOutput(bool high_order_output_)
void SetLevelsOfDetail(int levels_of_detail_)
virtual void Save() override
void SetDataFormat(VTKFormat fmt)
Scalar coefficient defined as the product of two scalar coefficients or a scalar and a scalar coeffic...
Solid isotropic material penalization (SIMP) coefficient.
Definition ex37.hpp:98
Data type sparse matrix.
Definition sparsemat.hpp:51
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Strain energy density coefficient.
Definition ex37.hpp:122
Vector data type.
Definition vector.hpp:80
real_t Sum() const
Return the sum of the vector entries.
Definition vector.cpp:1283
Vector & Add(const real_t a, const Vector &Va)
(*this) += a * Va
Definition vector.cpp:247
Volumetric force for linear elasticity.
Definition ex37.hpp:168
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
const real_t alpha
Definition ex15.cpp:369
int dim
Definition ex24.cpp:53
real_t mu
Definition ex25.cpp:140
real_t epsilon
Definition ex25.cpp:141
real_t proj(GridFunction &psi, real_t target_volume, real_t tol=1e-12, int max_its=10)
Bregman projection of ρ = sigmoid(ψ) onto the subspace ∫_Ω ρ dx = θ vol(Ω) as follows:
Definition ex37.cpp:72
int main()
const int visport
real_t sigmoid(real_t x)
Sigmoid function.
Definition ex37.hpp:20
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void mfem_warning(const char *msg)
Definition error.cpp:187
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
real_t inv_sigmoid(real_t x)
Inverse sigmoid function.
Definition ex37.hpp:12
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
real_t der_sigmoid(real_t x)
Derivative of sigmoid function.
Definition ex37.hpp:33
const char vishost[]