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