MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex2p.cpp
Go to the documentation of this file.
1// MFEM Example 2 - Parallel Version
2//
3// Compile with: make ex2p
4//
5// Sample runs: mpirun -np 4 ex2p -m ../data/beam-tri.mesh
6// mpirun -np 4 ex2p -m ../data/beam-quad.mesh
7// mpirun -np 4 ex2p -m ../data/beam-tet.mesh
8// mpirun -np 4 ex2p -m ../data/beam-hex.mesh
9// mpirun -np 4 ex2p -m ../data/beam-wedge.mesh
10// mpirun -np 4 ex2p -m ../data/beam-tri.mesh -o 2 -sys
11// mpirun -np 4 ex2p -m ../data/beam-quad.mesh -o 3 -elast
12// mpirun -np 4 ex2p -m ../data/beam-quad.mesh -o 3 -sc
13// mpirun -np 4 ex2p -m ../data/beam-quad-nurbs.mesh
14// mpirun -np 4 ex2p -m ../data/beam-hex-nurbs.mesh
15//
16// Description: This example code solves a simple linear elasticity problem
17// describing a multi-material cantilever beam.
18//
19// Specifically, we approximate the weak form of -div(sigma(u))=0
21// tensor corresponding to displacement field u, and lambda and mu
22// are the material Lame constants. The boundary conditions are
23// u=0 on the fixed part of the boundary with attribute 1, and
24// sigma(u).n=f on the remainder with f being a constant pull down
25// vector on boundary elements with attribute 2, and zero
26// otherwise. The geometry of the domain is assumed to be as
27// follows:
28//
29// +----------+----------+
30// boundary --->| material | material |<--- boundary
31// attribute 1 | 1 | 2 | attribute 2
32// (fixed) +----------+----------+ (pull down)
33//
34// The example demonstrates the use of high-order and NURBS vector
35// finite element spaces with the linear elasticity bilinear form,
36// meshes with curved elements, and the definition of piece-wise
37// constant and vector coefficient objects. Static condensation is
38// also illustrated.
39//
40// We recommend viewing Example 1 before viewing this example.
41
42#include "mfem.hpp"
43#include <fstream>
44#include <iostream>
45
46using namespace std;
47using namespace mfem;
48
49int main(int argc, char *argv[])
50{
51 // 1. Initialize MPI and HYPRE.
52 Mpi::Init(argc, argv);
53 int num_procs = Mpi::WorldSize();
54 int myid = Mpi::WorldRank();
56
57 // 2. Parse command-line options.
58 const char *mesh_file = "../data/beam-tri.mesh";
59 int order = 1;
60 bool static_cond = false;
61 bool visualization = 1;
62 bool amg_elast = 0;
63 bool reorder_space = false;
64 const char *device_config = "cpu";
65
66 OptionsParser args(argc, argv);
68 "Mesh file to use.");
70 "Finite element order (polynomial degree).");
72 "--amg-for-systems",
73 "Use the special AMG elasticity solver (GM/LN approaches), "
74 "or standard AMG for systems (unknown approach).");
76 "--no-static-condensation", "Enable static condensation.");
78 "--no-visualization",
79 "Enable or disable GLVis visualization.");
80 args.AddOption(&reorder_space, "-nodes", "--by-nodes", "-vdim", "--by-vdim",
81 "Use byNODES ordering of vector space instead of byVDIM");
83 "Device configuration string, see Device::Configure().");
84 args.Parse();
85 if (!args.Good())
86 {
87 if (myid == 0)
88 {
89 args.PrintUsage(cout);
90 }
91 return 1;
92 }
93 if (myid == 0)
94 {
95 args.PrintOptions(cout);
96 }
97
98 // 3. Enable hardware devices such as GPUs, and programming models such as
99 // CUDA, OCCA, RAJA and OpenMP based on command line options.
100 Device device(device_config);
101 if (myid == 0) { device.Print(); }
102
103 // 4. Read the (serial) mesh from the given mesh file on all processors. We
104 // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
105 // and volume meshes with the same code.
106 Mesh *mesh = new Mesh(mesh_file, 1, 1);
107 int dim = mesh->Dimension();
108
109 if (mesh->attributes.Max() < 2 || mesh->bdr_attributes.Max() < 2)
110 {
111 if (myid == 0)
112 cerr << "\nInput mesh should have at least two materials and "
113 << "two boundary attributes! (See schematic in ex2.cpp)\n"
114 << endl;
115 return 3;
116 }
117
118 // 5. Select the order of the finite element discretization space. For NURBS
119 // meshes, we increase the order by degree elevation.
120 if (mesh->NURBSext)
121 {
122 mesh->DegreeElevate(order, order);
123 }
124
125 // 6. Refine the serial mesh on all processors to increase the resolution. In
126 // this example we do 'ref_levels' of uniform refinement. We choose
127 // 'ref_levels' to be the largest number that gives a final mesh with no
128 // more than 1,000 elements.
129 {
130 int ref_levels =
131 (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
132 for (int l = 0; l < ref_levels; l++)
133 {
134 mesh->UniformRefinement();
135 }
136 }
137
138 // 7. Define a parallel mesh by a partitioning of the serial mesh. Refine
139 // this mesh further in parallel to increase the resolution. Once the
140 // parallel mesh is defined, the serial mesh can be deleted.
141 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
142 delete mesh;
143 {
144 int par_ref_levels = 1;
145 for (int l = 0; l < par_ref_levels; l++)
146 {
147 pmesh->UniformRefinement();
148 }
149 }
150
151 // 8. Define a parallel finite element space on the parallel mesh. Here we
152 // use vector finite elements, i.e. dim copies of a scalar finite element
153 // space. We use the ordering by vector dimension (the last argument of
154 // the FiniteElementSpace constructor) which is expected in the systems
155 // version of BoomerAMG preconditioner. For NURBS meshes, we use the
156 // (degree elevated) NURBS space associated with the mesh nodes.
158 ParFiniteElementSpace *fespace;
159 const bool use_nodal_fespace = pmesh->NURBSext && !amg_elast;
160 if (use_nodal_fespace)
161 {
162 fec = NULL;
163 fespace = (ParFiniteElementSpace *)pmesh->GetNodes()->FESpace();
164 }
165 else
166 {
167 fec = new H1_FECollection(order, dim);
168 if (reorder_space)
169 {
170 fespace = new ParFiniteElementSpace(pmesh, fec, dim, Ordering::byNODES);
171 }
172 else
173 {
174 fespace = new ParFiniteElementSpace(pmesh, fec, dim, Ordering::byVDIM);
175 }
176 }
177 HYPRE_BigInt size = fespace->GlobalTrueVSize();
178 if (myid == 0)
179 {
180 cout << "Number of finite element unknowns: " << size << endl
181 << "Assembling: " << flush;
182 }
183
184 // 9. Determine the list of true (i.e. parallel conforming) essential
185 // boundary dofs. In this example, the boundary conditions are defined by
186 // marking only boundary attribute 1 from the mesh as essential and
187 // converting it to a list of true dofs.
188 Array<int> ess_tdof_list, ess_bdr(pmesh->bdr_attributes.Max());
189 ess_bdr = 0;
190 ess_bdr[0] = 1;
191 fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
192
193 // 10. Set up the parallel linear form b(.) which corresponds to the
194 // right-hand side of the FEM linear system. In this case, b_i equals the
195 // boundary integral of f*phi_i where f represents a "pull down" force on
196 // the Neumann part of the boundary and phi_i are the basis functions in
197 // the finite element fespace. The force is defined by the object f, which
198 // is a vector of Coefficient objects. The fact that f is non-zero on
199 // boundary attribute 2 is indicated by the use of piece-wise constants
200 // coefficient for its last component.
202 for (int i = 0; i < dim-1; i++)
203 {
204 f.Set(i, new ConstantCoefficient(0.0));
205 }
206 {
207 Vector pull_force(pmesh->bdr_attributes.Max());
208 pull_force = 0.0;
209 pull_force(1) = -1.0e-2;
210 f.Set(dim-1, new PWConstCoefficient(pull_force));
211 }
212
213 ParLinearForm *b = new ParLinearForm(fespace);
215 if (myid == 0)
216 {
217 cout << "r.h.s. ... " << flush;
218 }
219 b->Assemble();
220
221 // 11. Define the solution vector x as a parallel finite element grid
222 // function corresponding to fespace. Initialize x with initial guess of
223 // zero, which satisfies the boundary conditions.
224 ParGridFunction x(fespace);
225 x = 0.0;
226
227 // 12. Set up the parallel bilinear form a(.,.) on the finite element space
228 // corresponding to the linear elasticity integrator with piece-wise
229 // constants coefficient lambda and mu.
230 Vector lambda(pmesh->attributes.Max());
231 lambda = 1.0;
232 lambda(0) = lambda(1)*50;
233 PWConstCoefficient lambda_func(lambda);
234 Vector mu(pmesh->attributes.Max());
235 mu = 1.0;
236 mu(0) = mu(1)*50;
237 PWConstCoefficient mu_func(mu);
238
239 ParBilinearForm *a = new ParBilinearForm(fespace);
241
242 // 13. Assemble the parallel bilinear form and the corresponding linear
243 // system, applying any necessary transformations such as: parallel
244 // assembly, eliminating boundary conditions, applying conforming
245 // constraints for non-conforming AMR, static condensation, etc.
246 if (myid == 0) { cout << "matrix ... " << flush; }
247 if (static_cond) { a->EnableStaticCondensation(); }
248 a->Assemble();
249
251 Vector B, X;
252 a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
253 if (myid == 0)
254 {
255 cout << "done." << endl;
256 cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
257 }
258
259 // 14. Define and apply a parallel PCG solver for A X = B with the BoomerAMG
260 // preconditioner from hypre.
261 HypreBoomerAMG *amg = new HypreBoomerAMG(A);
262 if (amg_elast && !a->StaticCondensationIsEnabled())
263 {
264 amg->SetElasticityOptions(fespace);
265 }
266 else
267 {
268 amg->SetSystemsOptions(dim, reorder_space);
269 }
270 HyprePCG *pcg = new HyprePCG(A);
271 pcg->SetTol(1e-8);
272 pcg->SetMaxIter(500);
273 pcg->SetPrintLevel(2);
274 pcg->SetPreconditioner(*amg);
275 pcg->Mult(B, X);
276
277 // 15. Recover the parallel grid function corresponding to X. This is the
278 // local finite element solution on each processor.
279 a->RecoverFEMSolution(X, *b, x);
280
281 // 16. For non-NURBS meshes, make the mesh curved based on the finite element
282 // space. This means that we define the mesh elements through a fespace
283 // based transformation of the reference element. This allows us to save
284 // the displaced mesh as a curved mesh when using high-order finite
285 // element displacement field. We assume that the initial mesh (read from
286 // the file) is not higher order curved mesh compared to the chosen FE
287 // space.
288 if (!use_nodal_fespace)
289 {
290 pmesh->SetNodalFESpace(fespace);
291 }
292
293 // 17. Save in parallel the displaced mesh and the inverted solution (which
294 // gives the backward displacements to the original grid). This output
295 // can be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
296 {
297 GridFunction *nodes = pmesh->GetNodes();
298 *nodes += x;
299 x *= -1;
300
301 ostringstream mesh_name, sol_name;
302 mesh_name << "mesh." << setfill('0') << setw(6) << myid;
303 sol_name << "sol." << setfill('0') << setw(6) << myid;
304
305 ofstream mesh_ofs(mesh_name.str().c_str());
306 mesh_ofs.precision(8);
307 pmesh->Print(mesh_ofs);
308
309 ofstream sol_ofs(sol_name.str().c_str());
310 sol_ofs.precision(8);
311 x.Save(sol_ofs);
312 }
313
314 // 18. Send the above data by socket to a GLVis server. Use the "n" and "b"
315 // keys in GLVis to visualize the displacements.
316 if (visualization)
317 {
318 char vishost[] = "localhost";
319 int visport = 19916;
320 socketstream sol_sock(vishost, visport);
321 sol_sock << "parallel " << num_procs << " " << myid << "\n";
322 sol_sock.precision(8);
323 sol_sock << "solution\n" << *pmesh << x << flush;
324 }
325
326 // 19. Free the used memory.
327 delete pcg;
328 delete amg;
329 delete a;
330 delete b;
331 if (fec)
332 {
333 delete fespace;
334 delete fec;
335 }
336 delete pmesh;
337
338 return 0;
339}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
A coefficient that is constant across space and time.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition device.hpp:123
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition device.cpp:286
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
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 SetSystemsOptions(int dim, bool order_bynodes=false)
Definition hypre.cpp:5111
void SetElasticityOptions(ParFiniteElementSpace *fespace, bool interp_refine=true)
Definition hypre.cpp:5238
PCG solver in hypre.
Definition hypre.hpp:1275
void SetPrintLevel(int print_lvl)
Definition hypre.cpp:4156
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition hypre.cpp:4161
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre's PCG.
Definition hypre.cpp:4184
void SetMaxIter(int max_iter)
Definition hypre.cpp:4146
void SetTol(real_t tol)
Definition hypre.cpp:4136
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
HYPRE_BigInt GetGlobalNumRows() const
Return the global number of rows.
Definition hypre.hpp:679
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
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
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
void DegreeElevate(int rel_degree, int degree=16)
Definition mesh.cpp:5779
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
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
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Save(std::ostream &out) const override
Class for parallel linear form.
Class for parallel meshes.
Definition pmesh.hpp:34
void SetNodalFESpace(FiniteElementSpace *nfes) override
Definition pmesh.cpp:2028
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Vector coefficient defined by an array of scalar coefficients. Coefficients that are not set will eva...
Vector data type.
Definition vector.hpp:80
Vector & Set(const real_t a, const Vector &x)
(*this) = a * x
Definition vector.cpp:262
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[]