MFEM v4.7.0 Finite element discretization library
Loading...
Searching...
No Matches
ex37.hpp
Go to the documentation of this file.
1// MFEM Example 37 - Serial/Parallel Shared Code
2
3#include "mfem.hpp"
4#include <fstream>
5#include <iostream>
6#include <functional>
7
8namespace mfem
9{
10
11/// @brief Inverse sigmoid function
13{
14 real_t tol = 1e-12;
15 x = std::min(std::max(tol,x), real_t(1.0)-tol);
16 return std::log(x/(1.0-x));
17}
18
19/// @brief Sigmoid function
21{
22 if (x >= 0)
23 {
24 return 1.0/(1.0+std::exp(-x));
25 }
26 else
27 {
28 return std::exp(x)/(1.0+std::exp(x));
29 }
30}
31
32/// @brief Derivative of sigmoid function
34{
35 real_t tmp = sigmoid(-x);
36 return tmp - std::pow(tmp,2);
37}
38
39/// @brief Returns f(u(x)) where u is a scalar GridFunction and f:R → R
41{
42protected:
43 std::function<real_t(const real_t)> fun; // f:R → R
44public:
47 fun([](real_t x) {return x;}) {}
49 std::function<real_t(const real_t)> fun_,
50 int comp=1)
51 :GridFunctionCoefficient(gf, comp),
52 fun(fun_) {}
53
54
56 const IntegrationPoint &ip)
57 {
59 }
60 void SetFunction(std::function<real_t(const real_t)> fun_) { fun = fun_; }
61};
62
63
64/// @brief Returns f(u(x)) - f(v(x)) where u, v are scalar GridFunctions and f:R → R
66{
67protected:
70 std::function<real_t(const real_t)> fun; // f:R → R
71public:
74 OtherGridF(nullptr),
76 fun([](real_t x) {return x;}) {}
78 const GridFunction *other_gf,
79 std::function<real_t(const real_t)> fun_,
80 int comp=1)
81 :GridFunctionCoefficient(gf, comp),
82 OtherGridF(other_gf),
84 fun(fun_) {}
85
87 const IntegrationPoint &ip)
88 {
89 const real_t value1 = fun(GridFunctionCoefficient::Eval(T, ip));
90 const real_t value2 = fun(OtherGridF_cf.Eval(T, ip));
91 return value1 - value2;
92 }
93 void SetFunction(std::function<real_t(const real_t)> fun_) { fun = fun_; }
94};
95
96/// @brief Solid isotropic material penalization (SIMP) coefficient
98{
99protected:
104
105public:
107 real_t max_val_ = 1.0, real_t exponent_ = 3)
108 : rho_filter(rho_filter_), min_val(min_val_), max_val(max_val_),
109 exponent(exponent_) { }
110
112 {
113 real_t val = rho_filter->GetValue(T, ip);
114 real_t coeff = min_val + pow(val,exponent)*(max_val-min_val);
115 return coeff;
116 }
117};
118
119
120/// @brief Strain energy density coefficient
122{
123protected:
125 Coefficient * mu=nullptr;
126 GridFunction *u = nullptr; // displacement
127 GridFunction *rho_filter = nullptr; // filter density
128 DenseMatrix grad; // auxiliary matrix, used in Eval
131
132public:
134 GridFunction * u_, GridFunction * rho_filter_, real_t rho_min_=1e-6,
135 real_t exponent_ = 3.0)
136 : lambda(lambda_), mu(mu_), u(u_), rho_filter(rho_filter_),
137 exponent(exponent_), rho_min(rho_min_)
138 {
139 MFEM_ASSERT(rho_min_ >= 0.0, "rho_min must be >= 0");
140 MFEM_ASSERT(rho_min_ < 1.0, "rho_min must be > 1");
141 MFEM_ASSERT(u, "displacement field is not set");
142 MFEM_ASSERT(rho_filter, "density field is not set");
143 }
144
145 virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)
146 {
147 real_t L = lambda->Eval(T, ip);
148 real_t M = mu->Eval(T, ip);
149 u->GetVectorGradient(T, grad);
150 real_t div_u = grad.Trace();
151 real_t density = L*div_u*div_u;
152 int dim = T.GetSpaceDim();
153 for (int i=0; i<dim; i++)
154 {
155 for (int j=0; j<dim; j++)
156 {
157 density += M*grad(i,j)*(grad(i,j)+grad(j,i));
158 }
159 }
160 real_t val = rho_filter->GetValue(T,ip);
161
162 return -exponent * pow(val, exponent-1.0) * (1-rho_min) * density;
163 }
164};
165
166/// @brief Volumetric force for linear elasticity
167class VolumeForceCoefficient : public VectorCoefficient
168{
169private:
170 real_t r;
171 Vector center;
172 Vector force;
173public:
174 VolumeForceCoefficient(real_t r_,Vector & center_, Vector & force_) :
175 VectorCoefficient(center_.Size()), r(r_), center(center_), force(force_) { }
176
177 using VectorCoefficient::Eval;
178
179 virtual void Eval(Vector &V, ElementTransformation &T,
180 const IntegrationPoint &ip)
181 {
182 Vector xx; xx.SetSize(T.GetDimension());
183 T.Transform(ip,xx);
184 for (int i=0; i<xx.Size(); i++)
185 {
186 xx[i]=xx[i]-center[i];
187 }
188
189 real_t cr=xx.Norml2();
190 V.SetSize(T.GetDimension());
191 if (cr <= r)
192 {
193 V = force;
194 }
195 else
196 {
197 V = 0.0;
198 }
199 }
200
201 void Set(real_t r_,Vector & center_, Vector & force_)
202 {
203 r=r_;
204 center = center_;
205 force = force_;
206 }
207};
208
209/**
210 * @brief Class for solving Poisson's equation:
211 *
212 * - ∇ ⋅(κ ∇ u) = f in Ω
213 *
214 */
215class DiffusionSolver
216{
217private:
218 Mesh * mesh = nullptr;
219 int order = 1;
220 // diffusion coefficient
221 Coefficient * diffcf = nullptr;
222 // mass coefficient
223 Coefficient * masscf = nullptr;
224 Coefficient * rhscf = nullptr;
225 Coefficient * essbdr_cf = nullptr;
226 Coefficient * neumann_cf = nullptr;
227 VectorCoefficient * gradient_cf = nullptr;
228
229 // FEM solver
230 int dim;
231 FiniteElementCollection * fec = nullptr;
232 FiniteElementSpace * fes = nullptr;
233 Array<int> ess_bdr;
234 Array<int> neumann_bdr;
235 GridFunction * u = nullptr;
236 LinearForm * b = nullptr;
237 bool parallel;
238#ifdef MFEM_USE_MPI
239 ParMesh * pmesh = nullptr;
240 ParFiniteElementSpace * pfes = nullptr;
241#endif
242
243public:
244 DiffusionSolver() { }
245 DiffusionSolver(Mesh * mesh_, int order_, Coefficient * diffcf_,
246 Coefficient * cf_);
247
248 void SetMesh(Mesh * mesh_)
249 {
250 mesh = mesh_;
251 parallel = false;
252#ifdef MFEM_USE_MPI
253 pmesh = dynamic_cast<ParMesh *>(mesh);
254 if (pmesh) { parallel = true; }
255#endif
256 }
257 void SetOrder(int order_) { order = order_ ; }
258 void SetDiffusionCoefficient(Coefficient * diffcf_) { diffcf = diffcf_; }
259 void SetMassCoefficient(Coefficient * masscf_) { masscf = masscf_; }
260 void SetRHSCoefficient(Coefficient * rhscf_) { rhscf = rhscf_; }
261 void SetEssentialBoundary(const Array<int> & ess_bdr_) { ess_bdr = ess_bdr_;};
262 void SetNeumannBoundary(const Array<int> & neumann_bdr_) { neumann_bdr = neumann_bdr_;};
263 void SetNeumannData(Coefficient * neumann_cf_) {neumann_cf = neumann_cf_;}
264 void SetEssBdrData(Coefficient * essbdr_cf_) {essbdr_cf = essbdr_cf_;}
265 void SetGradientData(VectorCoefficient * gradient_cf_) {gradient_cf = gradient_cf_;}
266
267 void ResetFEM();
268 void SetupFEM();
269
270 void Solve();
271 GridFunction * GetFEMSolution();
272 LinearForm * GetLinearForm() {return b;}
273#ifdef MFEM_USE_MPI
274 ParGridFunction * GetParFEMSolution();
275 ParLinearForm * GetParLinearForm()
276 {
277 if (parallel)
278 {
279 return dynamic_cast<ParLinearForm *>(b);
280 }
281 else
282 {
283 MFEM_ABORT("Wrong code path. Call GetLinearForm");
284 return nullptr;
285 }
286 }
287#endif
288
289 ~DiffusionSolver();
290
291};
292
293/**
294 * @brief Class for solving linear elasticity:
295 *
296 * -∇ ⋅ σ(u) = f in Ω + BCs
297 *
298 * where
299 *
300 * σ(u) = λ ∇⋅u I + μ (∇ u + ∇uᵀ)
301 *
302 */
303class LinearElasticitySolver
304{
305private:
306 Mesh * mesh = nullptr;
307 int order = 1;
308 Coefficient * lambda_cf = nullptr;
309 Coefficient * mu_cf = nullptr;
310 VectorCoefficient * essbdr_cf = nullptr;
311 VectorCoefficient * rhs_cf = nullptr;
312
313 // FEM solver
314 int dim;
315 FiniteElementCollection * fec = nullptr;
316 FiniteElementSpace * fes = nullptr;
317 Array<int> ess_bdr;
318 Array<int> neumann_bdr;
319 GridFunction * u = nullptr;
320 LinearForm * b = nullptr;
321 bool parallel;
322#ifdef MFEM_USE_MPI
323 ParMesh * pmesh = nullptr;
324 ParFiniteElementSpace * pfes = nullptr;
325#endif
326
327public:
328 LinearElasticitySolver() { }
329 LinearElasticitySolver(Mesh * mesh_, int order_,
330 Coefficient * lambda_cf_, Coefficient * mu_cf_);
331
332 void SetMesh(Mesh * mesh_)
333 {
334 mesh = mesh_;
335 parallel = false;
336#ifdef MFEM_USE_MPI
337 pmesh = dynamic_cast<ParMesh *>(mesh);
338 if (pmesh) { parallel = true; }
339#endif
340 }
341 void SetOrder(int order_) { order = order_ ; }
342 void SetLameCoefficients(Coefficient * lambda_cf_, Coefficient * mu_cf_) { lambda_cf = lambda_cf_; mu_cf = mu_cf_; }
343 void SetRHSCoefficient(VectorCoefficient * rhs_cf_) { rhs_cf = rhs_cf_; }
344 void SetEssentialBoundary(const Array<int> & ess_bdr_) { ess_bdr = ess_bdr_;};
345 void SetNeumannBoundary(const Array<int> & neumann_bdr_) { neumann_bdr = neumann_bdr_;};
346 void SetEssBdrData(VectorCoefficient * essbdr_cf_) {essbdr_cf = essbdr_cf_;}
347
348 void ResetFEM();
349 void SetupFEM();
350
351 void Solve();
352 GridFunction * GetFEMSolution();
353 LinearForm * GetLinearForm() {return b;}
354#ifdef MFEM_USE_MPI
355 ParGridFunction * GetParFEMSolution();
356 ParLinearForm * GetParLinearForm()
357 {
358 if (parallel)
359 {
360 return dynamic_cast<ParLinearForm *>(b);
361 }
362 else
363 {
364 MFEM_ABORT("Wrong code path. Call GetLinearForm");
365 return nullptr;
366 }
367 }
368#endif
369
370 ~LinearElasticitySolver();
371
372};
373
374
375// Poisson solver
376
377DiffusionSolver::DiffusionSolver(Mesh * mesh_, int order_,
378 Coefficient * diffcf_, Coefficient * rhscf_)
379 : mesh(mesh_), order(order_), diffcf(diffcf_), rhscf(rhscf_)
380{
381
382#ifdef MFEM_USE_MPI
383 pmesh = dynamic_cast<ParMesh *>(mesh);
384 if (pmesh) { parallel = true; }
385#endif
386
387 SetupFEM();
388}
389
390void DiffusionSolver::SetupFEM()
391{
392 dim = mesh->Dimension();
393 fec = new H1_FECollection(order, dim);
394
395#ifdef MFEM_USE_MPI
396 if (parallel)
397 {
398 pfes = new ParFiniteElementSpace(pmesh, fec);
399 u = new ParGridFunction(pfes);
400 b = new ParLinearForm(pfes);
401 }
402 else
403 {
404 fes = new FiniteElementSpace(mesh, fec);
405 u = new GridFunction(fes);
406 b = new LinearForm(fes);
407 }
408#else
409 fes = new FiniteElementSpace(mesh, fec);
410 u = new GridFunction(fes);
411 b = new LinearForm(fes);
412#endif
413 *u=0.0;
414
415 if (!ess_bdr.Size())
416 {
417 if (mesh->bdr_attributes.Size())
418 {
419 ess_bdr.SetSize(mesh->bdr_attributes.Max());
420 ess_bdr = 1;
421 }
422 }
423}
424
425void DiffusionSolver::Solve()
426{
427 OperatorPtr A;
428 Vector B, X;
429 Array<int> ess_tdof_list;
430
431#ifdef MFEM_USE_MPI
432 if (parallel)
433 {
434 pfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
435 }
436 else
437 {
438 fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
439 }
440#else
441 fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
442#endif
443 *u=0.0;
444 if (b)
445 {
446 delete b;
447#ifdef MFEM_USE_MPI
448 if (parallel)
449 {
450 b = new ParLinearForm(pfes);
451 }
452 else
453 {
454 b = new LinearForm(fes);
455 }
456#else
457 b = new LinearForm(fes);
458#endif
459 }
460 if (rhscf)
461 {
462 b->AddDomainIntegrator(new DomainLFIntegrator(*rhscf));
463 }
464 if (neumann_cf)
465 {
466 MFEM_VERIFY(neumann_bdr.Size(), "neumann_bdr attributes not provided");
467 b->AddBoundaryIntegrator(new BoundaryLFIntegrator(*neumann_cf),neumann_bdr);
468 }
469 else if (gradient_cf)
470 {
471 MFEM_VERIFY(neumann_bdr.Size(), "neumann_bdr attributes not provided");
472 b->AddBoundaryIntegrator(new BoundaryNormalLFIntegrator(*gradient_cf),
473 neumann_bdr);
474 }
475
476 b->Assemble();
477
478 BilinearForm * a = nullptr;
479
480#ifdef MFEM_USE_MPI
481 if (parallel)
482 {
483 a = new ParBilinearForm(pfes);
484 }
485 else
486 {
487 a = new BilinearForm(fes);
488 }
489#else
490 a = new BilinearForm(fes);
491#endif
492 a->AddDomainIntegrator(new DiffusionIntegrator(*diffcf));
493 if (masscf)
494 {
495 a->AddDomainIntegrator(new MassIntegrator(*masscf));
496 }
497 a->Assemble();
498 if (essbdr_cf)
499 {
500 u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
501 }
502 a->FormLinearSystem(ess_tdof_list, *u, *b, A, X, B);
503
504 CGSolver * cg = nullptr;
505 Solver * M = nullptr;
506#ifdef MFEM_USE_MPI
507 if (parallel)
508 {
509 M = new HypreBoomerAMG;
510 dynamic_cast<HypreBoomerAMG*>(M)->SetPrintLevel(0);
511 cg = new CGSolver(pmesh->GetComm());
512 }
513 else
514 {
515 M = new GSSmoother((SparseMatrix&)(*A));
516 cg = new CGSolver;
517 }
518#else
519 M = new GSSmoother((SparseMatrix&)(*A));
520 cg = new CGSolver;
521#endif
522 cg->SetRelTol(1e-12);
523 cg->SetMaxIter(10000);
524 cg->SetPrintLevel(0);
525 cg->SetPreconditioner(*M);
526 cg->SetOperator(*A);
527 cg->Mult(B, X);
528 delete M;
529 delete cg;
530 a->RecoverFEMSolution(X, *b, *u);
531 delete a;
532}
533
534GridFunction * DiffusionSolver::GetFEMSolution()
535{
536 return u;
537}
538
539#ifdef MFEM_USE_MPI
540ParGridFunction * DiffusionSolver::GetParFEMSolution()
541{
542 if (parallel)
543 {
544 return dynamic_cast<ParGridFunction*>(u);
545 }
546 else
547 {
548 MFEM_ABORT("Wrong code path. Call GetFEMSolution");
549 return nullptr;
550 }
551}
552#endif
553
554DiffusionSolver::~DiffusionSolver()
555{
556 delete u; u = nullptr;
557 delete fes; fes = nullptr;
558#ifdef MFEM_USE_MPI
559 delete pfes; pfes=nullptr;
560#endif
561 delete fec; fec = nullptr;
562 delete b;
563}
564
565
566// Elasticity solver
567
568LinearElasticitySolver::LinearElasticitySolver(Mesh * mesh_, int order_,
569 Coefficient * lambda_cf_, Coefficient * mu_cf_)
570 : mesh(mesh_), order(order_), lambda_cf(lambda_cf_), mu_cf(mu_cf_)
571{
572#ifdef MFEM_USE_MPI
573 pmesh = dynamic_cast<ParMesh *>(mesh);
574 if (pmesh) { parallel = true; }
575#endif
576 SetupFEM();
577}
578
579void LinearElasticitySolver::SetupFEM()
580{
581 dim = mesh->Dimension();
582 fec = new H1_FECollection(order, dim,BasisType::Positive);
583
584#ifdef MFEM_USE_MPI
585 if (parallel)
586 {
587 pfes = new ParFiniteElementSpace(pmesh, fec, dim);
588 u = new ParGridFunction(pfes);
589 b = new ParLinearForm(pfes);
590 }
591 else
592 {
593 fes = new FiniteElementSpace(mesh, fec,dim);
594 u = new GridFunction(fes);
595 b = new LinearForm(fes);
596 }
597#else
598 fes = new FiniteElementSpace(mesh, fec, dim);
599 u = new GridFunction(fes);
600 b = new LinearForm(fes);
601#endif
602 *u=0.0;
603
604 if (!ess_bdr.Size())
605 {
606 if (mesh->bdr_attributes.Size())
607 {
608 ess_bdr.SetSize(mesh->bdr_attributes.Max());
609 ess_bdr = 1;
610 }
611 }
612}
613
614void LinearElasticitySolver::Solve()
615{
616 GridFunction * x = nullptr;
617 OperatorPtr A;
618 Vector B, X;
619 Array<int> ess_tdof_list;
620
621#ifdef MFEM_USE_MPI
622 if (parallel)
623 {
624 x = new ParGridFunction(pfes);
625 pfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
626 }
627 else
628 {
629 x = new GridFunction(fes);
630 fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
631 }
632#else
633 x = new GridFunction(fes);
634 fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
635#endif
636 *u=0.0;
637 if (b)
638 {
639 delete b;
640#ifdef MFEM_USE_MPI
641 if (parallel)
642 {
643 b = new ParLinearForm(pfes);
644 }
645 else
646 {
647 b = new LinearForm(fes);
648 }
649#else
650 b = new LinearForm(fes);
651#endif
652 }
653 if (rhs_cf)
654 {
655 b->AddDomainIntegrator(new VectorDomainLFIntegrator(*rhs_cf));
656 }
657
658 b->Assemble();
659
660 *x = 0.0;
661
662 BilinearForm * a = nullptr;
663
664#ifdef MFEM_USE_MPI
665 if (parallel)
666 {
667 a = new ParBilinearForm(pfes);
668 }
669 else
670 {
671 a = new BilinearForm(fes);
672 }
673#else
674 a = new BilinearForm(fes);
675#endif
676 a->AddDomainIntegrator(new ElasticityIntegrator(*lambda_cf, *mu_cf));
677 a->Assemble();
678 if (essbdr_cf)
679 {
680 u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
681 }
682 a->FormLinearSystem(ess_tdof_list, *x, *b, A, X, B);
683
684 CGSolver * cg = nullptr;
685 Solver * M = nullptr;
686#ifdef MFEM_USE_MPI
687 if (parallel)
688 {
689 M = new HypreBoomerAMG;
690 dynamic_cast<HypreBoomerAMG*>(M)->SetPrintLevel(0);
691 cg = new CGSolver(pmesh->GetComm());
692 }
693 else
694 {
695 M = new GSSmoother((SparseMatrix&)(*A));
696 cg = new CGSolver;
697 }
698#else
699 M = new GSSmoother((SparseMatrix&)(*A));
700 cg = new CGSolver;
701#endif
702 cg->SetRelTol(1e-10);
703 cg->SetMaxIter(10000);
704 cg->SetPrintLevel(0);
705 cg->SetPreconditioner(*M);
706 cg->SetOperator(*A);
707 cg->Mult(B, X);
708 delete M;
709 delete cg;
710 a->RecoverFEMSolution(X, *b, *x);
711 *u+=*x;
712 delete a;
713 delete x;
714}
715
716GridFunction * LinearElasticitySolver::GetFEMSolution()
717{
718 return u;
719}
720
721#ifdef MFEM_USE_MPI
722ParGridFunction * LinearElasticitySolver::GetParFEMSolution()
723{
724 if (parallel)
725 {
726 return dynamic_cast<ParGridFunction*>(u);
727 }
728 else
729 {
730 MFEM_ABORT("Wrong code path. Call GetFEMSolution");
731 return nullptr;
732 }
733}
734#endif
735
736LinearElasticitySolver::~LinearElasticitySolver()
737{
738 delete u; u = nullptr;
739 delete fes; fes = nullptr;
740#ifdef MFEM_USE_MPI
741 delete pfes; pfes=nullptr;
742#endif
743 delete fec; fec = nullptr;
744 delete b;
745}
746
747} // namespace mfem
748
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:42
Data type dense matrix using column-major storage.
Definition: densemat.hpp:24
Returns f(u(x)) - f(v(x)) where u, v are scalar GridFunctions and f:R → R.
Definition: ex37.hpp:66
GridFunctionCoefficient OtherGridF_cf
Definition: ex37.hpp:69
void SetFunction(std::function< real_t(const real_t)> fun_)
Definition: ex37.hpp:93
DiffMappedGridFunctionCoefficient(const GridFunction *gf, const GridFunction *other_gf, std::function< real_t(const real_t)> fun_, int comp=1)
Definition: ex37.hpp:77
const GridFunction * OtherGridF
Definition: ex37.hpp:68
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
Definition: ex37.hpp:86
std::function< real_t(const real_t)> fun
Definition: ex37.hpp:70
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:31
virtual real_t GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: gridfunc.cpp:449
Class for integration point with weight.
Definition: intrules.hpp:35
Returns f(u(x)) where u is a scalar GridFunction and f:R → R.
Definition: ex37.hpp:41
void SetFunction(std::function< real_t(const real_t)> fun_)
Definition: ex37.hpp:60
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
Definition: ex37.hpp:55
std::function< real_t(const real_t)> fun
Definition: ex37.hpp:43
MappedGridFunctionCoefficient(const GridFunction *gf, std::function< real_t(const real_t)> fun_, int comp=1)
Definition: ex37.hpp:48
Solid isotropic material penalization (SIMP) coefficient.
Definition: ex37.hpp:98
SIMPInterpolationCoefficient(GridFunction *rho_filter_, real_t min_val_=1e-6, real_t max_val_=1.0, real_t exponent_=3)
Definition: ex37.hpp:106
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
Definition: ex37.hpp:111
Strain energy density coefficient.
Definition: ex37.hpp:122
StrainEnergyDensityCoefficient(Coefficient *lambda_, Coefficient *mu_, GridFunction *u_, GridFunction *rho_filter_, real_t rho_min_=1e-6, real_t exponent_=3.0)
Definition: ex37.hpp:133
real_t sigmoid(real_t x)
Sigmoid function.
Definition: ex37.hpp:20
real_t inv_sigmoid(real_t x)
Inverse sigmoid function.
Definition: ex37.hpp:12
float real_t
Definition: config.hpp:43
real_t der_sigmoid(real_t x)
Derivative of sigmoid function.
Definition: ex37.hpp:33