MFEM  v4.5.2
Finite element discretization library
ex24.cpp
Go to the documentation of this file.
1 // MFEM Example 24
2 //
3 // Compile with: make ex24
4 //
5 // Sample runs: ex24 -m ../data/star.mesh
6 // ex24 -m ../data/square-disc.mesh -o 2
7 // ex24 -m ../data/beam-tet.mesh
8 // ex24 -m ../data/beam-hex.mesh -o 2 -pa
9 // ex24 -m ../data/beam-hex.mesh -o 2 -pa -p 1
10 // ex24 -m ../data/beam-hex.mesh -o 2 -pa -p 2
11 // ex24 -m ../data/escher.mesh
12 // ex24 -m ../data/escher.mesh -o 2
13 // ex24 -m ../data/fichera.mesh
14 // ex24 -m ../data/fichera-q2.vtk
15 // ex24 -m ../data/fichera-q3.mesh
16 // ex24 -m ../data/square-disc-nurbs.mesh
17 // ex24 -m ../data/beam-hex-nurbs.mesh
18 // ex24 -m ../data/amr-quad.mesh -o 2
19 // ex24 -m ../data/amr-hex.mesh
20 //
21 // Device sample runs:
22 // ex24 -m ../data/star.mesh -pa -d cuda
23 // ex24 -m ../data/star.mesh -pa -d raja-cuda
24 // ex24 -m ../data/star.mesh -pa -d raja-omp
25 // ex24 -m ../data/beam-hex.mesh -pa -d cuda
26 //
27 // Description: This example code illustrates usage of mixed finite element
28 // spaces, with three variants:
29 //
30 // 1) (grad p, u) for p in H^1 tested against u in H(curl)
31 // 2) (curl v, u) for v in H(curl) tested against u in H(div), 3D
32 // 3) (div v, q) for v in H(div) tested against q in L_2
33 //
34 // Using different approaches, we project the gradient, curl, or
35 // divergence to the appropriate space.
36 //
37 // We recommend viewing examples 1, 3, and 5 before viewing this
38 // example.
39 
40 #include "mfem.hpp"
41 #include <fstream>
42 #include <iostream>
43 
44 using namespace std;
45 using namespace mfem;
46 
47 double p_exact(const Vector &x);
48 void gradp_exact(const Vector &, Vector &);
49 double div_gradp_exact(const Vector &x);
50 void v_exact(const Vector &x, Vector &v);
51 void curlv_exact(const Vector &x, Vector &cv);
52 
53 int dim;
54 double freq = 1.0, kappa;
55 
56 int main(int argc, char *argv[])
57 {
58  // 1. Parse command-line options.
59  const char *mesh_file = "../data/beam-hex.mesh";
60  int order = 1;
61  int prob = 0;
62  bool static_cond = false;
63  bool pa = false;
64  const char *device_config = "cpu";
65  bool visualization = 1;
66 
67  OptionsParser args(argc, argv);
68  args.AddOption(&mesh_file, "-m", "--mesh",
69  "Mesh file to use.");
70  args.AddOption(&order, "-o", "--order",
71  "Finite element order (polynomial degree).");
72  args.AddOption(&prob, "-p", "--problem-type",
73  "Choose between 0: grad, 1: curl, 2: div");
74  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
75  "--no-static-condensation", "Enable static condensation.");
76  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
77  "--no-partial-assembly", "Enable Partial Assembly.");
78  args.AddOption(&device_config, "-d", "--device",
79  "Device configuration string, see Device::Configure().");
80  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
81  "--no-visualization",
82  "Enable or disable GLVis visualization.");
83 
84  args.Parse();
85  if (!args.Good())
86  {
87  args.PrintUsage(cout);
88  return 1;
89  }
90  args.PrintOptions(cout);
91  kappa = freq * M_PI;
92 
93  // 2. Enable hardware devices such as GPUs, and programming models such as
94  // CUDA, OCCA, RAJA and OpenMP based on command line options.
95  Device device(device_config);
96  device.Print();
97 
98  // 3. Read the mesh from the given mesh file. We can handle triangular,
99  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
100  // the same code.
101  Mesh *mesh = new Mesh(mesh_file, 1, 1);
102  dim = mesh->Dimension();
103  int sdim = mesh->SpaceDimension();
104 
105  // 4. Refine the mesh to increase the resolution. In this example we do
106  // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
107  // largest number that gives a final mesh with no more than 50,000
108  // elements.
109  {
110  int ref_levels = (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
111  for (int l = 0; l < ref_levels; l++)
112  {
113  mesh->UniformRefinement();
114  }
115  }
116 
117  // 5. Define a finite element space on the mesh. Here we use Nedelec or
118  // Raviart-Thomas finite elements of the specified order.
119  FiniteElementCollection *trial_fec = NULL;
120  FiniteElementCollection *test_fec = NULL;
121 
122  if (prob == 0)
123  {
124  trial_fec = new H1_FECollection(order, dim);
125  test_fec = new ND_FECollection(order, dim);
126  }
127  else if (prob == 1)
128  {
129  trial_fec = new ND_FECollection(order, dim);
130  test_fec = new RT_FECollection(order-1, dim);
131  }
132  else
133  {
134  trial_fec = new RT_FECollection(order-1, dim);
135  test_fec = new L2_FECollection(order-1, dim);
136  }
137 
138  FiniteElementSpace trial_fes(mesh, trial_fec);
139  FiniteElementSpace test_fes(mesh, test_fec);
140 
141  int trial_size = trial_fes.GetTrueVSize();
142  int test_size = test_fes.GetTrueVSize();
143 
144  if (prob == 0)
145  {
146  cout << "Number of Nedelec finite element unknowns: " << test_size << endl;
147  cout << "Number of H1 finite element unknowns: " << trial_size << endl;
148  }
149  else if (prob == 1)
150  {
151  cout << "Number of Nedelec finite element unknowns: " << trial_size << endl;
152  cout << "Number of Raviart-Thomas finite element unknowns: " << test_size <<
153  endl;
154  }
155  else
156  {
157  cout << "Number of Raviart-Thomas finite element unknowns: "
158  << trial_size << endl;
159  cout << "Number of L2 finite element unknowns: " << test_size << endl;
160  }
161 
162  // 6. Define the solution vector as a finite element grid function
163  // corresponding to the trial fespace.
164  GridFunction gftest(&test_fes);
165  GridFunction gftrial(&trial_fes);
166  GridFunction x(&test_fes);
168  VectorFunctionCoefficient gradp_coef(sdim, gradp_exact);
169  VectorFunctionCoefficient v_coef(sdim, v_exact);
170  VectorFunctionCoefficient curlv_coef(sdim, curlv_exact);
171  FunctionCoefficient divgradp_coef(div_gradp_exact);
172 
173  if (prob == 0)
174  {
175  gftrial.ProjectCoefficient(p_coef);
176  }
177  else if (prob == 1)
178  {
179  gftrial.ProjectCoefficient(v_coef);
180  }
181  else
182  {
183  gftrial.ProjectCoefficient(gradp_coef);
184  }
185 
186  gftrial.SetTrueVector();
187  gftrial.SetFromTrueVector();
188 
189  // 7. Set up the bilinear forms for L2 projection.
190  ConstantCoefficient one(1.0);
191  BilinearForm a(&test_fes);
192  MixedBilinearForm a_mixed(&trial_fes, &test_fes);
193  if (pa)
194  {
195  a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
196  a_mixed.SetAssemblyLevel(AssemblyLevel::PARTIAL);
197  }
198 
199  if (prob == 0)
200  {
201  a.AddDomainIntegrator(new VectorFEMassIntegrator(one));
203  }
204  else if (prob == 1)
205  {
206  a.AddDomainIntegrator(new VectorFEMassIntegrator(one));
208  }
209  else
210  {
211  a.AddDomainIntegrator(new MassIntegrator(one));
213  }
214 
215  // 8. Assemble the bilinear form and the corresponding linear system,
216  // applying any necessary transformations such as: eliminating boundary
217  // conditions, applying conforming constraints for non-conforming AMR,
218  // static condensation, etc.
219  if (static_cond) { a.EnableStaticCondensation(); }
220 
221  a.Assemble();
222  if (!pa) { a.Finalize(); }
223 
224  a_mixed.Assemble();
225  if (!pa) { a_mixed.Finalize(); }
226 
227  if (pa)
228  {
229  a_mixed.Mult(gftrial, x);
230  }
231  else
232  {
233  SparseMatrix& mixed = a_mixed.SpMat();
234  mixed.Mult(gftrial, x);
235  }
236 
237  // 9. Define and apply a PCG solver for Ax = b with Jacobi preconditioner.
238  {
239  GridFunction rhs(&test_fes);
240  rhs = x;
241  x = 0.0;
242 
243  CGSolver cg;
244  cg.SetRelTol(1e-12);
245  cg.SetMaxIter(1000);
246  cg.SetPrintLevel(1);
247  if (pa)
248  {
249  Array<int> ess_tdof_list; // empty
250  OperatorJacobiSmoother Jacobi(a, ess_tdof_list);
251 
252  cg.SetOperator(a);
253  cg.SetPreconditioner(Jacobi);
254  cg.Mult(rhs, x);
255  }
256  else
257  {
258  SparseMatrix& Amat = a.SpMat();
259  DSmoother Jacobi(Amat);
260 
261  cg.SetOperator(Amat);
262  cg.SetPreconditioner(Jacobi);
263  cg.Mult(rhs, x);
264  }
265  }
266 
267  // 10. Compute the same field by applying a DiscreteInterpolator.
268  GridFunction discreteInterpolant(&test_fes);
269  DiscreteLinearOperator dlo(&trial_fes, &test_fes);
270  if (prob == 0)
271  {
273  }
274  else if (prob == 1)
275  {
277  }
278  else
279  {
281  }
282 
283  dlo.Assemble();
284  dlo.Mult(gftrial, discreteInterpolant);
285 
286  // 11. Compute the projection of the exact field.
287  GridFunction exact_proj(&test_fes);
288  if (prob == 0)
289  {
290  exact_proj.ProjectCoefficient(gradp_coef);
291  }
292  else if (prob == 1)
293  {
294  exact_proj.ProjectCoefficient(curlv_coef);
295  }
296  else
297  {
298  exact_proj.ProjectCoefficient(divgradp_coef);
299  }
300 
301  exact_proj.SetTrueVector();
302  exact_proj.SetFromTrueVector();
303 
304  // 12. Compute and print the L_2 norm of the error.
305  if (prob == 0)
306  {
307  double errSol = x.ComputeL2Error(gradp_coef);
308  double errInterp = discreteInterpolant.ComputeL2Error(gradp_coef);
309  double errProj = exact_proj.ComputeL2Error(gradp_coef);
310 
311  cout << "\n Solution of (E_h,v) = (grad p_h,v) for E_h and v in H(curl): "
312  "|| E_h - grad p ||_{L_2} = " << errSol << '\n' << endl;
313  cout << " Gradient interpolant E_h = grad p_h in H(curl): || E_h - grad p"
314  " ||_{L_2} = " << errInterp << '\n' << endl;
315  cout << " Projection E_h of exact grad p in H(curl): || E_h - grad p "
316  "||_{L_2} = " << errProj << '\n' << endl;
317  }
318  else if (prob == 1)
319  {
320  double errSol = x.ComputeL2Error(curlv_coef);
321  double errInterp = discreteInterpolant.ComputeL2Error(curlv_coef);
322  double errProj = exact_proj.ComputeL2Error(curlv_coef);
323 
324  cout << "\n Solution of (E_h,w) = (curl v_h,w) for E_h and w in H(div): "
325  "|| E_h - curl v ||_{L_2} = " << errSol << '\n' << endl;
326  cout << " Curl interpolant E_h = curl v_h in H(div): || E_h - curl v "
327  "||_{L_2} = " << errInterp << '\n' << endl;
328  cout << " Projection E_h of exact curl v in H(div): || E_h - curl v "
329  "||_{L_2} = " << errProj << '\n' << endl;
330  }
331  else
332  {
333  int order_quad = max(2, 2*order+1);
334  const IntegrationRule *irs[Geometry::NumGeom];
335  for (int i=0; i < Geometry::NumGeom; ++i)
336  {
337  irs[i] = &(IntRules.Get(i, order_quad));
338  }
339 
340  double errSol = x.ComputeL2Error(divgradp_coef, irs);
341  double errInterp = discreteInterpolant.ComputeL2Error(divgradp_coef, irs);
342  double errProj = exact_proj.ComputeL2Error(divgradp_coef, irs);
343 
344  cout << "\n Solution of (f_h,q) = (div v_h,q) for f_h and q in L_2: "
345  "|| f_h - div v ||_{L_2} = " << errSol << '\n' << endl;
346  cout << " Divergence interpolant f_h = div v_h in L_2: || f_h - div v "
347  "||_{L_2} = " << errInterp << '\n' << endl;
348  cout << " Projection f_h of exact div v in L_2: || f_h - div v "
349  "||_{L_2} = " << errProj << '\n' << endl;
350  }
351 
352  // 13. Save the refined mesh and the solution. This output can be viewed
353  // later using GLVis: "glvis -m refined.mesh -g sol.gf".
354  ofstream mesh_ofs("refined.mesh");
355  mesh_ofs.precision(8);
356  mesh->Print(mesh_ofs);
357  ofstream sol_ofs("sol.gf");
358  sol_ofs.precision(8);
359  x.Save(sol_ofs);
360 
361  // 14. Send the solution by socket to a GLVis server.
362  if (visualization)
363  {
364  char vishost[] = "localhost";
365  int visport = 19916;
366  socketstream sol_sock(vishost, visport);
367  sol_sock.precision(8);
368  sol_sock << "solution\n" << *mesh << x << flush;
369  }
370 
371  // 15. Free the used memory.
372  delete trial_fec;
373  delete test_fec;
374  delete mesh;
375 
376  return 0;
377 }
378 
379 double p_exact(const Vector &x)
380 {
381  if (dim == 3)
382  {
383  return sin(x(0)) * sin(x(1)) * sin(x(2));
384  }
385  else if (dim == 2)
386  {
387  return sin(x(0)) * sin(x(1));
388  }
389 
390  return 0.0;
391 }
392 
393 void gradp_exact(const Vector &x, Vector &f)
394 {
395  if (dim == 3)
396  {
397  f(0) = cos(x(0)) * sin(x(1)) * sin(x(2));
398  f(1) = sin(x(0)) * cos(x(1)) * sin(x(2));
399  f(2) = sin(x(0)) * sin(x(1)) * cos(x(2));
400  }
401  else
402  {
403  f(0) = cos(x(0)) * sin(x(1));
404  f(1) = sin(x(0)) * cos(x(1));
405  if (x.Size() == 3) { f(2) = 0.0; }
406  }
407 }
408 
409 double div_gradp_exact(const Vector &x)
410 {
411  if (dim == 3)
412  {
413  return -3.0 * sin(x(0)) * sin(x(1)) * sin(x(2));
414  }
415  else if (dim == 2)
416  {
417  return -2.0 * sin(x(0)) * sin(x(1));
418  }
419 
420  return 0.0;
421 }
422 
423 void v_exact(const Vector &x, Vector &v)
424 {
425  if (dim == 3)
426  {
427  v(0) = sin(kappa * x(1));
428  v(1) = sin(kappa * x(2));
429  v(2) = sin(kappa * x(0));
430  }
431  else
432  {
433  v(0) = sin(kappa * x(1));
434  v(1) = sin(kappa * x(0));
435  if (x.Size() == 3) { v(2) = 0.0; }
436  }
437 }
438 
439 void curlv_exact(const Vector &x, Vector &cv)
440 {
441  if (dim == 3)
442  {
443  cv(0) = -kappa * cos(kappa * x(2));
444  cv(1) = -kappa * cos(kappa * x(0));
445  cv(2) = -kappa * cos(kappa * x(1));
446  }
447  else
448  {
449  cv = 0.0;
450  }
451 }
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: gridfunc.cpp:2763
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:733
Conjugate gradient method.
Definition: solvers.hpp:493
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
Data type for scaled Jacobi-type smoother of sparse matrix.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:150
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:324
int Dimension() const
Definition: mesh.hpp:1047
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:712
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
double kappa
Definition: ex24.cpp:54
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
STL namespace.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
virtual void Mult(const Vector &x, Vector &y) const
Matrix multiplication: .
void AddDomainInterpolator(DiscreteInterpolator *di)
Adds a domain interpolator. Assumes ownership of di.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds a domain integrator. Assumes ownership of bfi.
int main(int argc, char *argv[])
Definition: ex24.cpp:56
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:144
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
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:50
void curlv_exact(const Vector &x, Vector &cv)
Definition: ex24.cpp:439
constexpr char vishost[]
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:302
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9878
constexpr int visport
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
void Assemble(int skip_zeros=1)
void v_exact(const Vector &x, Vector &v)
Definition: ex24.cpp:423
prob_type prob
Definition: ex25.cpp:153
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:373
double p_exact(const Vector &x)
Definition: ex24.cpp:379
A general vector function coefficient.
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::LEGACY.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:590
void SetRelTol(double rtol)
Definition: solvers.hpp:199
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
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
int SpaceDimension() const
Definition: mesh.hpp:1048
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:936
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
double freq
Definition: ex24.cpp:54
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2396
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
A general function coefficient.
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:447
Vector data type.
Definition: vector.hpp:60
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1749
double div_gradp_exact(const Vector &x)
Definition: ex24.cpp:409
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3673
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:320
double f(const Vector &p)
void gradp_exact(const Vector &, Vector &)
Definition: ex24.cpp:393