MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex1.cpp
Go to the documentation of this file.
1 // MFEM Example 1 - High-Performance Version
2 //
3 // Compile with: make ex1
4 //
5 // Sample runs: ex1 -m ../../data/fichera.mesh -perf -mf -pc lor
6 // ex1 -m ../../data/fichera.mesh -perf -asm -pc ho
7 // ex1 -m ../../data/fichera.mesh -perf -asm -pc ho -sc
8 // ex1 -m ../../data/fichera.mesh -std -asm -pc ho
9 // ex1 -m ../../data/fichera.mesh -std -asm -pc ho -sc
10 // ex1 -m ../../data/amr-hex.mesh -perf -asm -pc ho -sc
11 // ex1 -m ../../data/amr-hex.mesh -std -asm -pc ho -sc
12 // ex1 -m ../../data/ball-nurbs.mesh -perf -asm -pc ho -sc
13 // ex1 -m ../../data/ball-nurbs.mesh -std -asm -pc ho -sc
14 // ex1 -m ../../data/pipe-nurbs.mesh -perf -mf -pc lor
15 // ex1 -m ../../data/pipe-nurbs.mesh -std -asm -pc ho -sc
16 // ex1 -m ../../data/star.mesh -perf -mf -pc lor
17 // ex1 -m ../../data/star.mesh -perf -asm -pc ho
18 // ex1 -m ../../data/star.mesh -perf -asm -pc ho -sc
19 // ex1 -m ../../data/star.mesh -std -asm -pc ho
20 // ex1 -m ../../data/star.mesh -std -asm -pc ho -sc
21 // ex1 -m ../../data/amr-quad.mesh -perf -asm -pc ho -sc
22 // ex1 -m ../../data/amr-quad.mesh -std -asm -pc ho -sc
23 // ex1 -m ../../data/disc-nurbs.mesh -perf -asm -pc ho -sc
24 // ex1 -m ../../data/disc-nurbs.mesh -std -asm -pc ho -sc
25 //
26 // Description: This example code demonstrates the use of MFEM to define a
27 // simple finite element discretization of the Laplace problem
28 // -Delta u = 1 with homogeneous Dirichlet boundary conditions.
29 // Specifically, we discretize using a FE space of the specified
30 // order, or if order < 1 using an isoparametric/isogeometric
31 // space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
32 // NURBS mesh, etc.)
33 //
34 // The example highlights the use of mesh refinement, finite
35 // element grid functions, as well as linear and bilinear forms
36 // corresponding to the left-hand side and right-hand side of the
37 // discrete linear system. We also cover the explicit elimination
38 // of essential boundary conditions, static condensation, and the
39 // optional connection to the GLVis tool for visualization.
40 
41 #include "mfem-performance.hpp"
42 #include <fstream>
43 #include <iostream>
44 
45 using namespace std;
46 using namespace mfem;
47 
48 enum class PCType { NONE, LOR, HO };
49 
50 // Define template parameters for optimized build.
51 template <int dim> struct geom_t { };
52 template <>
53 struct geom_t<2> { static const Geometry::Type value = Geometry::SQUARE; };
54 template <>
55 struct geom_t<3> { static const Geometry::Type value = Geometry::CUBE; };
56 
57 const int mesh_p = 3; // mesh curvature (default: 3)
58 const int sol_p = 3; // solution order (default: 3)
59 
60 template <int dim>
61 struct ex1_t
62 {
63  static const Geometry::Type geom = geom_t<dim>::value;
64  static const int rdim = Geometry::Constants<geom>::Dimension;
65  static const int ir_order = 2*sol_p+rdim-1;
66 
67  // Static mesh type
68  using mesh_fe_t = H1_FiniteElement<geom,mesh_p>;
69  using mesh_fes_t = H1_FiniteElementSpace<mesh_fe_t>;
70  using mesh_t = TMesh<mesh_fes_t>;
71 
72  // Static solution finite element space type
73  using sol_fe_t = H1_FiniteElement<geom,sol_p>;
74  using sol_fes_t = H1_FiniteElementSpace<sol_fe_t>;
75 
76  // Static quadrature, coefficient and integrator types
77  using int_rule_t = TIntegrationRule<geom,ir_order>;
78  using coeff_t = TConstantCoefficient<>;
80 
82 
83  static int run(Mesh *mesh, int ref_levels, int order, int basis,
84  bool static_cond, PCType pc_choice, bool perf,
85  bool matrix_free, bool visualization);
86 };
87 
88 int main(int argc, char *argv[])
89 {
90  // 1. Parse command-line options.
91  const char *mesh_file = "../../data/fichera.mesh";
92  int ref_levels = -1;
93  int order = sol_p;
94  const char *basis_type = "G"; // Gauss-Lobatto
95  bool static_cond = false;
96  const char *pc = "none";
97  bool perf = true;
98  bool matrix_free = true;
99  bool visualization = 1;
100 
101  OptionsParser args(argc, argv);
102  args.AddOption(&mesh_file, "-m", "--mesh",
103  "Mesh file to use.");
104  args.AddOption(&ref_levels, "-r", "--refine",
105  "Number of times to refine the mesh uniformly;"
106  " -1 = auto: <= 50,000 elements.");
107  args.AddOption(&order, "-o", "--order",
108  "Finite element order (polynomial degree) or -1 for"
109  " isoparametric space.");
110  args.AddOption(&basis_type, "-b", "--basis-type",
111  "Basis: G - Gauss-Lobatto, P - Positive, U - Uniform");
112  args.AddOption(&perf, "-perf", "--hpc-version", "-std", "--standard-version",
113  "Enable high-performance, tensor-based, assembly/evaluation.");
114  args.AddOption(&matrix_free, "-mf", "--matrix-free", "-asm", "--assembly",
115  "Use matrix-free evaluation or efficient matrix assembly in "
116  "the high-performance version.");
117  args.AddOption(&pc, "-pc", "--preconditioner",
118  "Preconditioner: lor - low-order-refined (matrix-free) GS, "
119  "ho - high-order (assembled) GS, none.");
120  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
121  "--no-static-condensation", "Enable static condensation.");
122  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
123  "--no-visualization",
124  "Enable or disable GLVis visualization.");
125  args.Parse();
126  if (!args.Good())
127  {
128  args.PrintUsage(cout);
129  return 1;
130  }
131  if (static_cond && perf && matrix_free)
132  {
133  cout << "\nStatic condensation can not be used with matrix-free"
134  " evaluation!\n" << endl;
135  return 2;
136  }
137  MFEM_VERIFY(perf || !matrix_free,
138  "--standard-version is not compatible with --matrix-free");
139  args.PrintOptions(cout);
140 
141  PCType pc_choice;
142  if (!strcmp(pc, "ho")) { pc_choice = PCType::HO; }
143  else if (!strcmp(pc, "lor")) { pc_choice = PCType::LOR; }
144  else if (!strcmp(pc, "none")) { pc_choice = PCType::NONE; }
145  else
146  {
147  mfem_error("Invalid Preconditioner specified");
148  return 3;
149  }
150 
151  cout << "\nMFEM SIMD width: " << MFEM_SIMD_BYTES/sizeof(double)
152  << " doubles\n" << endl;
153 
154  // See class BasisType in fem/fe_coll.hpp for available basis types
155  int basis = BasisType::GetType(basis_type[0]);
156  cout << "Using " << BasisType::Name(basis) << " basis ..." << endl;
157 
158  // 2. Read the mesh from the given mesh file. We can handle triangular,
159  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
160  // the same code.
161  Mesh *mesh = new Mesh(mesh_file, 1, 1);
162  int dim = mesh->Dimension();
163 
164  if (dim == 2)
165  {
166  return ex1_t<2>::run(mesh, ref_levels, order, basis, static_cond,
167  pc_choice, perf, matrix_free, visualization);
168  }
169  else if (dim == 3)
170  {
171  return ex1_t<3>::run(mesh, ref_levels, order, basis, static_cond,
172  pc_choice, perf, matrix_free, visualization);
173  }
174  else
175  {
176  MFEM_ABORT("Dimension must be 2 or 3.")
177  }
178 
179  return 0;
180 }
181 
182 template <int dim>
183 int ex1_t<dim>::run(Mesh *mesh, int ref_levels, int order, int basis,
184  bool static_cond, PCType pc_choice, bool perf,
185  bool matrix_free, bool visualization)
186 {
187  // 3. Check if the optimized version matches the given mesh
188  if (perf)
189  {
190  cout << "High-performance version using integration rule with "
191  << int_rule_t::qpts << " points ..." << endl;
192  if (!mesh_t::MatchesGeometry(*mesh))
193  {
194  cout << "The given mesh does not match the optimized 'geom' parameter.\n"
195  << "Recompile with suitable 'geom' value." << endl;
196  delete mesh;
197  return 4;
198  }
199  else if (!mesh_t::MatchesNodes(*mesh))
200  {
201  cout << "Switching the mesh curvature to match the "
202  << "optimized value (order " << mesh_p << ") ..." << endl;
203  mesh->SetCurvature(mesh_p, false, -1, Ordering::byNODES);
204  }
205  }
206 
207  // 4. Refine the mesh to increase the resolution. In this example we do
208  // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
209  // largest number that gives a final mesh with no more than 50,000
210  // elements, or as specified on the command line with the option
211  // '--refine'.
212  {
213  ref_levels = (ref_levels != -1) ? ref_levels :
214  (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
215  for (int l = 0; l < ref_levels; l++)
216  {
217  mesh->UniformRefinement();
218  }
219  }
220  if (mesh->MeshGenerator() & 1) // simplex mesh
221  {
222  MFEM_VERIFY(pc_choice != PCType::LOR, "triangle and tet meshes do not "
223  " support the LOR preconditioner yet");
224  }
225 
226  // 5. Define a finite element space on the mesh. Here we use continuous
227  // Lagrange finite elements of the specified order. If order < 1, we
228  // instead use an isoparametric/isogeometric space.
230  if (order > 0)
231  {
232  fec = new H1_FECollection(order, dim, basis);
233  }
234  else if (mesh->GetNodes())
235  {
236  fec = mesh->GetNodes()->OwnFEC();
237  cout << "Using isoparametric FEs: " << fec->Name() << endl;
238  }
239  else
240  {
241  fec = new H1_FECollection(order = 1, dim, basis);
242  }
243  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
244  cout << "Number of finite element unknowns: "
245  << fespace->GetTrueVSize() << endl;
246 
247  // Create the LOR mesh and finite element space. In the settings of this
248  // example, we can transfer between HO and LOR with the identity operator.
249  Mesh mesh_lor;
250  FiniteElementCollection *fec_lor = NULL;
251  FiniteElementSpace *fespace_lor = NULL;
252  if (pc_choice == PCType::LOR)
253  {
254  int basis_lor = basis;
255  if (basis == BasisType::Positive) { basis_lor=BasisType::ClosedUniform; }
256  mesh_lor = Mesh::MakeRefined(*mesh, order, basis_lor);
257  fec_lor = new H1_FECollection(1, dim);
258  fespace_lor = new FiniteElementSpace(&mesh_lor, fec_lor);
259  }
260 
261  // 6. Check if the optimized version matches the given space
262  if (perf && !sol_fes_t::Matches(*fespace))
263  {
264  cout << "The given order does not match the optimized parameter.\n"
265  << "Recompile with suitable 'sol_p' value." << endl;
266  delete fespace;
267  delete fec;
268  delete mesh;
269  return 5;
270  }
271 
272  // 7. Determine the list of true (i.e. conforming) essential boundary dofs.
273  // In this example, the boundary conditions are defined by marking all
274  // the boundary attributes from the mesh as essential (Dirichlet) and
275  // converting them to a list of true dofs.
276  Array<int> ess_tdof_list;
277  if (mesh->bdr_attributes.Size())
278  {
279  Array<int> ess_bdr(mesh->bdr_attributes.Max());
280  ess_bdr = 1;
281  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
282  }
283 
284  // 8. Set up the linear form b(.) which corresponds to the right-hand side of
285  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
286  // the basis functions in the finite element fespace.
287  LinearForm *b = new LinearForm(fespace);
288  ConstantCoefficient one(1.0);
290  b->Assemble();
291 
292  // 9. Define the solution vector x as a finite element grid function
293  // corresponding to fespace. Initialize x with initial guess of zero,
294  // which satisfies the boundary conditions.
295  GridFunction x(fespace);
296  x = 0.0;
297 
298  // 10. Set up the bilinear form a(.,.) on the finite element space that will
299  // hold the matrix corresponding to the Laplacian operator -Delta.
300  // Optionally setup a form to be assembled for preconditioning (a_pc).
301  BilinearForm *a = new BilinearForm(fespace);
302  BilinearForm *a_pc = NULL;
303  if (pc_choice == PCType::LOR) { a_pc = new BilinearForm(fespace_lor); }
304  if (pc_choice == PCType::HO) { a_pc = new BilinearForm(fespace); }
305 
306  // 11. Assemble the bilinear form and the corresponding linear system,
307  // applying any necessary transformations such as: eliminating boundary
308  // conditions, applying conforming constraints for non-conforming AMR,
309  // static condensation, etc.
310  if (static_cond)
311  {
313  MFEM_VERIFY(pc_choice != PCType::LOR,
314  "cannot use LOR preconditioner with static condensation");
315  }
316 
317  cout << "Assembling the bilinear form ..." << flush;
318  tic_toc.Clear();
319  tic_toc.Start();
320  // Pre-allocate sparsity assuming dense element matrices
322 
323  HPCBilinearForm *a_hpc = NULL;
324  Operator *a_oper = NULL;
325 
326  if (!perf)
327  {
328  // Standard assembly using a diffusion domain integrator
330  a->Assemble();
331  }
332  else
333  {
334  // High-performance assembly/evaluation using the templated operator type
335  a_hpc = new HPCBilinearForm(integ_t(coeff_t(1.0)), *fespace);
336  if (matrix_free)
337  {
338  a_hpc->Assemble(); // partial assembly
339  }
340  else
341  {
342  a_hpc->AssembleBilinearForm(*a); // full matrix assembly
343  }
344  }
345  tic_toc.Stop();
346  cout << " done, " << tic_toc.RealTime() << "s." << endl;
347 
348  // 12. Solve the system A X = B with CG. In the standard case, use a simple
349  // symmetric Gauss-Seidel preconditioner.
350 
351  // Setup the operator matrix (if applicable)
352  SparseMatrix A;
353  Vector B, X;
354  if (perf && matrix_free)
355  {
356  a_hpc->FormLinearSystem(ess_tdof_list, x, *b, a_oper, X, B);
357  cout << "Size of linear system: " << a_hpc->Height() << endl;
358  }
359  else
360  {
361  a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
362  cout << "Size of linear system: " << A.Height() << endl;
363  a_oper = &A;
364  }
365 
366  // Setup the matrix used for preconditioning
367  cout << "Assembling the preconditioning matrix ..." << flush;
368  tic_toc.Clear();
369  tic_toc.Start();
370 
371  SparseMatrix A_pc;
372  if (pc_choice == PCType::LOR)
373  {
374  // TODO: assemble the LOR matrix using the performance code
376  a_pc->UsePrecomputedSparsity();
377  a_pc->Assemble();
378  a_pc->FormSystemMatrix(ess_tdof_list, A_pc);
379  }
380  else if (pc_choice == PCType::HO)
381  {
382  if (!matrix_free)
383  {
384  A_pc.MakeRef(A); // matrix already assembled, reuse it
385  }
386  else
387  {
388  a_pc->UsePrecomputedSparsity();
389  a_hpc->AssembleBilinearForm(*a_pc);
390  a_pc->FormSystemMatrix(ess_tdof_list, A_pc);
391  }
392  }
393 
394  tic_toc.Stop();
395  cout << " done, " << tic_toc.RealTime() << "s." << endl;
396 
397  // Solve with CG or PCG, depending if the matrix A_pc is available
398  if (pc_choice != PCType::NONE)
399  {
400  GSSmoother M(A_pc);
401  PCG(*a_oper, M, B, X, 1, 500, 1e-12, 0.0);
402  }
403  else
404  {
405  CG(*a_oper, B, X, 1, 500, 1e-12, 0.0);
406  }
407 
408  // 13. Recover the solution as a finite element grid function.
409  if (perf && matrix_free)
410  {
411  a_hpc->RecoverFEMSolution(X, *b, x);
412  }
413  else
414  {
415  a->RecoverFEMSolution(X, *b, x);
416  }
417 
418  // 14. Save the refined mesh and the solution. This output can be viewed later
419  // using GLVis: "glvis -m refined.mesh -g sol.gf".
420  ofstream mesh_ofs("refined.mesh");
421  mesh_ofs.precision(8);
422  mesh->Print(mesh_ofs);
423  ofstream sol_ofs("sol.gf");
424  sol_ofs.precision(8);
425  x.Save(sol_ofs);
426 
427  // 15. Send the solution by socket to a GLVis server.
428  if (visualization)
429  {
430  char vishost[] = "localhost";
431  int visport = 19916;
432  socketstream sol_sock(vishost, visport);
433  sol_sock.precision(8);
434  sol_sock << "solution\n" << *mesh << x << flush;
435  }
436 
437  // 16. Free the used memory.
438  delete a;
439  delete a_hpc;
440  if (a_oper != &A) { delete a_oper; }
441  delete a_pc;
442  delete b;
443  delete fespace;
444  delete fespace_lor;
445  delete fec_lor;
446  if (order > 0) { delete fec; }
447  delete mesh;
448 
449  return 0;
450 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
Templated bilinear form class, cf. bilinearform.?pp.
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
void MakeRef(const SparseMatrix &master)
Clear the contents of the SparseMatrix and make it a reference to master.
Definition: sparsemat.cpp:305
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:102
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
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:564
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
StopWatch tic_toc
Definition: tic_toc.cpp:447
double RealTime()
Definition: tic_toc.cpp:426
PCType
Definition: ex1.cpp:48
The Integrator class combines a kernel and a coefficient.
Definition: tbilininteg.hpp:26
Data type for Gauss-Seidel smoother of sparse matrix.
void Stop()
Stop the stopwatch.
Definition: tic_toc.cpp:416
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
Data type sparse matrix.
Definition: sparsemat.hpp:46
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
constexpr char vishost[]
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
constexpr int visport
const int sol_p
Definition: ex1.cpp:58
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:5312
int MeshGenerator()
Get the mesh generator/type.
Definition: mesh.hpp:921
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
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:889
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
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:904
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:566
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition: tic_toc.cpp:411
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
void UsePrecomputedSparsity(int ps=1)
For scalar FE spaces, precompute the sparsity pattern of the matrix (assuming dense element matrices)...
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
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
FDualNumber< tbase > log(const FDualNumber< tbase > &f)
log([dual number])
Definition: fdual.hpp:515
double a
Definition: lissajous.cpp:41
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
virtual const char * Name() const
Definition: fe_coll.hpp:61
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1646
const int mesh_p
Definition: ex1.cpp:57
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7841
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:23
Abstract operator.
Definition: operator.hpp:24
int main()
void Clear()
Clear the elapsed time on the stopwatch and restart it if it&#39;s running.
Definition: tic_toc.cpp:406
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150