MFEM  v4.1.0 Finite element discretization library
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
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);
78  "Mesh file to use.");
80  "Number of times to refine the mesh uniformly;"
81  " -1 = auto: <= 50,000 elements.");
83  "Finite element order (polynomial degree) or -1 for"
84  " isoparametric space.");
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.");
93  "Preconditioner: lor - low-order-refined (matrix-free) GS, "
94  "ho - high-order (assembled) GS, none.");
96  "--no-static-condensation", "Enable static condensation.");
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  // See class BasisType in fem/fe_coll.hpp for available basis types
128  int basis = BasisType::GetType(basis_type[0]);
129  cout << "Using " << BasisType::Name(basis) << " basis ..." << endl;
130
131  // 2. Read the mesh from the given mesh file. We can handle triangular,
132  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
133  // the same code.
134  Mesh *mesh = new Mesh(mesh_file, 1, 1);
135  int dim = mesh->Dimension();
136
137  // 3. Check if the optimized version matches the given mesh
138  if (perf)
139  {
140  cout << "High-performance version using integration rule with "
141  << int_rule_t::qpts << " points ..." << endl;
142  if (!mesh_t::MatchesGeometry(*mesh))
143  {
144  cout << "The given mesh does not match the optimized 'geom' parameter.\n"
145  << "Recompile with suitable 'geom' value." << endl;
146  delete mesh;
147  return 4;
148  }
149  else if (!mesh_t::MatchesNodes(*mesh))
150  {
151  cout << "Switching the mesh curvature to match the "
152  << "optimized value (order " << mesh_p << ") ..." << endl;
153  mesh->SetCurvature(mesh_p, false, -1, Ordering::byNODES);
154  }
155  }
156
157  // 4. Refine the mesh to increase the resolution. In this example we do
158  // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
159  // largest number that gives a final mesh with no more than 50,000
160  // elements, or as specified on the command line with the option
161  // '--refine'.
162  {
163  ref_levels = (ref_levels != -1) ? ref_levels :
164  (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
165  for (int l = 0; l < ref_levels; l++)
166  {
167  mesh->UniformRefinement();
168  }
169  }
170  if (mesh->MeshGenerator() & 1) // simplex mesh
171  {
172  MFEM_VERIFY(pc_choice != LOR, "triangle and tet meshes do not support"
173  " the LOR preconditioner yet");
174  }
175
176  // 5. Define a finite element space on the mesh. Here we use continuous
177  // Lagrange finite elements of the specified order. If order < 1, we
178  // instead use an isoparametric/isogeometric space.
180  if (order > 0)
181  {
182  fec = new H1_FECollection(order, dim, basis);
183  }
184  else if (mesh->GetNodes())
185  {
186  fec = mesh->GetNodes()->OwnFEC();
187  cout << "Using isoparametric FEs: " << fec->Name() << endl;
188  }
189  else
190  {
191  fec = new H1_FECollection(order = 1, dim, basis);
192  }
193  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
194  cout << "Number of finite element unknowns: "
195  << fespace->GetTrueVSize() << endl;
196
197  // Create the LOR mesh and finite element space. In the settings of this
198  // example, we can transfer between HO and LOR with the identity operator.
199  Mesh *mesh_lor = NULL;
200  FiniteElementCollection *fec_lor = NULL;
201  FiniteElementSpace *fespace_lor = NULL;
202  if (pc_choice == LOR)
203  {
204  int basis_lor = basis;
205  if (basis == BasisType::Positive) { basis_lor=BasisType::ClosedUniform; }
206  mesh_lor = new Mesh(mesh, order, basis_lor);
207  fec_lor = new H1_FECollection(1, dim);
208  fespace_lor = new FiniteElementSpace(mesh_lor, fec_lor);
209  }
210
211  // 6. Check if the optimized version matches the given space
212  if (perf && !sol_fes_t::Matches(*fespace))
213  {
214  cout << "The given order does not match the optimized parameter.\n"
215  << "Recompile with suitable 'sol_p' value." << endl;
216  delete fespace;
217  delete fec;
218  delete mesh;
219  return 5;
220  }
221
222  // 7. Determine the list of true (i.e. conforming) essential boundary dofs.
223  // In this example, the boundary conditions are defined by marking all
224  // the boundary attributes from the mesh as essential (Dirichlet) and
225  // converting them to a list of true dofs.
226  Array<int> ess_tdof_list;
227  if (mesh->bdr_attributes.Size())
228  {
229  Array<int> ess_bdr(mesh->bdr_attributes.Max());
230  ess_bdr = 1;
231  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
232  }
233
234  // 8. Set up the linear form b(.) which corresponds to the right-hand side of
235  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
236  // the basis functions in the finite element fespace.
237  LinearForm *b = new LinearForm(fespace);
238  ConstantCoefficient one(1.0);
240  b->Assemble();
241
242  // 9. Define the solution vector x as a finite element grid function
243  // corresponding to fespace. Initialize x with initial guess of zero,
244  // which satisfies the boundary conditions.
245  GridFunction x(fespace);
246  x = 0.0;
247
248  // 10. Set up the bilinear form a(.,.) on the finite element space that will
249  // hold the matrix corresponding to the Laplacian operator -Delta.
250  // Optionally setup a form to be assembled for preconditioning (a_pc).
251  BilinearForm *a = new BilinearForm(fespace);
252  BilinearForm *a_pc = NULL;
253  if (pc_choice == LOR) { a_pc = new BilinearForm(fespace_lor); }
254  if (pc_choice == HO) { a_pc = new BilinearForm(fespace); }
255
256  // 11. Assemble the bilinear form and the corresponding linear system,
257  // applying any necessary transformations such as: eliminating boundary
258  // conditions, applying conforming constraints for non-conforming AMR,
259  // static condensation, etc.
260  if (static_cond)
261  {
263  MFEM_VERIFY(pc_choice != LOR,
264  "cannot use LOR preconditioner with static condensation");
265  }
266
267  cout << "Assembling the bilinear form ..." << flush;
268  tic_toc.Clear();
269  tic_toc.Start();
270  // Pre-allocate sparsity assuming dense element matrices
272
273  HPCBilinearForm *a_hpc = NULL;
274  Operator *a_oper = NULL;
275
276  if (!perf)
277  {
278  // Standard assembly using a diffusion domain integrator
280  a->Assemble();
281  }
282  else
283  {
284  // High-performance assembly/evaluation using the templated operator type
285  a_hpc = new HPCBilinearForm(integ_t(coeff_t(1.0)), *fespace);
286  if (matrix_free)
287  {
288  a_hpc->Assemble(); // partial assembly
289  }
290  else
291  {
292  a_hpc->AssembleBilinearForm(*a); // full matrix assembly
293  }
294  }
295  tic_toc.Stop();
296  cout << " done, " << tic_toc.RealTime() << "s." << endl;
297
298  // 12. Solve the system A X = B with CG. In the standard case, use a simple
299  // symmetric Gauss-Seidel preconditioner.
300
301  // Setup the operator matrix (if applicable)
302  SparseMatrix A;
303  Vector B, X;
304  if (perf && matrix_free)
305  {
306  a_hpc->FormLinearSystem(ess_tdof_list, x, *b, a_oper, X, B);
307  cout << "Size of linear system: " << a_hpc->Height() << endl;
308  }
309  else
310  {
311  a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
312  cout << "Size of linear system: " << A.Height() << endl;
313  a_oper = &A;
314  }
315
316  // Setup the matrix used for preconditioning
317  cout << "Assembling the preconditioning matrix ..." << flush;
318  tic_toc.Clear();
319  tic_toc.Start();
320
321  SparseMatrix A_pc;
322  if (pc_choice == LOR)
323  {
324  // TODO: assemble the LOR matrix using the performance code
326  a_pc->UsePrecomputedSparsity();
327  a_pc->Assemble();
328  a_pc->FormSystemMatrix(ess_tdof_list, A_pc);
329  }
330  else if (pc_choice == HO)
331  {
332  if (!matrix_free)
333  {
334  A_pc.MakeRef(A); // matrix already assembled, reuse it
335  }
336  else
337  {
338  a_pc->UsePrecomputedSparsity();
339  a_hpc->AssembleBilinearForm(*a_pc);
340  a_pc->FormSystemMatrix(ess_tdof_list, A_pc);
341  }
342  }
343
344  tic_toc.Stop();
345  cout << " done, " << tic_toc.RealTime() << "s." << endl;
346
347  // Solve with CG or PCG, depending if the matrix A_pc is available
348  if (pc_choice != NONE)
349  {
350  GSSmoother M(A_pc);
351  PCG(*a_oper, M, B, X, 1, 500, 1e-12, 0.0);
352  }
353  else
354  {
355  CG(*a_oper, B, X, 1, 500, 1e-12, 0.0);
356  }
357
358  // 13. Recover the solution as a finite element grid function.
359  if (perf && matrix_free)
360  {
361  a_hpc->RecoverFEMSolution(X, *b, x);
362  }
363  else
364  {
365  a->RecoverFEMSolution(X, *b, x);
366  }
367
368  // 14. Save the refined mesh and the solution. This output can be viewed later
369  // using GLVis: "glvis -m refined.mesh -g sol.gf".
370  ofstream mesh_ofs("refined.mesh");
371  mesh_ofs.precision(8);
372  mesh->Print(mesh_ofs);
373  ofstream sol_ofs("sol.gf");
374  sol_ofs.precision(8);
375  x.Save(sol_ofs);
376
377  // 15. Send the solution by socket to a GLVis server.
378  if (visualization)
379  {
380  char vishost[] = "localhost";
381  int visport = 19916;
382  socketstream sol_sock(vishost, visport);
383  sol_sock.precision(8);
384  sol_sock << "solution\n" << *mesh << x << flush;
385  }
386
387  // 16. Free the used memory.
388  delete a;
389  delete a_hpc;
390  if (a_oper != &A) { delete a_oper; }
391  delete a_pc;
392  delete b;
393  delete fespace;
394  delete fespace_lor;
395  delete fec_lor;
396  delete mesh_lor;
397  if (order > 0) { delete fec; }
398  delete mesh;
399
400  return 0;
401 }
int Size() const
Logical size of the array.
Definition: array.hpp:124
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:93
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1188
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
Subclass constant coefficient.
Definition: coefficient.hpp:67
void MakeRef(const SparseMatrix &master)
Clear the contents of the SparseMatrix and make it a reference to master.
Definition: sparsemat.cpp:226
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:79
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().
void AssembleBilinearForm(BilinearForm &a) const
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)
Definition: fespace.cpp:366
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:693
StopWatch tic_toc
Definition: tic_toc.cpp:447
double RealTime()
Definition: tic_toc.cpp:426
int main(int argc, char *argv[])
Definition: ex1.cpp:62
Data type for Gauss-Seidel smoother of sparse matrix.
H1_FiniteElement< geom, sol_p > sol_fe_t
Definition: ex1.cpp:52
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2731
TConstantCoefficient coeff_t
Definition: ex1.cpp:57
Data type sparse matrix.
Definition: sparsemat.hpp:40
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:57
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:7982
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:3924
int MeshGenerator()
Get the mesh generator/type.
Definition: mesh.hpp:686
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:516
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:529
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:384
int Dimension() const
Definition: mesh.hpp:744
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
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:189
H1_FiniteElementSpace< mesh_fe_t > mesh_fes_t
Definition: ex1.cpp:48
void UsePrecomputedSparsity(int ps=1)
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
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)
Definition: optparser.hpp:76
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
virtual const char * Name() const
Definition: fe_coll.hpp:53
const int mesh_p
Definition: ex1.cpp:41
int dim
Definition: ex24.cpp:43
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
Vector data type.
Definition: vector.hpp:48
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6203
TIntegrationRule< geom, ir_order > int_rule_t
Definition: ex1.cpp:56
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
Class for linear form - Vector with associated FE space and LFIntegrators.
Definition: linearform.hpp:23
TIntegrator< coeff_t, TDiffusionKernel > integ_t
Definition: ex1.cpp:58
Abstract operator.
Definition: operator.hpp:24
void EnableStaticCondensation()
const int ir_order
Definition: ex1.cpp:44
bool Good() const
Definition: optparser.hpp:125