MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex29p.cpp
Go to the documentation of this file.
1// MFEM Example 29 - Parallel Version
2//
3// Compile with: make ex29p
4//
5// Sample runs: mpirun -np 4 ex29p
6// mpirun -np 4 ex29p -sc
7// mpirun -np 4 ex29p -mt 3 -o 3 -sc
8// mpirun -np 4 ex29p -mt 3 -rs 1 -o 4 -sc
9//
10// Description: This example code demonstrates the use of MFEM to define a
11// finite element discretization of a PDE on a 2 dimensional
12// surface embedded in a 3 dimensional domain. In this case we
13// solve the Laplace problem -Div(sigma Grad u) = 1, with
14// homogeneous Dirichlet boundary conditions, where sigma is an
15// anisotropic diffusion constant defined as a 3x3 matrix
16// coefficient.
17//
18// This example demonstrates the use of finite element integrators
19// on 2D domains with 3D coefficients.
20//
21// We recommend viewing examples 1 and 7 before viewing this
22// example.
23
24#include "mfem.hpp"
25#include <fstream>
26#include <iostream>
27
28using namespace std;
29using namespace mfem;
30
31Mesh * GetMesh(int type);
32
33void trans(const Vector &x, Vector &r);
34
35void sigmaFunc(const Vector &x, DenseMatrix &s);
36
38{
39 return (0.25 * (2.0 + x[0]) - x[2]) * (x[2] + 0.25 * (2.0 + x[0]));
40}
41
42void duExact(const Vector &x, Vector &du)
43{
44 du.SetSize(3);
45 du[0] = 0.125 * (2.0 + x[0]) * x[1] * x[1];
46 du[1] = -0.125 * (2.0 + x[0]) * x[0] * x[1];
47 du[2] = -2.0 * x[2];
48}
49
50void fluxExact(const Vector &x, Vector &f)
51{
52 f.SetSize(3);
53
54 DenseMatrix s(3);
55 sigmaFunc(x, s);
56
57 Vector du(3);
58 duExact(x, du);
59
60 s.Mult(du, f);
61 f *= -1.0;
62}
63
64int main(int argc, char *argv[])
65{
66 // 1. Initialize MPI and HYPRE.
67 Mpi::Init(argc, argv);
68 int num_procs = Mpi::WorldSize();
69 int myid = Mpi::WorldRank();
71
72 // 2. Parse command-line options.
73 int order = 3;
74 int mesh_type = 4; // Default to Quadrilateral mesh
75 int mesh_order = 3;
76 int ser_ref_levels = 2;
77 int par_ref_levels = 1;
78 bool static_cond = false;
79 bool visualization = true;
80
81 OptionsParser args(argc, argv);
83 "Mesh type: 3 - Triangular, 4 - Quadrilateral.");
85 "Geometric order of the curved mesh.");
87 "Number of times to refine the mesh uniformly in serial.");
89 "Number of times to refine the mesh uniformly in parallel.");
91 "Finite element order (polynomial degree) or -1 for"
92 " isoparametric space.");
94 "--no-static-condensation", "Enable static condensation.");
96 "--no-visualization",
97 "Enable or disable GLVis visualization.");
98 args.ParseCheck();
99
100 // 3. Construct a quadrilateral or triangular mesh with the topology of a
101 // cylindrical surface.
102 Mesh *mesh = GetMesh(mesh_type);
103 int dim = mesh->Dimension();
104
105 // 4. Refine the mesh to increase the resolution. In this example we do
106 // 'ser_ref_levels' of uniform refinement.
107 for (int l = 0; l < ser_ref_levels; l++)
108 {
109 mesh->UniformRefinement();
110 }
111
112 // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
113 // this mesh further in parallel to increase the resolution. Once the
114 // parallel mesh is defined, the serial mesh can be deleted.
115 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
116 delete mesh;
117 for (int l = 0; l < par_ref_levels; l++)
118 {
119 pmesh.UniformRefinement();
120 }
121
122 // 6. Transform the mesh so that it has a more interesting geometry.
123 pmesh.SetCurvature(mesh_order);
124 pmesh.Transform(trans);
125
126 // 7. Define a finite element space on the mesh. Here we use continuous
127 // Lagrange finite elements of the specified order.
128 H1_FECollection fec(order, dim);
129 ParFiniteElementSpace fespace(&pmesh, &fec);
130 HYPRE_Int total_num_dofs = fespace.GlobalTrueVSize();
131 if (Mpi::Root())
132 {
133 cout << "Number of unknowns: " << total_num_dofs << endl;
134 }
135
136 // 8. Determine the list of true (i.e. conforming) essential boundary dofs.
137 // In this example, the boundary conditions are defined by marking all
138 // the boundary attributes from the mesh as essential (Dirichlet) and
139 // converting them to a list of true dofs.
140 Array<int> ess_tdof_list;
141 if (pmesh.bdr_attributes.Size())
142 {
143 Array<int> ess_bdr(pmesh.bdr_attributes.Max());
144 ess_bdr = 1;
145 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
146 }
147
148 // 9. Set up the linear form b(.) which corresponds to the right-hand side of
149 // the FEM linear system, which in this case is (1,phi_i) where phi_i are
150 // the basis functions in the finite element fespace.
151 ParLinearForm b(&fespace);
152 ConstantCoefficient one(1.0);
154 b.Assemble();
155
156 // 10. Define the solution vector x as a finite element grid function
157 // corresponding to fespace. Initialize x with initial guess of zero,
158 // which satisfies the boundary conditions.
159 ParGridFunction x(&fespace);
160 x = 0.0;
161
162 // 11. Set up the bilinear form a(.,.) on the finite element space
163 // corresponding to the Laplacian operator -Delta, by adding the
164 // Diffusion domain integrator.
165 ParBilinearForm a(&fespace);
169
170 // 12. Assemble the bilinear form and the corresponding linear system,
171 // applying any necessary transformations such as: eliminating boundary
172 // conditions, applying conforming constraints for non-conforming AMR,
173 // static condensation, etc.
174 if (static_cond) { a.EnableStaticCondensation(); }
175 a.Assemble();
176
177 OperatorPtr A;
178 Vector B, X;
179 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
180
181 if (myid == 0)
182 {
183 cout << "Size of linear system: "
184 << A.As<HypreParMatrix>()->GetGlobalNumRows() << endl;
185 }
186
187 // 13. Define and apply a parallel PCG solver for A X = B with the BoomerAMG
188 // preconditioner from hypre.
190 CGSolver cg(MPI_COMM_WORLD);
191 cg.SetRelTol(1e-12);
192 cg.SetMaxIter(2000);
193 cg.SetPrintLevel(1);
194 cg.SetPreconditioner(*amg);
195 cg.SetOperator(*A);
196 cg.Mult(B, X);
197 delete amg;
198
199 // 14. Recover the solution as a finite element grid function.
200 a.RecoverFEMSolution(X, b, x);
201
202 // 15. Compute error in the solution and its flux
204 real_t error = x.ComputeL2Error(uCoef);
205
206 if (myid == 0) { cout << "|u - u_h|_2 = " << error << endl; }
207
208 ParFiniteElementSpace flux_fespace(&pmesh, &fec, 3);
209 ParGridFunction flux(&flux_fespace);
210 x.ComputeFlux(*integ, flux); flux *= -1.0;
211
213 real_t flux_err = flux.ComputeL2Error(fluxCoef);
214
215 if (myid == 0) { cout << "|f - f_h|_2 = " << flux_err << endl; }
216
217 // 16. Save the refined mesh and the solution. This output can be viewed
218 // later using GLVis: "glvis -np <np> -m mesh -g sol".
219 {
220 ostringstream mesh_name, sol_name, flux_name;
221 mesh_name << "mesh." << setfill('0') << setw(6) << myid;
222 sol_name << "sol." << setfill('0') << setw(6) << myid;
223 flux_name << "flux." << setfill('0') << setw(6) << myid;
224
225 ofstream mesh_ofs(mesh_name.str().c_str());
226 mesh_ofs.precision(8);
227 pmesh.Print(mesh_ofs);
228
229 ofstream sol_ofs(sol_name.str().c_str());
230 sol_ofs.precision(8);
231 x.Save(sol_ofs);
232
233 ofstream flux_ofs(flux_name.str().c_str());
234 flux_ofs.precision(8);
235 flux.Save(flux_ofs);
236 }
237
238 // 17. Send the solution by socket to a GLVis server.
239 if (visualization)
240 {
241 char vishost[] = "localhost";
242 int visport = 19916;
243 socketstream sol_sock(vishost, visport);
244 sol_sock << "parallel " << num_procs << " " << myid << "\n";
245 sol_sock.precision(8);
246 sol_sock << "solution\n" << pmesh << x
247 << "window_title 'Solution'\n" << flush;
248
249 socketstream flux_sock(vishost, visport);
250 flux_sock << "parallel " << num_procs << " " << myid << "\n";
251 flux_sock.precision(8);
252 flux_sock << "solution\n" << pmesh << flux
253 << "keys vvv\n"
254 << "window_geometry 402 0 400 350\n"
255 << "window_title 'Flux'\n" << flush;
256 }
257
258 return 0;
259}
260
261// Defines a mesh consisting of four flat rectangular surfaces connected to form
262// a loop.
263Mesh * GetMesh(int type)
264{
265 Mesh * mesh = NULL;
266
267 if (type == 3)
268 {
269 mesh = new Mesh(2, 12, 16, 8, 3);
270
283
300
309 }
310 else if (type == 4)
311 {
312 mesh = new Mesh(2, 8, 4, 8, 3);
313
322
327
336 }
337 else
338 {
339 MFEM_ABORT("Unrecognized mesh type " << type << "!");
340 }
341 mesh->FinalizeTopology();
342
343 return mesh;
344}
345
346// Transforms the four-sided loop into a curved cylinder with skewed top and
347// base.
348void trans(const Vector &x, Vector &r)
349{
350 r.SetSize(3);
351
352 real_t tol = 1e-6;
353 real_t theta = 0.0;
354 if (fabs(x[1] + 1.0) < tol)
355 {
356 theta = 0.25 * M_PI * (x[0] - 2.0);
357 }
358 else if (fabs(x[0] - 1.0) < tol)
359 {
360 theta = 0.25 * M_PI * x[1];
361 }
362 else if (fabs(x[1] - 1.0) < tol)
363 {
364 theta = 0.25 * M_PI * (2.0 - x[0]);
365 }
366 else if (fabs(x[0] + 1.0) < tol)
367 {
368 theta = 0.25 * M_PI * (4.0 - x[1]);
369 }
370 else
371 {
372 cerr << "side not recognized "
373 << x[0] << " " << x[1] << " " << x[2] << endl;
374 }
375
376 r[0] = cos(theta);
377 r[1] = sin(theta);
378 r[2] = 0.25 * (2.0 * x[2] - 1.0) * (r[0] + 2.0);
379}
380
381// Anisotropic diffusion coefficient
383{
384 s.SetSize(3);
385 real_t a = 17.0 - 2.0 * x[0] * (1.0 + x[0]);
386 s(0,0) = 0.5 + x[0] * x[0] * (8.0 / a - 0.5);
387 s(0,1) = x[0] * x[1] * (8.0 / a - 0.5);
388 s(0,2) = 0.0;
389 s(1,0) = s(0,1);
390 s(1,1) = 0.5 * x[0] * x[0] + 8.0 * x[1] * x[1] / a;
391 s(1,2) = 0.0;
392 s(2,0) = 0.0;
393 s(2,1) = 0.0;
394 s(2,2) = a / 32.0;
395}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
int Size() const
Return the logical size of the array.
Definition array.hpp:144
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.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
Class for domain integration .
Definition lininteg.hpp:109
A general function coefficient.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
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
A matrix coefficient with an optional scalar coefficient multiplier q. The matrix function can either...
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
Adds a quadrilateral to the mesh given by 4 vertices v1 through v4.
Definition mesh.cpp:1743
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
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition mesh.cpp:3135
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
int AddBdrSegment(int v1, int v2, int attr=1)
Definition mesh.cpp:2035
void Transform(void(*f)(const Vector &, Vector &))
Definition mesh.cpp:12821
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
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).
Pointer to an Operator of a specified type.
Definition handle.hpp:34
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
Definition handle.hpp:104
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Reconstruct a solution vector x (e.g. a GridFunction) from the solution X of a constrained linear sys...
Definition operator.cpp:148
void ParseCheck(std::ostream &out=mfem::out)
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
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
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &flux, bool wcoef=true, int subdomain=-1) override
Class for parallel linear form.
Class for parallel meshes.
Definition pmesh.hpp:34
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1) override
Set the curvature of the mesh nodes using the given polynomial degree.
Definition pmesh.cpp:2007
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
real_t sigma(const Vector &x)
Definition maxwell.cpp:102
int dim
Definition ex24.cpp:53
void sigmaFunc(const Vector &x, DenseMatrix &s)
Definition ex29p.cpp:382
void duExact(const Vector &x, Vector &du)
Definition ex29p.cpp:42
void fluxExact(const Vector &x, Vector &f)
Definition ex29p.cpp:50
real_t uExact(const Vector &x)
Definition ex29p.cpp:37
Mesh * GetMesh(int type)
Definition ex29p.cpp:263
void trans(const Vector &x, Vector &r)
Definition ex29p.cpp:348
int main()
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
const char vishost[]
RefCoord s[3]