 MFEM  v4.4.0 Finite element discretization library
ex21.cpp
Go to the documentation of this file.
1 // MFEM Example 21
2 //
3 // Compile with: make ex21
4 //
5 // Sample runs: ex21
6 // ex21 -o 3
8 // ex21 -m ../data/beam-quad.mesh -o 3
9 // ex21 -m ../data/beam-quad.mesh -o 3 -f 1
10 // ex21 -m ../data/beam-tet.mesh
11 // ex21 -m ../data/beam-tet.mesh -o 2
12 // ex21 -m ../data/beam-hex.mesh
13 // ex21 -m ../data/beam-hex.mesh -o 2
14 //
15 // Description: This is a version of Example 2 with a simple adaptive mesh
16 // refinement loop. The problem being solved is again the linear
17 // elasticity describing a multi-material cantilever beam.
18 // The problem is solved on a sequence of meshes which
19 // are locally refined in a conforming (triangles, tetrahedrons)
20 // or non-conforming (quadrilaterals, hexahedra) manner according
21 // to a simple ZZ error estimator.
22 //
23 // The example demonstrates MFEM's capability to work with both
24 // conforming and nonconforming refinements, in 2D and 3D, on
25 // linear and curved meshes. Interpolation of functions from
26 // coarse to fine meshes, as well as persistent GLVis
27 // visualization are also illustrated.
28 //
29 // We recommend viewing Examples 2 and 6 before viewing this
30 // example.
31
32 #include "mfem.hpp"
33 #include <fstream>
34 #include <iostream>
35
36 using namespace std;
37 using namespace mfem;
38
39 int main(int argc, char *argv[])
40 {
41  // 1. Parse command-line options.
42  const char *mesh_file = "../data/beam-tri.mesh";
43  int order = 1;
44  bool static_cond = false;
45  int flux_averaging = 0;
46  bool visualization = 1;
47
48  OptionsParser args(argc, argv);
50  "Mesh file to use.");
52  "Finite element order (polynomial degree).");
54  "--no-static-condensation", "Enable static condensation.");
56  "Flux averaging: 0 - global, 1 - by mesh attribute.");
58  "--no-visualization",
59  "Enable or disable GLVis visualization.");
60  args.Parse();
61  if (!args.Good())
62  {
63  args.PrintUsage(cout);
64  return 1;
65  }
66  args.PrintOptions(cout);
67
68  // 2. Read the mesh from the given mesh file. We can handle triangular,
69  // quadrilateral, tetrahedral, and hexahedral meshes with the same code.
70  Mesh mesh(mesh_file, 1, 1);
71  int dim = mesh.Dimension();
72  MFEM_VERIFY(mesh.SpaceDimension() == dim, "invalid mesh");
73
74  if (mesh.attributes.Max() < 2 || mesh.bdr_attributes.Max() < 2)
75  {
76  cerr << "\nInput mesh should have at least two materials and "
77  << "two boundary attributes! (See schematic in ex2.cpp)\n"
78  << endl;
79  return 3;
80  }
81
82  // 3. Since a NURBS mesh can currently only be refined uniformly, we need to
83  // convert it to a piecewise-polynomial curved mesh. First we refine the
84  // NURBS mesh a bit more and then project the curvature to quadratic Nodes.
85  if (mesh.NURBSext)
86  {
87  for (int i = 0; i < 2; i++)
88  {
89  mesh.UniformRefinement();
90  }
91  mesh.SetCurvature(2);
92  }
93
94  // 4. Define a finite element space on the mesh. The polynomial order is
95  // one (linear) by default, but this can be changed on the command line.
96  H1_FECollection fec(order, dim);
97  FiniteElementSpace fespace(&mesh, &fec, dim);
98
99  // 5. As in Example 2, we set up the linear form b(.) which corresponds to
100  // the right-hand side of the FEM linear system. In this case, b_i equals
101  // the boundary integral of f*phi_i where f represents a "pull down"
102  // force on the Neumann part of the boundary and phi_i are the basis
103  // functions in the finite element fespace. The force is defined by the
104  // VectorArrayCoefficient object f, which is a vector of Coefficient
105  // objects. The fact that f is non-zero on boundary attribute 2 is
106  // indicated by the use of piece-wise constants coefficient for its last
107  // component. We don't assemble the discrete problem yet, this will be
108  // done in the main loop.
110  for (int i = 0; i < dim-1; i++)
111  {
112  f.Set(i, new ConstantCoefficient(0.0));
113  }
114  {
115  Vector pull_force(mesh.bdr_attributes.Max());
116  pull_force = 0.0;
117  pull_force(1) = -1.0e-2;
118  f.Set(dim-1, new PWConstCoefficient(pull_force));
119  }
120
121  LinearForm b(&fespace);
123
124  // 6. Set up the bilinear form a(.,.) on the finite element space
125  // corresponding to the linear elasticity integrator with piece-wise
126  // constants coefficient lambda and mu.
127  Vector lambda(mesh.attributes.Max());
128  lambda = 1.0;
129  lambda(0) = lambda(1)*50;
130  PWConstCoefficient lambda_func(lambda);
131  Vector mu(mesh.attributes.Max());
132  mu = 1.0;
133  mu(0) = mu(1)*50;
134  PWConstCoefficient mu_func(mu);
135
136  BilinearForm a(&fespace);
137  BilinearFormIntegrator *integ =
138  new ElasticityIntegrator(lambda_func,mu_func);
140  if (static_cond) { a.EnableStaticCondensation(); }
141
142  // 7. The solution vector x and the associated finite element grid function
143  // will be maintained over the AMR iterations. We initialize it to zero.
144  Vector zero_vec(dim);
145  zero_vec = 0.0;
146  VectorConstantCoefficient zero_vec_coeff(zero_vec);
147  GridFunction x(&fespace);
148  x = 0.0;
149
150  // 8. Determine the list of true (i.e. conforming) essential boundary dofs.
151  // In this example, the boundary conditions are defined by marking only
152  // boundary attribute 1 from the mesh as essential and converting it to a
153  // list of true dofs. The conversion to true dofs will be done in the
154  // main loop.
155  Array<int> ess_bdr(mesh.bdr_attributes.Max());
156  ess_bdr = 0;
157  ess_bdr = 1;
158
159  // 9. Connect to GLVis.
160  char vishost[] = "localhost";
161  int visport = 19916;
162  socketstream sol_sock;
163  if (visualization)
164  {
165  sol_sock.open(vishost, visport);
166  sol_sock.precision(8);
167  }
168
169  // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
170  // that uses the ComputeElementFlux method of the ElasticityIntegrator to
171  // recover a smoothed flux (stress) that is subtracted from the element
172  // flux to get an error indicator. We need to supply the space for the
173  // smoothed flux: an (H1)^tdim (i.e., vector-valued) space is used here.
174  // Here, tdim represents the number of components for a symmetric (dim x
175  // dim) tensor.
176  const int tdim = dim*(dim+1)/2;
177  FiniteElementSpace flux_fespace(&mesh, &fec, tdim);
178  ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
179  estimator.SetFluxAveraging(flux_averaging);
180
181  // 11. A refiner selects and refines elements based on a refinement strategy.
182  // The strategy here is to refine elements with errors larger than a
183  // fraction of the maximum element error. Other strategies are possible.
184  // The refiner will call the given error estimator.
185  ThresholdRefiner refiner(estimator);
186  refiner.SetTotalErrorFraction(0.7);
187
188  // 12. The main AMR loop. In each iteration we solve the problem on the
189  // current mesh, visualize the solution, and refine the mesh.
190  const int max_dofs = 50000;
191  const int max_amr_itr = 20;
192  for (int it = 0; it <= max_amr_itr; it++)
193  {
194  int cdofs = fespace.GetTrueVSize();
195  cout << "\nAMR iteration " << it << endl;
196  cout << "Number of unknowns: " << cdofs << endl;
197
198  // 13. Assemble the stiffness matrix and the right-hand side.
199  a.Assemble();
200  b.Assemble();
201
202  // 14. Set Dirichlet boundary values in the GridFunction x.
203  // Determine the list of Dirichlet true DOFs in the linear system.
204  Array<int> ess_tdof_list;
205  x.ProjectBdrCoefficient(zero_vec_coeff, ess_bdr);
206  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
207
208  // 15. Create the linear system: eliminate boundary conditions, constrain
209  // hanging nodes and possibly apply other transformations. The system
210  // will be solved for true (unconstrained) DOFs only.
211  SparseMatrix A;
212  Vector B, X;
213  const int copy_interior = 1;
214  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
215
216 #ifndef MFEM_USE_SUITESPARSE
217  // 16. Define a simple symmetric Gauss-Seidel preconditioner and use it to
218  // solve the linear system with PCG.
219  GSSmoother M(A);
220  PCG(A, M, B, X, 3, 2000, 1e-12, 0.0);
221 #else
222  // 16. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the
223  // the linear system.
224  UMFPackSolver umf_solver;
225  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
226  umf_solver.SetOperator(A);
227  umf_solver.Mult(B, X);
228 #endif
229
230  // 17. After solving the linear system, reconstruct the solution as a
231  // finite element GridFunction. Constrained nodes are interpolated
232  // from true DOFs (it may therefore happen that x.Size() >= X.Size()).
233  a.RecoverFEMSolution(X, b, x);
234
235  // 18. Send solution by socket to the GLVis server.
236  if (visualization && sol_sock.good())
237  {
238  GridFunction nodes(&fespace), *nodes_p = &nodes;
239  mesh.GetNodes(nodes);
240  nodes += x;
241  int own_nodes = 0;
242  mesh.SwapNodes(nodes_p, own_nodes);
243  x.Neg(); // visualize the backward displacement
244  sol_sock << "solution\n" << mesh << x << flush;
245  x.Neg();
246  mesh.SwapNodes(nodes_p, own_nodes);
247  if (it == 0)
248  {
249  sol_sock << "keys '" << ((dim == 2) ? "Rjl" : "") << "m'" << endl;
250  }
251  sol_sock << "window_title 'AMR iteration: " << it << "'\n"
252  << "pause" << endl;
253  cout << "Visualization paused. "
254  "Press <space> in the GLVis window to continue." << endl;
255  }
256
257  if (cdofs > max_dofs)
258  {
259  cout << "Reached the maximum number of dofs. Stop." << endl;
260  break;
261  }
262
263  // 19. Call the refiner to modify the mesh. The refiner calls the error
264  // estimator to obtain element errors, then it selects elements to be
265  // refined and finally it modifies the mesh. The Stop() method can be
266  // used to determine if a stopping criterion was met.
267  refiner.Apply(mesh);
268  if (refiner.Stop())
269  {
270  cout << "Stopping criterion satisfied. Stop." << endl;
271  break;
272  }
273
274  // 20. Update the space to reflect the new state of the mesh. Also,
275  // interpolate the solution x so that it lies in the new space but
276  // represents the same function. This saves solver iterations later
277  // since we'll have a good initial guess of x in the next step.
278  // Internally, FiniteElementSpace::Update() calculates an
279  // interpolation matrix which is then used by GridFunction::Update().
280  fespace.Update();
281  x.Update();
282
283  // 21. Inform also the bilinear and linear forms that the space has
284  // changed.
285  a.Update();
286  b.Update();
287  }
288
289  {
290  ofstream mesh_ref_out("ex21_reference.mesh");
291  mesh_ref_out.precision(16);
292  mesh.Print(mesh_ref_out);
293
294  ofstream mesh_out("ex21_deformed.mesh");
295  mesh_out.precision(16);
296  GridFunction nodes(&fespace), *nodes_p = &nodes;
297  mesh.GetNodes(nodes);
298  nodes += x;
299  int own_nodes = 0;
300  mesh.SwapNodes(nodes_p, own_nodes);
301  mesh.Print(mesh_out);
302  mesh.SwapNodes(nodes_p, own_nodes);
303
304  ofstream x_out("ex21_displacement.sol");
305  x_out.precision(16);
306  x.Save(x_out);
307  }
308
309  return 0;
310 }
Vector coefficient defined by an array of scalar coefficients. Coefficients that are not set will eva...
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3336
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
Vector coefficient that is constant in space and time.
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Definition: mesh.cpp:7884
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:102
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition: fespace.cpp:564
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:1025
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3619
bool Apply(Mesh &mesh)
Perform the mesh operation.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
Data type sparse matrix.
Definition: sparsemat.hpp:46
constexpr char vishost[]
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
constexpr int visport
Mesh refinement operator using an error threshold.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:5312
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void SetFluxAveraging(int fa)
Set the way the flux is averaged (smoothed) across elements.
Definition: estimators.hpp:173
void Set(int i, Coefficient *c, bool own=true)
Sets coefficient in the vector.
virtual void Update(FiniteElementSpace *nfes=NULL)
Update the FiniteElementSpace and delete all data associated with the old one.
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition: solvers.cpp:904
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:566
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
int SpaceDimension() const
Definition: mesh.hpp:1000
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:34
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:1036
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
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
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:164
void Update()
Update the object according to the associated FE space fes.
Definition: linearform.hpp:188
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
A piecewise constant coefficient with the constants keyed off the element attribute numbers...
Definition: coefficient.hpp:94
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1646
double mu
Definition: ex25.cpp:139
int dim
Definition: ex24.cpp:53
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
The ZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu error estimation procedure...
Definition: estimators.hpp:88
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7841
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:23
void SetTotalErrorFraction(double fraction)
Set the total fraction used in the computation of the threshold. The default value is 1/2...
int main()
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:268
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3179
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
Definition: gridfunc.hpp:447
double f(const Vector &p)
void Neg()
(*this) = -(*this)
Definition: vector.cpp:292
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:3084
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150