MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex6.cpp
Go to the documentation of this file.
1 // MFEM Example 6
2 //
3 // Compile with: make ex6
4 //
5 // Sample runs: ex6 -m ../data/square-disc.mesh -o 1
6 // ex6 -m ../data/square-disc.mesh -o 2
7 // ex6 -m ../data/square-disc-nurbs.mesh -o 2
8 // ex6 -m ../data/star.mesh -o 3
9 // ex6 -m ../data/escher.mesh -o 2
10 // ex6 -m ../data/fichera.mesh -o 2
11 // ex6 -m ../data/disc-nurbs.mesh -o 2
12 // ex6 -m ../data/ball-nurbs.mesh
13 // ex6 -m ../data/pipe-nurbs.mesh
14 // ex6 -m ../data/star-surf.mesh -o 2
15 // ex6 -m ../data/square-disc-surf.mesh -o 2
16 // ex6 -m ../data/amr-quad.mesh
17 //
18 // Device sample runs:
19 // ex6 -pa -d cuda
20 // ex6 -pa -d occa-cuda
21 // ex6 -pa -d raja-omp
22 // ex6 -pa -d ceed-cpu
23 // * ex6 -pa -d ceed-cuda
24 // ex6 -pa -d ceed-cuda:/gpu/cuda/shared
25 //
26 // Description: This is a version of Example 1 with a simple adaptive mesh
27 // refinement loop. The problem being solved is again the Laplace
28 // equation -Delta u = 1 with homogeneous Dirichlet boundary
29 // conditions. The problem is solved on a sequence of meshes which
30 // are locally refined in a conforming (triangles, tetrahedrons)
31 // or non-conforming (quadrilaterals, hexahedra) manner according
32 // to a simple ZZ error estimator.
33 //
34 // The example demonstrates MFEM's capability to work with both
35 // conforming and nonconforming refinements, in 2D and 3D, on
36 // linear, curved and surface meshes. Interpolation of functions
37 // from coarse to fine meshes, as well as persistent GLVis
38 // visualization are also illustrated.
39 //
40 // We recommend viewing Example 1 before viewing this example.
41 
42 #include "mfem.hpp"
43 #include <fstream>
44 #include <iostream>
45 
46 using namespace std;
47 using namespace mfem;
48 
49 int main(int argc, char *argv[])
50 {
51  // 1. Parse command-line options.
52  const char *mesh_file = "../data/star.mesh";
53  int order = 1;
54  bool pa = false;
55  const char *device_config = "cpu";
56  bool visualization = true;
57 
58  OptionsParser args(argc, argv);
59  args.AddOption(&mesh_file, "-m", "--mesh",
60  "Mesh file to use.");
61  args.AddOption(&order, "-o", "--order",
62  "Finite element order (polynomial degree).");
63  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
64  "--no-partial-assembly", "Enable Partial Assembly.");
65  args.AddOption(&device_config, "-d", "--device",
66  "Device configuration string, see Device::Configure().");
67  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
68  "--no-visualization",
69  "Enable or disable GLVis visualization.");
70  args.Parse();
71  if (!args.Good())
72  {
73  args.PrintUsage(cout);
74  return 1;
75  }
76  args.PrintOptions(cout);
77 
78  // 2. Enable hardware devices such as GPUs, and programming models such as
79  // CUDA, OCCA, RAJA and OpenMP based on command line options.
80  Device device(device_config);
81  device.Print();
82 
83  // 3. Read the mesh from the given mesh file. We can handle triangular,
84  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
85  // the same code.
86  Mesh mesh(mesh_file, 1, 1);
87  int dim = mesh.Dimension();
88  int sdim = mesh.SpaceDimension();
89 
90  // 4. Since a NURBS mesh can currently only be refined uniformly, we need to
91  // convert it to a piecewise-polynomial curved mesh. First we refine the
92  // NURBS mesh a bit more and then project the curvature to quadratic Nodes.
93  if (mesh.NURBSext)
94  {
95  for (int i = 0; i < 2; i++)
96  {
97  mesh.UniformRefinement();
98  }
99  mesh.SetCurvature(2);
100  }
101 
102  // 5. Define a finite element space on the mesh. The polynomial order is
103  // one (linear) by default, but this can be changed on the command line.
104  H1_FECollection fec(order, dim);
105  FiniteElementSpace fespace(&mesh, &fec);
106 
107  // 6. As in Example 1, we set up bilinear and linear forms corresponding to
108  // the Laplace problem -\Delta u = 1. We don't assemble the discrete
109  // problem yet, this will be done in the main loop.
110  BilinearForm a(&fespace);
111  if (pa)
112  {
113  a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
114  a.SetDiagonalPolicy(Operator::DIAG_ONE);
115  }
116  LinearForm b(&fespace);
117 
118  ConstantCoefficient one(1.0);
119  ConstantCoefficient zero(0.0);
120 
122  a.AddDomainIntegrator(integ);
124 
125  // 7. The solution vector x and the associated finite element grid function
126  // will be maintained over the AMR iterations. We initialize it to zero.
127  GridFunction x(&fespace);
128  x = 0.0;
129 
130  // 8. All boundary attributes will be used for essential (Dirichlet) BC.
131  MFEM_VERIFY(mesh.bdr_attributes.Size() > 0,
132  "Boundary attributes required in the mesh.");
133  Array<int> ess_bdr(mesh.bdr_attributes.Max());
134  ess_bdr = 1;
135 
136  // 9. Connect to GLVis.
137  char vishost[] = "localhost";
138  int visport = 19916;
139  socketstream sol_sock;
140  if (visualization)
141  {
142  sol_sock.open(vishost, visport);
143  }
144 
145  // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
146  // that uses the ComputeElementFlux method of the DiffusionIntegrator to
147  // recover a smoothed flux (gradient) that is subtracted from the element
148  // flux to get an error indicator. We need to supply the space for the
149  // smoothed flux: an (H1)^sdim (i.e., vector-valued) space is used here.
150  FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
151  ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
152  estimator.SetAnisotropic();
153 
154  // 11. A refiner selects and refines elements based on a refinement strategy.
155  // The strategy here is to refine elements with errors larger than a
156  // fraction of the maximum element error. Other strategies are possible.
157  // The refiner will call the given error estimator.
158  ThresholdRefiner refiner(estimator);
159  refiner.SetTotalErrorFraction(0.7);
160 
161  // 12. The main AMR loop. In each iteration we solve the problem on the
162  // current mesh, visualize the solution, and refine the mesh.
163  const int max_dofs = 50000;
164  for (int it = 0; ; it++)
165  {
166  int cdofs = fespace.GetTrueVSize();
167  cout << "\nAMR iteration " << it << endl;
168  cout << "Number of unknowns: " << cdofs << endl;
169 
170  // 13. Assemble the right-hand side.
171  b.Assemble();
172 
173  // 14. Set Dirichlet boundary values in the GridFunction x.
174  // Determine the list of Dirichlet true DOFs in the linear system.
175  Array<int> ess_tdof_list;
176  x.ProjectBdrCoefficient(zero, ess_bdr);
177  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
178 
179  // 15. Assemble the stiffness matrix.
180  a.Assemble();
181 
182  // 16. Create the linear system: eliminate boundary conditions, constrain
183  // hanging nodes and possibly apply other transformations. The system
184  // will be solved for true (unconstrained) DOFs only.
185  OperatorPtr A;
186  Vector B, X;
187 
188  const int copy_interior = 1;
189  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
190 
191  // 17. Solve the linear system A X = B.
192  if (!pa)
193  {
194 #ifndef MFEM_USE_SUITESPARSE
195  // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
196  GSSmoother M((SparseMatrix&)(*A));
197  PCG(*A, M, B, X, 3, 200, 1e-12, 0.0);
198 #else
199  // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
200  UMFPackSolver umf_solver;
201  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
202  umf_solver.SetOperator(*A);
203  umf_solver.Mult(B, X);
204 #endif
205  }
206  else // Diagonal preconditioning in partial assembly mode.
207  {
208  OperatorJacobiSmoother M(a, ess_tdof_list);
209  PCG(*A, M, B, X, 3, 2000, 1e-12, 0.0);
210  }
211 
212  // 18. After solving the linear system, reconstruct the solution as a
213  // finite element GridFunction. Constrained nodes are interpolated
214  // from true DOFs (it may therefore happen that x.Size() >= X.Size()).
215  a.RecoverFEMSolution(X, b, x);
216 
217  // 19. Send solution by socket to the GLVis server.
218  if (visualization && sol_sock.good())
219  {
220  sol_sock.precision(8);
221  sol_sock << "solution\n" << mesh << x << flush;
222  }
223 
224  if (cdofs > max_dofs)
225  {
226  cout << "Reached the maximum number of dofs. Stop." << endl;
227  break;
228  }
229 
230  // 20. Call the refiner to modify the mesh. The refiner calls the error
231  // estimator to obtain element errors, then it selects elements to be
232  // refined and finally it modifies the mesh. The Stop() method can be
233  // used to determine if a stopping criterion was met.
234  refiner.Apply(mesh);
235  if (refiner.Stop())
236  {
237  cout << "Stopping criterion satisfied. Stop." << endl;
238  break;
239  }
240 
241  // 21. Update the space to reflect the new state of the mesh. Also,
242  // interpolate the solution x so that it lies in the new space but
243  // represents the same function. This saves solver iterations later
244  // since we'll have a good initial guess of x in the next step.
245  // Internally, FiniteElementSpace::Update() calculates an
246  // interpolation matrix which is then used by GridFunction::Update().
247  fespace.Update();
248  x.Update();
249 
250  // 22. Inform also the bilinear and linear forms that the space has
251  // changed.
252  a.Update();
253  b.Update();
254  }
255 
256  return 0;
257 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:124
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:93
void SetAnisotropic(bool aniso=true)
Enable/disable anisotropic estimates. To enable this option, the BilinearFormIntegrator must support ...
Definition: estimators.hpp:158
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:2233
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 Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:79
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
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:415
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:261
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
int main(int argc, char *argv[])
Definition: ex1.cpp:66
void SetDiagonalPolicy(DiagonalPolicy policy)
Sets diagonal policy used upon construction of the linear system.
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:731
bool Apply(Mesh &mesh)
Perform the mesh operation.
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:46
constexpr char vishost[]
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:128
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:8382
constexpr int visport
Mesh refinement operator using an error threshold.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:4165
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
virtual void Update(FiniteElementSpace *nfes=NULL)
Update the FiniteElementSpace and delete all data associated with the old one.
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:733
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:403
int Dimension() const
Definition: mesh.hpp:788
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:434
int SpaceDimension() const
Definition: mesh.hpp:789
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:31
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:201
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:742
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)
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 Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:152
void Update()
Update the object according to the associated FE space fes.
Definition: linearform.hpp:166
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:203
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
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:304
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
The ZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu error estimation procedure...
Definition: estimators.hpp:81
Vector data type.
Definition: vector.hpp:51
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:159
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:23
void SetTotalErrorFraction(double fraction)
Set the total fraction used in the computation of the threshold. The default value is 1/2...
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:118
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:2818
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
Definition: gridfunc.hpp:399
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:2721
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:145