MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex39p.cpp
Go to the documentation of this file.
1// MFEM Example 39 - Parallel Version
2//
3// Compile with: make ex39p
4//
5// Sample runs: mpirun -np 4 ex39p
6// mpirun -np 4 ex39p -ess "Southern Boundary"
7// mpirun -np 4 ex39p -src Base
8//
9// Description: This example code demonstrates the use of named attribute
10// sets in MFEM to specify material regions, boundary regions,
11// or source regions by name rather than attribute numbers. It
12// also demonstrates how new named attribute sets may be created
13// from arbitrary groupings of attribute numbers and used as a
14// convenient shorthand to refer to those groupings in other
15// portions of the application or through the command line.
16//
17// The particular problem being solved here is nearly the same
18// as that in example 1 i.e. a simple finite element
19// discretization of the Laplace problem -Delta u = 1 with
20// homogeneous Dirichlet boundary conditions and, in this case,
21// an inhomogeneous diffusion coefficient. The diffusion
22// coefficient is given a small default value throughout the
23// domain which is increased by two separate amounts in two named
24// regions.
25//
26// This example makes use of a specific input mesh, "compass.msh",
27// containing named domain and boundary regions generated by Gmsh
28// and stored in their "msh" format (version 2.2). This file
29// defines eight boundary regions corresponding to eight compass
30// headings; "ENE", "NNE", "NNW", "WSW", "SSW", "SSE", and "ESE".
31// It also defines nine domain regions; "Base", "N Even", "N Odd",
32// "W Even", "W Odd", "S Even", "S Odd", "E Even", and "E Odd".
33// These regions split the four compass pointers into two halves
34// each and also label the remaining elements as "Base". Starting
35// with these named regions we test the construction of named
36// sets as well as reading and writing these named groupings from
37// and to mesh files.
38//
39// The example highlights the use of named attribute sets for
40// both subdomains and boundaries in different contexts as well
41// as basic methods to create named sets from existing attributes.
42
43#include "mfem.hpp"
44#include <fstream>
45#include <iostream>
46
47using namespace std;
48using namespace mfem;
49
50int main(int argc, char *argv[])
51{
52 // 1. Initialize MPI and HYPRE.
53 Mpi::Init();
55
56 // 2. Parse command-line options.
57 const char *mesh_file = "../data/compass.msh";
58 int order = 1;
59 string source_name = "Rose Even";
60 string ess_name = "Boundary";
61 bool visualization = true;
62
63 OptionsParser args(argc, argv);
65 "Mesh file to use.");
67 "Finite element order (polynomial degree) or -1 for"
68 " isoparametric space.");
70 "Name of attribute set containing source.");
72 "Name of attribute set containing essential BC.");
74 "--no-visualization",
75 "Enable or disable GLVis visualization.");
76 args.ParseCheck();
77
78 // 3. Read the serial mesh from the given mesh file.
79 Mesh mesh(mesh_file, 1, 1);
80 int dim = mesh.Dimension();
81
82 // 4. Refine the serial mesh on all processors to increase the resolution. In
83 // this example we do 'ref_levels' of uniform refinement. We choose
84 // 'ref_levels' to be the largest number that gives a final mesh with no
85 // more than 10,000 elements.
86 {
87 int ref_levels =
88 (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
89 for (int l = 0; l < ref_levels; l++)
90 {
91 mesh.UniformRefinement();
92 }
93 }
94
95 // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
96 // this mesh further in parallel to increase the resolution. Once the
97 // parallel mesh is defined, the serial mesh can be deleted.
98 ParMesh pmesh(MPI_COMM_WORLD, mesh);
99 mesh.Clear();
100 {
101 int par_ref_levels = 2;
102 for (int l = 0; l < par_ref_levels; l++)
103 {
104 pmesh.UniformRefinement();
105 }
106 }
107
108 // 6a. Display attribute set names contained in the initial mesh
109 AttributeSets &attr_sets = pmesh.attribute_sets;
110 AttributeSets &bdr_attr_sets = pmesh.bdr_attribute_sets;
111 if (Mpi::Root())
112 {
113 std::set<string> names = attr_sets.GetAttributeSetNames();
114 cout << "Element Attribute Set Names: ";
115 for (auto const &set_name : names)
116 {
117 cout << " \"" << set_name << "\"";
118 }
119 cout << endl;
120
121 std::set<string> bdr_names = bdr_attr_sets.GetAttributeSetNames();
122 cout << "Boundary Attribute Set Names: ";
123 for (auto const &bdr_set_name : bdr_names)
124 {
125 cout << " \"" << bdr_set_name << "\"";
126 }
127 cout << endl;
128 }
129
130 // 6b. Define new regions based on existing attribute sets
131 {
132 Array<int> & Na = attr_sets.GetAttributeSet("N Even");
133 Array<int> & Nb = attr_sets.GetAttributeSet("N Odd");
134 Array<int> & Sa = attr_sets.GetAttributeSet("S Even");
135 Array<int> & Sb = attr_sets.GetAttributeSet("S Odd");
136 Array<int> & Ea = attr_sets.GetAttributeSet("E Even");
137 Array<int> & Eb = attr_sets.GetAttributeSet("E Odd");
138 Array<int> & Wa = attr_sets.GetAttributeSet("W Even");
139 Array<int> & Wb = attr_sets.GetAttributeSet("W Odd");
140
141 // Create a new set spanning the North point
142 attr_sets.SetAttributeSet("North", Na);
144
145 // Create a new set spanning the South point
146 attr_sets.SetAttributeSet("South", Sa);
148
149 // Create a new set spanning the East point
150 attr_sets.SetAttributeSet("East", Ea);
152
153 // Create a new set spanning the West point
154 attr_sets.SetAttributeSet("West", Wa);
156
157 // Create a new set consisting of the "a" sides of the compass rose
158 attr_sets.SetAttributeSet("Rose Even", Na);
162
163 // Create a new set consisting of the "b" sides of the compass rose
164 attr_sets.SetAttributeSet("Rose Odd", Nb);
168
169
170 // Create a new set consisting of the full compass rose
171 Array<int> & Ra = attr_sets.GetAttributeSet("Rose Even");
172 Array<int> & Rb = attr_sets.GetAttributeSet("Rose Odd");
173 attr_sets.SetAttributeSet("Rose", Ra);
175 }
176 // 6c. Define new boundary regions based on existing boundary attribute sets
177 {
178 Array<int> & NNE = bdr_attr_sets.GetAttributeSet("NNE");
179 Array<int> & NNW = bdr_attr_sets.GetAttributeSet("NNW");
180 Array<int> & ENE = bdr_attr_sets.GetAttributeSet("ENE");
181 Array<int> & ESE = bdr_attr_sets.GetAttributeSet("ESE");
182 Array<int> & SSE = bdr_attr_sets.GetAttributeSet("SSE");
183 Array<int> & SSW = bdr_attr_sets.GetAttributeSet("SSW");
184 Array<int> & WNW = bdr_attr_sets.GetAttributeSet("WNW");
185 Array<int> & WSW = bdr_attr_sets.GetAttributeSet("WSW");
186
187 bdr_attr_sets.SetAttributeSet("Northern Boundary", NNE);
189
190 bdr_attr_sets.SetAttributeSet("Southern Boundary", SSE);
192
193 bdr_attr_sets.SetAttributeSet("Eastern Boundary", ENE);
195
196 bdr_attr_sets.SetAttributeSet("Western Boundary", WNW);
198
199 bdr_attr_sets.SetAttributeSet("Boundary",
200 bdr_attr_sets.GetAttributeSet
201 ("Northern Boundary"));
203 bdr_attr_sets.GetAttributeSet
204 ("Southern Boundary"));
206 bdr_attr_sets.GetAttributeSet
207 ("Eastern Boundary"));
209 bdr_attr_sets.GetAttributeSet
210 ("Western Boundary"));
211 }
212
213 // 7. Define a parallel finite element space on the parallel mesh. Here we
214 // use continuous Lagrange finite elements of the specified order. If
215 // order < 1, we instead use an isoparametric/isogeometric space.
216 H1_FECollection fec(order, dim);
217 ParFiniteElementSpace fespace(&pmesh, &fec);
218 HYPRE_BigInt size = fespace.GlobalTrueVSize();
219 if (Mpi::Root())
220 {
221 cout << "Number of finite element unknowns: " << size << endl;
222 }
223
224 // 8. Determine the list of true (i.e. parallel conforming) essential
225 // boundary dofs. In this example, the boundary conditions are defined
226 // by marking all the boundary regions corresponding to the boundary
227 // attributes contained in the set named "ess_name" as essential
228 // (Dirichlet) and converting them to a list of true dofs.
229 Array<int> ess_tdof_list;
230 if (bdr_attr_sets.AttributeSetExists(ess_name))
231 {
232 Array<int> ess_bdr_marker = bdr_attr_sets.GetAttributeSetMarker(ess_name);
233 fespace.GetEssentialTrueDofs(ess_bdr_marker, ess_tdof_list);
234 }
235
236 // 9. Set up the parallel linear form b(.) which corresponds to the
237 // right-hand side of the FEM linear system, which in this case is
238 // (1_s,phi_i) where phi_i are the basis functions in fespace and 1_s
239 // is an indicator function equal to 1 on the region defined by the
240 // named set "source_name" and zero elsewhere.
241 Array<int> source_marker = attr_sets.GetAttributeSetMarker(source_name);
242
243 ParLinearForm b(&fespace);
244 ConstantCoefficient one(1.0);
246 b.Assemble();
247
248 // 10. Define the solution vector x as a parallel finite element grid
249 // function corresponding to fespace. Initialize x with initial guess of
250 // zero, which satisfies the boundary conditions.
251 ParGridFunction x(&fespace);
252 x = 0.0;
253
254 // 11. Set up the parallel bilinear form a(.,.) on the finite element space
255 // corresponding to the Laplacian operator -Delta, by adding the
256 // Diffusion domain integrator.
257 ParBilinearForm a(&fespace);
258
259 ConstantCoefficient defaultCoef(1.0e-6);
260 ConstantCoefficient baseCoef(1.0);
261 ConstantCoefficient roseCoef(2.0);
262
263 Array<int> base_marker = attr_sets.GetAttributeSetMarker("Base");
264 Array<int> rose_marker = attr_sets.GetAttributeSetMarker("Rose Even");
265
266 // Impose a very small diffusion coefficient across the entire mesh
268
269 // Impose an additional, stronger diffusion coefficient in select regions
272
273 // 12. Assemble the parallel bilinear form and the corresponding linear
274 // system, applying any necessary transformations.
275 a.Assemble();
276
278 Vector B, X;
279 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
280
281 // 13. Solve the system using PCG with hypre's BoomerAMG preconditioner.
282 HypreBoomerAMG M(A);
283 CGSolver cg(MPI_COMM_WORLD);
284 cg.SetRelTol(1e-12);
285 cg.SetMaxIter(2000);
286 cg.SetPrintLevel(1);
287 cg.SetPreconditioner(M);
288 cg.SetOperator(A);
289 cg.Mult(B, X);
290
291 // 14. Recover the parallel grid function corresponding to X. This is the
292 // local finite element solution on each processor.
293 a.RecoverFEMSolution(X, b, x);
294
295 // 15. Save the refined mesh and the solution in parallel. This output can
296 // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
297 pmesh.Save("mesh");
298 x.Save("sol");
299
300 // 16. Send the solution by socket to a GLVis server.
301 if (visualization)
302 {
303 char vishost[] = "localhost";
304 int visport = 19916;
305 int num_procs = Mpi::WorldSize();
306 int myid = Mpi::WorldRank();
307 socketstream sol_sock(vishost, visport);
308 sol_sock << "parallel " << num_procs << " " << myid << "\n";
309 sol_sock.precision(8);
310 sol_sock << "solution\n" << pmesh << x << "keys Rjmm" << flush;
311 }
312
313 return 0;
314}
void SetAttributeSet(const std::string &set_name, const Array< int > &attr)
Create a new attribute set.
Array< int > GetAttributeSetMarker(const std::string &set_name)
Return a marker array corresponding to a named attribute set.
bool AttributeSetExists(const std::string &name) const
Return true is the named attribute set is present.
Array< int > & GetAttributeSet(const std::string &set_name)
Access a named attribute set.
void AddToAttributeSet(const std::string &set_name, int attr)
Add a single entry to an existing attribute set.
std::set< std::string > GetAttributeSetNames() const
Return all attribute set names as an STL set.
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 domain integration .
Definition lininteg.hpp:109
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
Mesh data type.
Definition mesh.hpp:56
AttributeSets bdr_attribute_sets
Named sets of boundary element attributes.
Definition mesh.hpp:288
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
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
AttributeSets attribute_sets
Named sets of element attributes.
Definition mesh.hpp:285
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).
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
Class for parallel linear form.
Class for parallel meshes.
Definition pmesh.hpp:34
void Save(const std::string &fname, int precision=16) const override
Definition pmesh.cpp:4940
Vector data type.
Definition vector.hpp:80
int dim
Definition ex24.cpp:53
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
const char vishost[]