MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex21p.cpp
Go to the documentation of this file.
1// MFEM Example 21 - Parallel Version
2//
3// Compile with: make ex21p
4//
5// Sample runs: mpirun -np 4 ex21p
6// mpirun -np 4 ex21p -o 3
7// mpirun -np 4 ex21p -m ../data/beam-quad.mesh
8// mpirun -np 4 ex21p -m ../data/beam-quad.mesh -o 3
9// mpirun -np 4 ex21p -m ../data/beam-tet.mesh
10// mpirun -np 4 ex21p -m ../data/beam-tet.mesh -o 2
11// mpirun -np 4 ex21p -m ../data/beam-hex.mesh
12// mpirun -np 4 ex21p -m ../data/beam-hex.mesh -o 2
13//
14// Description: This is a version of Example 2p with a simple adaptive mesh
15// refinement loop. The problem being solved is again the linear
16// elasticity describing a multi-material cantilever beam.
17// The problem is solved on a sequence of meshes which
18// are locally refined in a conforming (triangles, tetrahedrons)
19// or non-conforming (quadrilaterals, hexahedra) manner according
20// to a simple ZZ error estimator.
21//
22// The example demonstrates MFEM's capability to work with both
23// conforming and nonconforming refinements, in 2D and 3D, on
24// linear and curved meshes. Interpolation of functions from
25// coarse to fine meshes, as well as persistent GLVis
26// visualization are also illustrated.
27//
28// We recommend viewing Examples 2p and 6p before viewing this
29// example.
30
31#include "mfem.hpp"
32#include <fstream>
33#include <iostream>
34
35using namespace std;
36using namespace mfem;
37
38int main(int argc, char *argv[])
39{
40 // 0. Initialize MPI and HYPRE.
41 Mpi::Init(argc, argv);
42 int num_procs = Mpi::WorldSize();
43 int myid = Mpi::WorldRank();
45
46 // 1. Parse command-line options.
47 const char *mesh_file = "../data/beam-tri.mesh";
48 int serial_ref_levels = 0;
49 int order = 1;
50 bool static_cond = false;
51 bool visualization = 1;
52
53 OptionsParser args(argc, argv);
55 "Mesh file to use.");
57 "Number of uniform serial refinements (before parallel"
58 " partitioning)");
60 "Finite element order (polynomial degree).");
62 "--no-static-condensation", "Enable static condensation.");
64 "--no-visualization",
65 "Enable or disable GLVis visualization.");
66 args.Parse();
67 if (!args.Good())
68 {
69 if (myid == 0)
70 {
71 args.PrintUsage(cout);
72 }
73 return 1;
74 }
75 if (myid == 0)
76 {
77 args.PrintOptions(cout);
78 }
79
80 // 2. Read the mesh from the given mesh file. We can handle triangular,
81 // quadrilateral, tetrahedral, and hexahedral meshes with the same code.
82 Mesh mesh(mesh_file, 1, 1);
83 int dim = mesh.Dimension();
84 MFEM_VERIFY(mesh.SpaceDimension() == dim, "invalid mesh");
85
86 if (mesh.attributes.Max() < 2 || mesh.bdr_attributes.Max() < 2)
87 {
88 cerr << "\nInput mesh should have at least two materials and "
89 << "two boundary attributes! (See schematic in ex2.cpp)\n"
90 << endl;
91 return 3;
92 }
93
94 // 3. Refine the mesh before parallel partitioning. Since a NURBS mesh can
95 // currently only be refined uniformly, we need to convert it to a
96 // piecewise-polynomial curved mesh. First we refine the NURBS mesh a bit
97 // more and then project the curvature to quadratic Nodes.
98 if (mesh.NURBSext && serial_ref_levels == 0)
99 {
100 serial_ref_levels = 2;
101 }
102 for (int i = 0; i < serial_ref_levels; i++)
103 {
104 mesh.UniformRefinement();
105 }
106 if (mesh.NURBSext)
107 {
108 mesh.SetCurvature(2);
109 }
110 mesh.EnsureNCMesh();
111
112 ParMesh pmesh(MPI_COMM_WORLD, mesh);
113 mesh.Clear();
114
115 // 4. Define a finite element space on the mesh. The polynomial order is
116 // one (linear) by default, but this can be changed on the command line.
117 H1_FECollection fec(order, dim);
118 ParFiniteElementSpace fespace(&pmesh, &fec, dim);
119
120 // 5. As in Example 2, we set up the linear form b(.) which corresponds to
121 // the right-hand side of the FEM linear system. In this case, b_i equals
122 // the boundary integral of f*phi_i where f represents a "pull down"
123 // force on the Neumann part of the boundary and phi_i are the basis
124 // functions in the finite element fespace. The force is defined by the
125 // VectorArrayCoefficient object f, which is a vector of Coefficient
126 // objects. The fact that f is non-zero on boundary attribute 2 is
127 // indicated by the use of piece-wise constants coefficient for its last
128 // component. We don't assemble the discrete problem yet, this will be
129 // done in the main loop.
131 for (int i = 0; i < dim-1; i++)
132 {
133 f.Set(i, new ConstantCoefficient(0.0));
134 }
135 {
136 Vector pull_force(pmesh.bdr_attributes.Max());
137 pull_force = 0.0;
138 pull_force(1) = -1.0e-2;
139 f.Set(dim-1, new PWConstCoefficient(pull_force));
140 }
141
142 ParLinearForm b(&fespace);
144
145 // 6. Set up the bilinear form a(.,.) on the finite element space
146 // corresponding to the linear elasticity integrator with piece-wise
147 // constants coefficient lambda and mu.
148 Vector lambda(pmesh.attributes.Max());
149 lambda = 1.0;
150 lambda(0) = lambda(1)*50;
151 PWConstCoefficient lambda_func(lambda);
152 Vector mu(pmesh.attributes.Max());
153 mu = 1.0;
154 mu(0) = mu(1)*50;
155 PWConstCoefficient mu_func(mu);
156
157 ParBilinearForm a(&fespace);
159 new ElasticityIntegrator(lambda_func,mu_func);
161 if (static_cond) { a.EnableStaticCondensation(); }
162
163 // 7. The solution vector x and the associated finite element grid function
164 // will be maintained over the AMR iterations. We initialize it to zero.
165 Vector zero_vec(dim);
166 zero_vec = 0.0;
167 VectorConstantCoefficient zero_vec_coeff(zero_vec);
168 ParGridFunction x(&fespace);
169 x = 0.0;
170
171 // 8. Determine the list of true (i.e. conforming) essential boundary dofs.
172 // In this example, the boundary conditions are defined by marking only
173 // boundary attribute 1 from the mesh as essential and converting it to a
174 // list of true dofs. The conversion to true dofs will be done in the
175 // main loop.
176 Array<int> ess_bdr(pmesh.bdr_attributes.Max());
177 ess_bdr = 0;
178 ess_bdr[0] = 1;
179
180 // 9. GLVis visualization.
181 char vishost[] = "localhost";
182 int visport = 19916;
183 socketstream sol_sock;
184
185 // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
186 // that uses the ComputeElementFlux method of the ElasticityIntegrator to
187 // recover a smoothed flux (stress) that is subtracted from the element
188 // flux to get an error indicator. We need to supply the space for the
189 // smoothed flux: an (H1)^tdim (i.e., vector-valued) space is used here.
190 // Here, tdim represents the number of components for a symmetric (dim x
191 // dim) tensor.
192 const int tdim = dim*(dim+1)/2;
193 L2_FECollection flux_fec(order, dim);
194 ParFiniteElementSpace flux_fespace(&pmesh, &flux_fec, tdim);
195 ParFiniteElementSpace smooth_flux_fespace(&pmesh, &fec, tdim);
196 L2ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace,
197 smooth_flux_fespace);
198
199 // 11. A refiner selects and refines elements based on a refinement strategy.
200 // The strategy here is to refine elements with errors larger than a
201 // fraction of the maximum element error. Other strategies are possible.
202 // The refiner will call the given error estimator.
203 ThresholdRefiner refiner(estimator);
204 refiner.SetTotalErrorFraction(0.7);
205
206 // 12. The main AMR loop. In each iteration we solve the problem on the
207 // current mesh, visualize the solution, and refine the mesh.
208 const int max_dofs = 50000;
209 const int max_amr_itr = 20;
210 for (int it = 0; it <= max_amr_itr; it++)
211 {
212 HYPRE_BigInt global_dofs = fespace.GlobalTrueVSize();
213 if (myid == 0)
214 {
215 cout << "\nAMR iteration " << it << endl;
216 cout << "Number of unknowns: " << global_dofs << endl;
217 }
218
219 // 13. Assemble the stiffness matrix and the right-hand side.
220 a.Assemble();
221 b.Assemble();
222
223 // 14. Set Dirichlet boundary values in the GridFunction x.
224 // Determine the list of Dirichlet true DOFs in the linear system.
225 Array<int> ess_tdof_list;
226 x.ProjectBdrCoefficient(zero_vec_coeff, ess_bdr);
227 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
228
229 // 15. Create the linear system: eliminate boundary conditions, constrain
230 // hanging nodes and possibly apply other transformations. The system
231 // will be solved for true (unconstrained) DOFs only.
232
234 Vector B, X;
235 const int copy_interior = 1;
236 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
237
238 // 16. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
239 // preconditioner from hypre.
240 HypreBoomerAMG amg;
241 amg.SetPrintLevel(0);
242 // amg.SetSystemsOptions(dim); // optional
243 CGSolver pcg(A.GetComm());
244 pcg.SetPreconditioner(amg);
245 pcg.SetOperator(A);
246 pcg.SetRelTol(1e-6);
247 pcg.SetMaxIter(500);
248 pcg.SetPrintLevel(3); // print the first and the last iterations only
249 pcg.Mult(B, X);
250
251 // 17. After solving the linear system, reconstruct the solution as a
252 // finite element GridFunction. Constrained nodes are interpolated
253 // from true DOFs (it may therefore happen that x.Size() >= X.Size()).
254 a.RecoverFEMSolution(X, b, x);
255
256 // 18. Send solution by socket to the GLVis server.
257 if (visualization && it == 0)
258 {
259 sol_sock.open(vishost, visport);
260 sol_sock.precision(8);
261 }
262 if (visualization && sol_sock.good())
263 {
264 GridFunction nodes(&fespace), *nodes_p = &nodes;
265 pmesh.GetNodes(nodes);
266 nodes += x;
267 int own_nodes = 0;
268 pmesh.SwapNodes(nodes_p, own_nodes);
269 x.Neg(); // visualize the backward displacement
270 sol_sock << "parallel " << num_procs << ' ' << myid << '\n';
271 sol_sock << "solution\n" << pmesh << x << flush;
272 x.Neg();
273 pmesh.SwapNodes(nodes_p, own_nodes);
274 if (it == 0)
275 {
276 sol_sock << "keys '" << ((dim == 2) ? "Rjl" : "") << "m'" << endl;
277 }
278 sol_sock << "window_title 'AMR iteration: " << it << "'\n"
279 << "pause" << endl;
280 if (myid == 0)
281 {
282 cout << "Visualization paused. "
283 "Press <space> in the GLVis window to continue." << endl;
284 }
285 }
286
287 if (global_dofs > max_dofs)
288 {
289 if (myid == 0)
290 {
291 cout << "Reached the maximum number of dofs. Stop." << endl;
292 }
293 break;
294 }
295
296 // 19. Call the refiner to modify the mesh. The refiner calls the error
297 // estimator to obtain element errors, then it selects elements to be
298 // refined and finally it modifies the mesh. The Stop() method can be
299 // used to determine if a stopping criterion was met.
300 refiner.Apply(pmesh);
301 if (refiner.Stop())
302 {
303 if (myid == 0)
304 {
305 cout << "Stopping criterion satisfied. Stop." << endl;
306 }
307 break;
308 }
309
310 // 20. Update the space to reflect the new state of the mesh. Also,
311 // interpolate the solution x so that it lies in the new space but
312 // represents the same function. This saves solver iterations later
313 // since we'll have a good initial guess of x in the next step.
314 // Internally, FiniteElementSpace::Update() calculates an
315 // interpolation matrix which is then used by GridFunction::Update().
316 fespace.Update();
317 x.Update();
318
319 // 21. Load balance the mesh, and update the space and solution. Currently
320 // available only for nonconforming meshes.
321 if (pmesh.Nonconforming())
322 {
323 pmesh.Rebalance();
324
325 // Update the space and the GridFunction. This time the update matrix
326 // redistributes the GridFunction among the processors.
327 fespace.Update();
328 x.Update();
329 }
330
331 // 21. Inform also the bilinear and linear forms that the space has
332 // changed.
333 a.Update();
334 b.Update();
335 }
336
337 {
338 ostringstream mref_name, mesh_name, sol_name;
339 mref_name << "ex21p_reference_mesh." << setfill('0') << setw(6) << myid;
340 mesh_name << "ex21p_deformed_mesh." << setfill('0') << setw(6) << myid;
341 sol_name << "ex21p_displacement." << setfill('0') << setw(6) << myid;
342
343 ofstream mesh_ref_out(mref_name.str().c_str());
344 mesh_ref_out.precision(16);
345 pmesh.Print(mesh_ref_out);
346
347 ofstream mesh_out(mesh_name.str().c_str());
348 mesh_out.precision(16);
349 GridFunction nodes(&fespace), *nodes_p = &nodes;
350 pmesh.GetNodes(nodes);
351 nodes += x;
352 int own_nodes = 0;
353 pmesh.SwapNodes(nodes_p, own_nodes);
354 pmesh.Print(mesh_out);
355 pmesh.SwapNodes(nodes_p, own_nodes);
356
357 ofstream x_out(sol_name.str().c_str());
358 x_out.precision(16);
359 x.Save(x_out);
360 }
361
362 return 0;
363}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
Abstract base class BilinearFormIntegrator.
Definition solvers.hpp:513
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:526
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:718
A coefficient that is constant across space and time.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
MPI_Comm GetComm() const
MPI communicator.
Definition hypre.hpp:578
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
The L2ZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu error estimation procedure where t...
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
bool Apply(Mesh &mesh)
Perform the mesh operation.
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
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
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:290
bool Nonconforming() const
Definition mesh.hpp:2229
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition mesh.cpp:10626
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition mesh.cpp:6211
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Swap the internal node GridFunction pointer and ownership flag members with the given ones.
Definition mesh.cpp:9022
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
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.
A piecewise constant coefficient with the constants keyed off the element attribute numbers.
Class for parallel bilinear form.
Abstract parallel finite element space.
Definition pfespace.hpp:29
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
void Update(bool want_transform=true) override
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Save(std::ostream &out) const override
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, const Array< int > &attr)
void Update() override
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition pgridfunc.cpp:90
Class for parallel linear form.
Class for parallel meshes.
Definition pmesh.hpp:34
void Rebalance()
Definition pmesh.cpp:4009
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Mesh refinement operator using an error threshold.
void SetTotalErrorFraction(real_t fraction)
Set the total fraction used in the computation of the threshold. The default value is 1/2.
Vector coefficient defined by an array of scalar coefficients. Coefficients that are not set will eva...
Vector coefficient that is constant in space and time.
Vector data type.
Definition vector.hpp:80
void Neg()
(*this) = -(*this)
Definition vector.cpp:300
Vector & Set(const real_t a, const Vector &x)
(*this) = a * x
Definition vector.cpp:262
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
int dim
Definition ex24.cpp:53
real_t mu
Definition ex25.cpp:140
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
const char vishost[]