MFEM  v4.6.0
Finite element discretization library
hooke.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 // This miniapp solves a quasistatic solid mechanics problem assuming an elastic
13 // material and no body forces.
14 //
15 // The equation
16 // ∇⋅σ(∇u) = 0
17 //
18 // with stress σ is solved for displacement u.
19 //
20 // +----------+----------+
21 // fixed --->| |<--- constant displacement
22 // | |
23 // +----------+----------+
24 //
25 // This miniapp uses an elasticity operator that allows for a custom material.
26 // By default the NeoHookeanMaterial is used. A linear elastic material is also
27 // provided. Based on these examples, other materials could be implemented.
28 //
29 // The implementation of NeoHookeanMaterial also demonstrates the use of
30 // automatic differentiation using either a native dual number forward mode
31 // implementation or leveraging the Enzyme third party library.
32 
33 #include <mfem.hpp>
34 
36 #include "materials/neohookean.hpp"
40 
41 using namespace std;
42 using namespace mfem;
43 
44 /// This example only works in 3D. Kernels for 2D are not implemented.
45 constexpr int dimension = 3;
46 
47 void display_banner(ostream& os)
48 {
49  os << R"(
50  ___ ___ ________ ________ ____ __.___________
51  / | \\_____ \ \_____ \ | |/ _|\_ _____/
52  / ~ \/ | \ / | \| < | __)_
53  \ Y / | \/ | \ | \ | \
54  \___|_ /\_______ /\_______ /____|__ \/_______ /
55  \/ \/ \/ \/ \/
56  )"
57  << endl << flush;
58 }
59 
60 int main(int argc, char *argv[])
61 {
62  Mpi::Init(argc, argv);
63  int num_procs = Mpi::WorldSize();
64  int myid = Mpi::WorldRank();
65 
66  int order = 1;
67  const char *device_config = "cpu";
68  int diagpc_type = ElasticityDiagonalPreconditioner::Type::Diagonal;
69  int serial_refinement_levels = 0;
70  bool visualization = true;
71  bool paraview = false;
72 
73  if (Mpi::Root())
74  {
76  }
77 
78  OptionsParser args(argc, argv);
79  args.AddOption(&order, "-o", "--order",
80  "Finite element order (polynomial degree).");
81  args.AddOption(&device_config, "-d", "--device",
82  "Device configuration string, see Device::Configure().");
83  args.AddOption(&diagpc_type, "-pc", "--pctype",
84  "Select diagonal preconditioner type"
85  " (0:Diagonal, 1:BlockDiagonal).");
86  args.AddOption(&serial_refinement_levels, "-rs", "--ref-serial",
87  "Number of uniform refinements on the serial mesh.");
88  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
89  "--no-visualization",
90  "Enable or disable GLVis visualization.");
91  args.AddOption(&paraview, "-pv", "--paraview", "-no-pv",
92  "--no-paraview",
93  "Enable or disable ParaView DataCollection output.");
94  args.ParseCheck();
95 
96  Device device(device_config);
97  if (Mpi::Root())
98  {
99  device.Print();
100  }
101 
102  auto mesh =
103  Mesh::MakeCartesian3D(8, 2, 2, Element::HEXAHEDRON, 8.0, 1.0, 1.0);
104  if (mesh.Dimension() != dimension)
105  {
106  MFEM_ABORT("This example only works in 3D.");
107  }
108  mesh.EnsureNodes();
109 
110  for (int l = 0; l < serial_refinement_levels; l++)
111  {
112  mesh.UniformRefinement();
113  }
114 
115  ParMesh pmesh(MPI_COMM_WORLD, mesh);
116  mesh.Clear();
117 
118  // Create the elasticity operator on the parallel mesh.
119  ElasticityOperator elasticity_op(pmesh, order);
120 
121  // Create and set the material type. We define its GradientType during
122  // instantiation.
123 
124  // As seen in materials/gradient_type.hpp there is a choice of the
125  // GradientType with either
126  // * Symbolic (Manually derived)
127  // * EnzymeFwd
128  // * EnzymeRev
129  // * FiniteDiff
130  // * InternalFwd
132  elasticity_op.SetMaterial(material);
133 
134  // Define all essential boundaries. In this specific example, this includes
135  // all fixed and statically displaced degrees of freedom on mesh entities in
136  // the defined attributes.
137  if (pmesh.bdr_attributes.Size())
138  {
139  Array<int> ess_attr(pmesh.bdr_attributes.Max());
140  ess_attr = 0;
141  ess_attr[4] = 1;
142  ess_attr[2] = 1;
143  elasticity_op.SetEssentialAttributes(ess_attr);
144  }
145 
146  // Define all statically displaced degrees of freedom on mesh entities in the
147  // defined attributes. On these degrees of freedom (determined from the mesh
148  // attributes), a fixed displacement is prescribed.
149  if (pmesh.bdr_attributes.Size())
150  {
151  Array<int> displaced_attr(pmesh.bdr_attributes.Max());
152  displaced_attr = 0;
153  displaced_attr[2] = 1;
154  elasticity_op.SetPrescribedDisplacement(displaced_attr);
155  }
156 
157  ParGridFunction U_gf(&elasticity_op.h1_fes_);
158  U_gf = 0.0;
159 
160  Vector U;
161  U_gf.GetTrueDofs(U);
162 
163  // Prescribe a fixed displacement to the displaced degrees of freedom.
164  U.SetSubVector(elasticity_op.GetPrescribedDisplacementTDofs(), 1.0e-2);
165 
166  // Define the type of preconditioner to use for the linear solver.
168  static_cast<ElasticityDiagonalPreconditioner::Type>(diagpc_type));
169 
170  CGSolver cg(MPI_COMM_WORLD);
171  cg.SetRelTol(1e-1);
172  cg.SetMaxIter(10000);
173  cg.SetPrintLevel(2);
174  cg.SetPreconditioner(diagonal_pc);
175 
176  NewtonSolver newton(MPI_COMM_WORLD);
177  newton.SetSolver(cg);
178  newton.SetOperator(elasticity_op);
179  newton.SetRelTol(1e-6);
180  newton.SetMaxIter(10);
181  newton.SetPrintLevel(1);
182 
183  Vector zero;
184  newton.Mult(zero, U);
185 
186  U_gf.Distribute(U);
187 
188  if (visualization)
189  {
190  char vishost[] = "localhost";
191  int visport = 19916;
192  socketstream sol_sock(vishost, visport);
193  sol_sock << "parallel " << num_procs << " " << myid << "\n";
194  sol_sock.precision(8);
195  sol_sock << "solution\n" << pmesh << U_gf << flush;
196  }
197 
198  if (paraview)
199  {
200  ParaViewDataCollection pd("elasticity_output", &pmesh);
201  pd.RegisterField("solution", &U_gf);
202  pd.SetLevelsOfDetail(order);
203  pd.SetDataFormat(VTKFormat::BINARY);
204  pd.SetHighOrderOutput(true);
205  pd.SetCycle(0);
206  pd.SetTime(0.0);
207  pd.Save();
208  }
209 
210  return 0;
211 }
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:605
int visport
Conjugate gradient method.
Definition: solvers.hpp:493
void SetMaterial(const material_type &material)
Set the material type.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetDataFormat(VTKFormat fmt)
Helper class for ParaView visualization data.
Neo-Hookean material.
Definition: neohookean.hpp:40
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: solvers.hpp:689
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
STL namespace.
int material(Vector &x, Vector &xmin, Vector &xmax)
Definition: shaper.cpp:53
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
char vishost[]
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:641
void SetHighOrderOutput(bool high_order_output_)
void SetTime(double t)
Set physical time (for time-dependent simulations)
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1819
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
void SetRelTol(double rtol)
Definition: solvers.hpp:199
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 Distribute(const Vector *tv)
Definition: pgridfunc.cpp:141
const Array< int > & GetPrescribedDisplacementTDofs()
Return the T-vector degrees of freedom that have been marked as displaced.
void display_banner(ostream &os)
Definition: hooke.cpp:47
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
Definition: hooke.cpp:45
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
int main(int argc, char *argv[])
Definition: hooke.cpp:60
void SetLevelsOfDetail(int levels_of_detail_)
void SetPrescribedDisplacement(const Array< int > attr)
Set the attributes which mark the degrees of freedom that have a fixed displacement.
void SetEssentialAttributes(const Array< int > attr)
Set the essential attributes which mark degrees of freedom for the solving process.
Vector data type.
Definition: vector.hpp:58
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:152
ElasticityDiagonalPreconditioner acts as a matrix-free preconditioner for ElasticityOperator.
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
virtual void Save() override
Class for parallel grid function.
Definition: pgridfunc.hpp:32
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
ParFiniteElementSpace h1_fes_
H1 finite element space.
Class for parallel meshes.
Definition: pmesh.hpp:32
void ParseCheck(std::ostream &out=mfem::out)
Definition: optparser.cpp:255
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1808