MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex7p.cpp
Go to the documentation of this file.
1// MFEM Example 7 - Parallel Version
2//
3// Compile with: make ex7p
4//
5// Sample runs: mpirun -np 4 ex7p -e 0 -o 2 -r 4
6// mpirun -np 4 ex7p -e 1 -o 2 -r 4 -snap
7// mpirun -np 4 ex7p -e 0 -amr 1
8// mpirun -np 4 ex7p -e 1 -amr 2 -o 2
9//
10// Description: This example code demonstrates the use of MFEM to define a
11// triangulation of a unit sphere and a simple isoparametric
12// finite element discretization of the Laplace problem with mass
13// term, -Delta u + u = f.
14//
15// The example highlights mesh generation, the use of mesh
16// refinement, high-order meshes and finite elements, as well as
17// surface-based linear and bilinear forms corresponding to the
18// left-hand side and right-hand side of the discrete linear
19// system. Simple local mesh refinement is also demonstrated.
20//
21// We recommend viewing Example 1 before viewing this example.
22
23#include "mfem.hpp"
24#include <fstream>
25#include <iostream>
26
27using namespace std;
28using namespace mfem;
29
30// Exact solution and r.h.s., see below for implementation.
32real_t analytic_rhs(const Vector &x);
33void SnapNodes(Mesh &mesh);
34
35int main(int argc, char *argv[])
36{
37 // 1. Initialize MPI and HYPRE.
38 Mpi::Init(argc, argv);
39 int num_procs = Mpi::WorldSize();
40 int myid = Mpi::WorldRank();
42
43 // 2. Parse command-line options.
44 int elem_type = 1;
45 int ref_levels = 2;
46 int amr = 0;
47 int order = 2;
48 bool always_snap = false;
49 bool visualization = 1;
50 const char *device_config = "cpu";
51
52 OptionsParser args(argc, argv);
54 "Type of elements to use: 0 - triangles, 1 - quads.");
56 "Finite element order (polynomial degree).");
58 "Number of times to refine the mesh uniformly.");
61 " 1 = refine around north pole, 2 = refine randomly.");
63 "--no-visualization",
64 "Enable or disable GLVis visualization.");
66 "--snap-at-the-end",
67 "If true, snap nodes to the sphere initially and after each refinement "
68 "otherwise, snap only after the last refinement");
70 "Device configuration string, see Device::Configure().");
71 args.Parse();
72 if (!args.Good())
73 {
74 if (myid == 0)
75 {
76 args.PrintUsage(cout);
77 }
78 return 1;
79 }
80 if (myid == 0)
81 {
82 args.PrintOptions(cout);
83 }
84
85 // 3. Enable hardware devices such as GPUs, and programming models such as
86 // CUDA, OCCA, RAJA and OpenMP based on command line options.
87 Device device(device_config);
88 if (myid == 0) { device.Print(); }
89
90 // 4. Generate an initial high-order (surface) mesh on the unit sphere. The
91 // Mesh object represents a 2D mesh in 3 spatial dimensions. We first add
92 // the elements and the vertices of the mesh, and then make it high-order
93 // by specifying a finite element space for its nodes.
94 int Nvert = 8, Nelem = 6;
95 if (elem_type == 0)
96 {
97 Nvert = 6;
98 Nelem = 8;
99 }
100 Mesh *mesh = new Mesh(2, Nvert, Nelem, 0, 3);
101
102 if (elem_type == 0) // inscribed octahedron
103 {
104 const real_t tri_v[6][3] =
105 {
106 { 1, 0, 0}, { 0, 1, 0}, {-1, 0, 0},
107 { 0, -1, 0}, { 0, 0, 1}, { 0, 0, -1}
108 };
109 const int tri_e[8][3] =
110 {
111 {0, 1, 4}, {1, 2, 4}, {2, 3, 4}, {3, 0, 4},
112 {1, 0, 5}, {2, 1, 5}, {3, 2, 5}, {0, 3, 5}
113 };
114
115 for (int j = 0; j < Nvert; j++)
116 {
118 }
119 for (int j = 0; j < Nelem; j++)
120 {
121 int attribute = j + 1;
123 }
124 mesh->FinalizeTriMesh(1, 1, true);
125 }
126 else // inscribed cube
127 {
129 {
130 {-1, -1, -1}, {+1, -1, -1}, {+1, +1, -1}, {-1, +1, -1},
131 {-1, -1, +1}, {+1, -1, +1}, {+1, +1, +1}, {-1, +1, +1}
132 };
134 {
135 {3, 2, 1, 0}, {0, 1, 5, 4}, {1, 2, 6, 5},
136 {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}
137 };
138
139 for (int j = 0; j < Nvert; j++)
140 {
142 }
143 for (int j = 0; j < Nelem; j++)
144 {
145 int attribute = j + 1;
147 }
149 }
150
151 // Set the space for the high-order mesh nodes.
152 H1_FECollection fec(order, mesh->Dimension());
153 FiniteElementSpace nodal_fes(mesh, &fec, mesh->SpaceDimension());
154 mesh->SetNodalFESpace(&nodal_fes);
155
156 // 5. Refine the mesh while snapping nodes to the sphere. Number of parallel
157 // refinements is fixed to 2.
158 for (int l = 0; l <= ref_levels; l++)
159 {
160 if (l > 0) // for l == 0 just perform snapping
161 {
162 mesh->UniformRefinement();
163 }
164
165 // Snap the nodes of the refined mesh back to sphere surface.
166 if (always_snap)
167 {
168 SnapNodes(*mesh);
169 }
170 }
171
172 if (amr == 1)
173 {
174 Vertex target(0.0, 0.0, 1.0);
175 for (int l = 0; l < 3; l++)
176 {
177 mesh->RefineAtVertex(target);
178 }
179 SnapNodes(*mesh);
180 }
181 else if (amr == 2)
182 {
183 for (int l = 0; l < 2; l++)
184 {
185 mesh->RandomRefinement(0.5); // 50% probability
186 }
187 SnapNodes(*mesh);
188 }
189
190 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
191 delete mesh;
192 {
193 int par_ref_levels = 2;
194 for (int l = 0; l < par_ref_levels; l++)
195 {
196 pmesh->UniformRefinement();
197
198 // Snap the nodes of the refined mesh back to sphere surface.
199 if (always_snap)
200 {
201 SnapNodes(*pmesh);
202 }
203 }
204 if (!always_snap || par_ref_levels < 1)
205 {
206 SnapNodes(*pmesh);
207 }
208 }
209
210 if (amr == 1)
211 {
212 Vertex target(0.0, 0.0, 1.0);
213 for (int l = 0; l < 2; l++)
214 {
215 pmesh->RefineAtVertex(target);
216 }
217 SnapNodes(*pmesh);
218 }
219 else if (amr == 2)
220 {
221 for (int l = 0; l < 2; l++)
222 {
223 pmesh->RandomRefinement(0.5); // 50% probability
224 }
225 SnapNodes(*pmesh);
226 }
227
228 // 6. Define a finite element space on the mesh. Here we use isoparametric
229 // finite elements -- the same as the mesh nodes.
230 ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, &fec);
231 HYPRE_BigInt size = fespace->GlobalTrueVSize();
232 if (myid == 0)
233 {
234 cout << "Number of unknowns: " << size << endl;
235 }
236
237 // 7. Set up the linear form b(.) which corresponds to the right-hand side of
238 // the FEM linear system, which in this case is (1,phi_i) where phi_i are
239 // the basis functions in the finite element fespace.
240 ParLinearForm *b = new ParLinearForm(fespace);
241 ConstantCoefficient one(1.0);
245 b->Assemble();
246
247 // 8. Define the solution vector x as a finite element grid function
248 // corresponding to fespace. Initialize x with initial guess of zero.
249 ParGridFunction x(fespace);
250 x = 0.0;
251
252 // 9. Set up the bilinear form a(.,.) on the finite element space
253 // corresponding to the Laplacian operator -Delta, by adding the Diffusion
254 // and Mass domain integrators.
255 ParBilinearForm *a = new ParBilinearForm(fespace);
258
259 // 10. Assemble the parallel linear system, applying any transformations
260 // such as: parallel assembly, applying conforming constraints, etc.
261 a->Assemble();
263 Vector B, X;
264 Array<int> empty_tdof_list;
265 a->FormLinearSystem(empty_tdof_list, x, *b, A, X, B);
266
267 // 11. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
268 // preconditioner from hypre. Extract the parallel grid function x
269 // corresponding to the finite element approximation X. This is the local
270 // solution on each processor.
271 HypreSolver *amg = new HypreBoomerAMG(A);
272 HyprePCG *pcg = new HyprePCG(A);
273 pcg->SetTol(1e-12);
274 pcg->SetMaxIter(200);
275 pcg->SetPrintLevel(2);
276 pcg->SetPreconditioner(*amg);
277 pcg->Mult(B, X);
278 a->RecoverFEMSolution(X, *b, x);
279
280 delete a;
281 delete b;
282
283 // 12. Compute and print the L^2 norm of the error.
284 real_t error = x.ComputeL2Error(sol_coef);
285 if (myid == 0)
286 {
287 cout << "\nL2 norm of error: " << error << endl;
288 }
289
290 // 13. Save the refined mesh and the solution. This output can be viewed
291 // later using GLVis: "glvis -np <np> -m sphere_refined -g sol".
292 {
293 ostringstream mesh_name, sol_name;
294 mesh_name << "sphere_refined." << setfill('0') << setw(6) << myid;
295 sol_name << "sol." << setfill('0') << setw(6) << myid;
296
297 ofstream mesh_ofs(mesh_name.str().c_str());
298 mesh_ofs.precision(8);
299 pmesh->Print(mesh_ofs);
300
301 ofstream sol_ofs(sol_name.str().c_str());
302 sol_ofs.precision(8);
303 x.Save(sol_ofs);
304 }
305
306 // 14. Send the solution by socket to a GLVis server.
307 if (visualization)
308 {
309 char vishost[] = "localhost";
310 int visport = 19916;
311 socketstream sol_sock(vishost, visport);
312 sol_sock << "parallel " << num_procs << " " << myid << "\n";
313 sol_sock.precision(8);
314 sol_sock << "solution\n" << *pmesh << x << flush;
315 }
316
317 // 15. Free the used memory.
318 delete pcg;
319 delete amg;
320 delete fespace;
321 delete pmesh;
322
323 return 0;
324}
325
327{
328 real_t l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
329 return x(0)*x(1)/l2;
330}
331
333{
334 real_t l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
335 return 7*x(0)*x(1)/l2;
336}
337
338void SnapNodes(Mesh &mesh)
339{
340 GridFunction &nodes = *mesh.GetNodes();
341 Vector node(mesh.SpaceDimension());
342 for (int i = 0; i < nodes.FESpace()->GetNDofs(); i++)
343 {
344 for (int d = 0; d < mesh.SpaceDimension(); d++)
345 {
346 node(d) = nodes(nodes.FESpace()->DofToVDof(i, d));
347 }
348
349 node /= node.Norml2();
350
351 for (int d = 0; d < mesh.SpaceDimension(); d++)
352 {
353 nodes(nodes.FESpace()->DofToVDof(i, d)) = node(d);
354 }
355 }
356 if (mesh.Nonconforming())
357 {
358 // Snap hanging nodes to the master side.
359 Vector tnodes;
360 nodes.GetTrueDofs(tnodes);
361 nodes.SetFromTrueDofs(tnodes);
362 }
363}
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
Class for domain integration .
Definition lininteg.hpp:109
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Definition fespace.cpp:249
A general function coefficient.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition gridfunc.cpp:381
FiniteElementSpace * FESpace()
Definition gridfunc.hpp:696
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction.
Definition gridfunc.cpp:366
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
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
Abstract class for hypre's solvers and preconditioners.
Definition hypre.hpp:1162
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
Adds a quadrilateral to the mesh given by 4 vertices v1 through v4.
Definition mesh.cpp:1743
bool Nonconforming() const
Definition mesh.hpp:2229
int AddTriangle(int v1, int v2, int v3, int attr=1)
Adds a triangle to the mesh given by 3 vertices v1 through v3.
Definition mesh.cpp:1729
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
Definition mesh.cpp:1658
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void RandomRefinement(real_t prob, bool aniso=false, int nonconforming=-1, int nc_limit=0)
Refine each element with given probability. Uses GeneralRefinement.
Definition mesh.cpp:10650
void RefineAtVertex(const Vertex &vert, real_t eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex. Uses GeneralRefinement.
Definition mesh.cpp:10669
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition mesh.cpp:2177
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 FinalizeTriMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a triangular Mesh.
Definition mesh.cpp:2148
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition mesh.cpp:6153
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
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.
Class for parallel bilinear form.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Save(std::ostream &out) const override
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
Class for parallel linear form.
Class for parallel meshes.
Definition pmesh.hpp:34
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Vector data type.
Definition vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
Data type for vertex.
Definition vertex.hpp:23
void SnapNodes(Mesh &mesh)
Definition ex7p.cpp:338
real_t analytic_rhs(const Vector &x)
Definition ex7p.cpp:332
real_t analytic_solution(const Vector &x)
Definition ex7p.cpp:326
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
float real_t
Definition config.hpp:43
const char vishost[]