MFEM  v4.6.0
Finite element discretization library
ex1.cpp
Go to the documentation of this file.
1 // MFEM Example 1
2 // Caliper Modification
3 //
4 // Compile with: make ex1
5 //
6 // Sample runs: ex1 -m ../data/square-disc.mesh
7 // ex1 -m ../data/star.mesh
8 // ex1 -m ../data/star-mixed.mesh
9 // ex1 -m ../data/escher.mesh
10 // ex1 -m ../data/fichera.mesh
11 // ex1 -m ../data/fichera-mixed.mesh
12 // ex1 -m ../data/toroid-wedge.mesh
13 // ex1 -m ../data/periodic-annulus-sector.msh
14 // ex1 -m ../data/periodic-torus-sector.msh
15 // ex1 -m ../data/square-disc-p2.vtk -o 2
16 // ex1 -m ../data/square-disc-p3.mesh -o 3
17 // ex1 -m ../data/square-disc-nurbs.mesh -o -1
18 // ex1 -m ../data/star-mixed-p2.mesh -o 2
19 // ex1 -m ../data/disc-nurbs.mesh -o -1
20 // ex1 -m ../data/pipe-nurbs.mesh -o -1
21 // ex1 -m ../data/fichera-mixed-p2.mesh -o 2
22 // ex1 -m ../data/star-surf.mesh
23 // ex1 -m ../data/square-disc-surf.mesh
24 // ex1 -m ../data/inline-segment.mesh
25 // ex1 -m ../data/amr-quad.mesh
26 // ex1 -m ../data/amr-hex.mesh
27 // ex1 -m ../data/fichera-amr.mesh
28 // ex1 -m ../data/mobius-strip.mesh
29 // ex1 -m ../data/mobius-strip.mesh -o -1 -sc
30 //
31 // Device sample runs:
32 // ex1 -pa -d cuda
33 // ex1 -pa -d raja-cuda
34 // ex1 -pa -d occa-cuda
35 // ex1 -pa -d raja-omp
36 // ex1 -pa -d occa-omp
37 // ex1 -pa -d ceed-cpu
38 // * ex1 -pa -d ceed-cuda
39 // ex1 -pa -d ceed-cuda:/gpu/cuda/shared
40 // ex1 -m ../data/beam-hex.mesh -pa -d cuda
41 // ex1 -m ../data/beam-tet.mesh -pa -d ceed-cpu
42 // ex1 -m ../data/beam-tet.mesh -pa -d ceed-cuda:/gpu/cuda/ref
43 //
44 // Description: This example is a copy of Example 1 instrumented with the
45 // Caliper performance profilinh library. Any option supported by
46 // the Caliper ConfigManager can be passed to the code using a
47 // configuration string after -p or --caliper flag. For more
48 // information, see the Caliper documentation.
49 //
50 // Examples: ex1 --caliper runtime-report
51 // ex1 --caliper runtime-report,mem.highwatermark
52 //
53 // The first run will return the default report. The second run will also output
54 // the memory high-water mark and time spent in MPI routines.
55 
56 #include "mfem.hpp"
57 #include <fstream>
58 #include <iostream>
59 
60 using namespace std;
61 using namespace mfem;
62 
63 int main(int argc, char *argv[])
64 {
65  // Define Caliper ConfigManager
66  cali::ConfigManager mgr;
67  // Caliper instrumentation
68  MFEM_PERF_FUNCTION;
69 
70  // 1. Parse command-line options.
71  const char *mesh_file = "../../data/star.mesh";
72  int order = 1;
73  bool static_cond = false;
74  bool pa = false;
75  const char *device_config = "cpu";
76  bool visualization = true;
77  const char* cali_config = "runtime-report";
78 
79  OptionsParser args(argc, argv);
80  args.AddOption(&mesh_file, "-m", "--mesh",
81  "Mesh file to use.");
82  args.AddOption(&order, "-o", "--order",
83  "Finite element order (polynomial degree) or -1 for"
84  " isoparametric space.");
85  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
86  "--no-static-condensation", "Enable static condensation.");
87  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
88  "--no-partial-assembly", "Enable Partial Assembly.");
89  args.AddOption(&device_config, "-d", "--device",
90  "Device configuration string, see Device::Configure().");
91  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
92  "--no-visualization",
93  "Enable or disable GLVis visualization.");
94  args.AddOption(&cali_config, "-p", "--caliper",
95  "Caliper configuration string.");
96 
97  args.Parse();
98  if (!args.Good())
99  {
100  args.PrintUsage(cout);
101  return 1;
102  }
103  args.PrintOptions(cout);
104 
105  // 2. Enable hardware devices such as GPUs, and programming models such as
106  // CUDA, OCCA, RAJA and OpenMP based on command line options.
107  Device device(device_config);
108  device.Print();
109 
110  // Caliper configuration
111  mgr.add(cali_config);
112  mgr.start();
113 
114  // 3. Read the mesh from the given mesh file. We can handle triangular,
115  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
116  // the same code.
117  Mesh mesh(mesh_file, 1, 1);
118  int dim = mesh.Dimension();
119 
120  // 4. Refine the mesh to increase the resolution. In this example we do
121  // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
122  // largest number that gives a final mesh with no more than 50,000
123  // elements.
124  {
125  int ref_levels =
126  (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
127  for (int l = 0; l < ref_levels; l++)
128  {
129  mesh.UniformRefinement();
130  }
131  }
132 
133  // 5. Define a finite element space on the mesh. Here we use continuous
134  // Lagrange finite elements of the specified order. If order < 1, we
135  // instead use an isoparametric/isogeometric space.
137  bool delete_fec;
138  if (order > 0)
139  {
140  fec = new H1_FECollection(order, dim);
141  delete_fec = true;
142  }
143  else if (mesh.GetNodes())
144  {
145  fec = mesh.GetNodes()->OwnFEC();
146  delete_fec = false;
147  cout << "Using isoparametric FEs: " << fec->Name() << endl;
148  }
149  else
150  {
151  fec = new H1_FECollection(order = 1, dim);
152  delete_fec = true;
153  }
154  FiniteElementSpace fespace(&mesh, fec);
155  cout << "Number of finite element unknowns: "
156  << fespace.GetTrueVSize() << endl;
157 
158  // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
159  // In this example, the boundary conditions are defined by marking all
160  // the boundary attributes from the mesh as essential (Dirichlet) and
161  // converting them to a list of true dofs.
162  Array<int> ess_tdof_list;
163  if (mesh.bdr_attributes.Size())
164  {
165  Array<int> ess_bdr(mesh.bdr_attributes.Max());
166  ess_bdr = 1;
167  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
168  }
169 
170  // 7. Set up the linear form b(.) which corresponds to the right-hand side of
171  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
172  // the basis functions in the finite element fespace.
173  MFEM_PERF_BEGIN("Set up the linear form");
174  LinearForm b(&fespace);
175  ConstantCoefficient one(1.0);
176  b.AddDomainIntegrator(new DomainLFIntegrator(one));
177  b.Assemble();
178  MFEM_PERF_END("Set up the linear form");
179 
180  // 8. Define the solution vector x as a finite element grid function
181  // corresponding to fespace. Initialize x with initial guess of zero,
182  // which satisfies the boundary conditions.
183  GridFunction x(&fespace);
184  x = 0.0;
185 
186  // 9. Set up the bilinear form a(.,.) on the finite element space
187  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
188  // domain integrator.
189  MFEM_PERF_BEGIN("Set up the bilinear form");
190  BilinearForm a(&fespace);
191  if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
192  a.AddDomainIntegrator(new DiffusionIntegrator(one));
193 
194  // 10. Assemble the bilinear form and the corresponding linear system,
195  // applying any necessary transformations such as: eliminating boundary
196  // conditions, applying conforming constraints for non-conforming AMR,
197  // static condensation, etc.
198  if (static_cond) { a.EnableStaticCondensation(); }
199  a.Assemble();
200 
201  OperatorPtr A;
202  Vector B, X;
203  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
204  MFEM_PERF_END("Set up the bilinear form");
205 
206  cout << "Size of linear system: " << A->Height() << endl;
207 
208  // 11. Solve the linear system A X = B.
209  if (!pa)
210  {
211  MFEM_PERF_SCOPE("Solve A X=B (FA)");
212 #ifndef MFEM_USE_SUITESPARSE
213  // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
214  GSSmoother M((SparseMatrix&)(*A));
215  PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
216 #else
217  // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
218  UMFPackSolver umf_solver;
219  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
220  umf_solver.SetOperator(*A);
221  umf_solver.Mult(B, X);
222 #endif
223  }
224  else // Jacobi preconditioning in partial assembly mode
225  {
226  MFEM_PERF_SCOPE("Solve A X=B (PA)");
227  if (UsesTensorBasis(fespace))
228  {
229  OperatorJacobiSmoother M(a, ess_tdof_list);
230  PCG(*A, M, B, X, 1, 400, 1e-12, 0.0);
231  }
232  else
233  {
234  CG(*A, B, X, 1, 400, 1e-12, 0.0);
235  }
236  }
237  // 12. Recover the solution as a finite element grid function.
238  a.RecoverFEMSolution(X, b, x);
239 
240  // 13. Save the refined mesh and the solution. This output can be viewed later
241  // using GLVis: "glvis -m refined.mesh -g sol.gf".
242  MFEM_PERF_BEGIN("Save the results");
243  ofstream mesh_ofs("refined.mesh");
244  mesh_ofs.precision(8);
245  mesh.Print(mesh_ofs);
246  ofstream sol_ofs("sol.gf");
247  sol_ofs.precision(8);
248  x.Save(sol_ofs);
249  MFEM_PERF_END("Save the results");
250  // 14. Send the solution by socket to a GLVis server.
251  if (visualization)
252  {
253  char vishost[] = "localhost";
254  int visport = 19916;
255  socketstream sol_sock(vishost, visport);
256  sol_sock.precision(8);
257  sol_sock << "solution\n" << mesh << x << flush;
258  }
259 
260  // 15. Free the used memory.
261  if (delete_fec)
262  {
263  delete fec;
264  }
265 
266  // Flush output
267  mgr.flush();
268 
269  return 0;
270 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
int visport
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition: fespace.cpp:587
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Definition: fespace.hpp:1306
int main(int argc, char *argv[])
Definition: ex1.cpp:74
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:1070
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
Data type sparse matrix.
Definition: sparsemat.hpp:50
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:302
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
virtual const char * Name() const
Definition: fe_coll.hpp:80
void CG(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Conjugate gradient method. (tolerances are squared)
Definition: solvers.cpp:898
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition: solvers.cpp:913
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:1081
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
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 &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3194
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
double a
Definition: lissajous.cpp:41
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
int dim
Definition: ex24.cpp:53
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3696
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:3099