MFEM  v3.3
Finite element discretization library
ex7.cpp
Go to the documentation of this file.
1 // MFEM Example 7
2 //
3 // Compile with: make ex7
4 //
5 // Sample runs: ex7 -e 0 -o 2 -r 4
6 // ex7 -e 1 -o 2 -r 4 -snap
7 // ex7 -e 0 -amr 1
8 // ex7 -e 1 -amr 2 -o 2
9 //
10 // Description: This example code demonstrates the use of MFEM to define a
11 // triangulation of a unit sphere and a simple isoparametric
12 // finite element discretization of the Laplace problem with mass
13 // term, -Delta u + u = f.
14 //
15 // The example highlights mesh generation, the use of mesh
16 // refinement, high-order meshes and finite elements, as well as
17 // surface-based linear and bilinear forms corresponding to the
18 // left-hand side and right-hand side of the discrete linear
19 // system. Simple local mesh refinement is also demonstrated.
20 //
21 // We recommend viewing Example 1 before viewing this example.
22 
23 #include "mfem.hpp"
24 #include <fstream>
25 #include <iostream>
26 
27 using namespace std;
28 using namespace mfem;
29 
30 // Exact solution and r.h.s., see below for implementation.
31 double analytic_solution(const Vector &x);
32 double analytic_rhs(const Vector &x);
33 void SnapNodes(Mesh &mesh);
34 
35 int main(int argc, char *argv[])
36 {
37  // 1. Parse command-line options.
38  int elem_type = 1;
39  int ref_levels = 2;
40  int amr = 0;
41  int order = 2;
42  bool always_snap = false;
43  bool visualization = 1;
44 
45  OptionsParser args(argc, argv);
46  args.AddOption(&elem_type, "-e", "--elem",
47  "Type of elements to use: 0 - triangles, 1 - quads.");
48  args.AddOption(&order, "-o", "--order",
49  "Finite element order (polynomial degree).");
50  args.AddOption(&ref_levels, "-r", "--refine",
51  "Number of times to refine the mesh uniformly.");
52  args.AddOption(&amr, "-amr", "--refine-locally",
53  "Additional local (non-conforming) refinement:"
54  " 1 = refine around north pole, 2 = refine randomly.");
55  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
56  "--no-visualization",
57  "Enable or disable GLVis visualization.");
58  args.AddOption(&always_snap, "-snap", "--always-snap", "-no-snap",
59  "--snap-at-the-end",
60  "If true, snap nodes to the sphere initially and after each refinement "
61  "otherwise, snap only after the last refinement");
62  args.Parse();
63  if (!args.Good())
64  {
65  args.PrintUsage(cout);
66  return 1;
67  }
68  args.PrintOptions(cout);
69 
70  // 2. Generate an initial high-order (surface) mesh on the unit sphere. The
71  // Mesh object represents a 2D mesh in 3 spatial dimensions. We first add
72  // the elements and the vertices of the mesh, and then make it high-order
73  // by specifying a finite element space for its nodes.
74  int Nvert = 8, Nelem = 6;
75  if (elem_type == 0)
76  {
77  Nvert = 6;
78  Nelem = 8;
79  }
80  Mesh *mesh = new Mesh(2, Nvert, Nelem, 0, 3);
81 
82  if (elem_type == 0) // inscribed octahedron
83  {
84  const double tri_v[6][3] =
85  {
86  { 1, 0, 0}, { 0, 1, 0}, {-1, 0, 0},
87  { 0, -1, 0}, { 0, 0, 1}, { 0, 0, -1}
88  };
89  const int tri_e[8][3] =
90  {
91  {0, 1, 4}, {1, 2, 4}, {2, 3, 4}, {3, 0, 4},
92  {1, 0, 5}, {2, 1, 5}, {3, 2, 5}, {0, 3, 5}
93  };
94 
95  for (int j = 0; j < Nvert; j++)
96  {
97  mesh->AddVertex(tri_v[j]);
98  }
99  for (int j = 0; j < Nelem; j++)
100  {
101  int attribute = j + 1;
102  mesh->AddTriangle(tri_e[j], attribute);
103  }
104  mesh->FinalizeTriMesh(1, 1, true);
105  }
106  else // inscribed cube
107  {
108  const double quad_v[8][3] =
109  {
110  {-1, -1, -1}, {+1, -1, -1}, {+1, +1, -1}, {-1, +1, -1},
111  {-1, -1, +1}, {+1, -1, +1}, {+1, +1, +1}, {-1, +1, +1}
112  };
113  const int quad_e[6][4] =
114  {
115  {3, 2, 1, 0}, {0, 1, 5, 4}, {1, 2, 6, 5},
116  {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}
117  };
118 
119  for (int j = 0; j < Nvert; j++)
120  {
121  mesh->AddVertex(quad_v[j]);
122  }
123  for (int j = 0; j < Nelem; j++)
124  {
125  int attribute = j + 1;
126  mesh->AddQuad(quad_e[j], attribute);
127  }
128  mesh->FinalizeQuadMesh(1, 1, true);
129  }
130 
131  // Set the space for the high-order mesh nodes.
132  H1_FECollection fec(order, mesh->Dimension());
133  FiniteElementSpace nodal_fes(mesh, &fec, mesh->SpaceDimension());
134  mesh->SetNodalFESpace(&nodal_fes);
135 
136  // 3. Refine the mesh while snapping nodes to the sphere.
137  for (int l = 0; l <= ref_levels; l++)
138  {
139  if (l > 0) // for l == 0 just perform snapping
140  {
141  mesh->UniformRefinement();
142  }
143 
144  // Snap the nodes of the refined mesh back to sphere surface.
145  if (always_snap || l == ref_levels)
146  {
147  SnapNodes(*mesh);
148  }
149  }
150 
151  if (amr == 1)
152  {
153  Vertex target(0.0, 0.0, 1.0);
154  for (int l = 0; l < 5; l++)
155  {
156  mesh->RefineAtVertex(target);
157  }
158  SnapNodes(*mesh);
159  }
160  else if (amr == 2)
161  {
162  for (int l = 0; l < 4; l++)
163  {
164  mesh->RandomRefinement(0.5); // 50% probability
165  }
166  SnapNodes(*mesh);
167  }
168 
169  // 4. Define a finite element space on the mesh. Here we use isoparametric
170  // finite elements -- the same as the mesh nodes.
171  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, &fec);
172  cout << "Number of unknowns: " << fespace->GetTrueVSize() << endl;
173 
174  // 5. Set up the linear form b(.) which corresponds to the right-hand side of
175  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
176  // the basis functions in the finite element fespace.
177  LinearForm *b = new LinearForm(fespace);
178  ConstantCoefficient one(1.0);
181  b->AddDomainIntegrator(new DomainLFIntegrator(rhs_coef));
182  b->Assemble();
183 
184  // 6. Define the solution vector x as a finite element grid function
185  // corresponding to fespace. Initialize x with initial guess of zero.
186  GridFunction x(fespace);
187  x = 0.0;
188 
189  // 7. Set up the bilinear form a(.,.) on the finite element space
190  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
191  // and Mass domain integrators.
192  BilinearForm *a = new BilinearForm(fespace);
194  a->AddDomainIntegrator(new MassIntegrator(one));
195 
196  // 8. Assemble the linear system, apply conforming constraints, etc.
197  a->Assemble();
198  SparseMatrix A;
199  Vector B, X;
200  Array<int> empty_tdof_list;
201  a->FormLinearSystem(empty_tdof_list, x, *b, A, X, B);
202 
203 #ifndef MFEM_USE_SUITESPARSE
204  // 9. Define a simple symmetric Gauss-Seidel preconditioner and use it to
205  // solve the system AX=B with PCG.
206  GSSmoother M(A);
207  PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
208 #else
209  // 9. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
210  UMFPackSolver umf_solver;
211  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
212  umf_solver.SetOperator(A);
213  umf_solver.Mult(B, X);
214 #endif
215 
216  // 10. Recover the solution as a finite element grid function.
217  a->RecoverFEMSolution(X, *b, x);
218 
219  // 11. Compute and print the L^2 norm of the error.
220  cout<<"\nL2 norm of error: " << x.ComputeL2Error(sol_coef) << endl;
221 
222  // 12. Save the refined mesh and the solution. This output can be viewed
223  // later using GLVis: "glvis -m sphere_refined.mesh -g sol.gf".
224  {
225  ofstream mesh_ofs("sphere_refined.mesh");
226  mesh_ofs.precision(8);
227  mesh->Print(mesh_ofs);
228  ofstream sol_ofs("sol.gf");
229  sol_ofs.precision(8);
230  x.Save(sol_ofs);
231  }
232 
233  // 13. Send the solution by socket to a GLVis server.
234  if (visualization)
235  {
236  char vishost[] = "localhost";
237  int visport = 19916;
238  socketstream sol_sock(vishost, visport);
239  sol_sock.precision(8);
240  sol_sock << "solution\n" << *mesh << x << flush;
241  }
242 
243  // 14. Free the used memory.
244  delete a;
245  delete b;
246  delete fespace;
247  delete mesh;
248 
249  return 0;
250 }
251 
252 double analytic_solution(const Vector &x)
253 {
254  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
255  return x(0)*x(1)/l2;
256 }
257 
258 double analytic_rhs(const Vector &x)
259 {
260  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
261  return 7*x(0)*x(1)/l2;
262 }
263 
264 void SnapNodes(Mesh &mesh)
265 {
266  GridFunction &nodes = *mesh.GetNodes();
267  Vector node(mesh.SpaceDimension());
268  for (int i = 0; i < nodes.FESpace()->GetNDofs(); i++)
269  {
270  for (int d = 0; d < mesh.SpaceDimension(); d++)
271  {
272  node(d) = nodes(nodes.FESpace()->DofToVDof(i, d));
273  }
274 
275  node /= node.Norml2();
276 
277  for (int d = 0; d < mesh.SpaceDimension(); d++)
278  {
279  nodes(nodes.FESpace()->DofToVDof(i, d)) = node(d);
280  }
281  }
282 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:161
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:104
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
Subclass constant coefficient.
Definition: coefficient.hpp:57
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:667
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:42
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Data type for vertex.
Definition: vertex.hpp:21
STL namespace.
Data type for Gauss-Seidel smoother of sparse matrix.
void AddVertex(const double *)
Definition: mesh.cpp:910
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:337
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2185
void FinalizeTriMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a triangular Mesh.
Definition: mesh.cpp:1054
double ComputeL2Error(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:217
Data type sparse matrix.
Definition: sparsemat.hpp:38
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6684
double analytic_solution(const Vector &x)
Definition: ex7.cpp:252
void RandomRefinement(double prob, bool aniso=false, int nonconforming=-1, int nc_limit=0)
Refine each element with given probability. Uses GeneralRefinement.
Definition: mesh.cpp:6367
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:456
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3261
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:293
int Dimension() const
Definition: mesh.hpp:611
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
double analytic_rhs(const Vector &x)
Definition: ex7.cpp:258
int SpaceDimension() const
Definition: mesh.hpp:612
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator.
Definition: linearform.cpp:19
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:348
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:74
void AddQuad(const int *vi, int attr=1)
Definition: mesh.cpp:931
void RefineAtVertex(const Vertex &vert, double eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex. Uses GeneralRefinement.
Definition: mesh.cpp:6386
virtual int GetTrueVSize()
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:166
void AddTriangle(const int *vi, int attr=1)
Definition: mesh.cpp:926
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition: mesh.cpp:1086
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
void SnapNodes(Mesh &mesh)
Definition: ex7.cpp:264
class for C-function coefficient
int main(int argc, char *argv[])
Definition: ex7.cpp:35
Vector data type.
Definition: vector.hpp:36
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5368
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:146
Class for linear form - Vector with associated FE space and LFIntegrators.
Definition: linearform.hpp:23
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, SparseMatrix &A, Vector &X, Vector &B, int copy_interior=0)
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1774
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:1674
virtual void Print(std::ostream &out=std::cout) const
Definition: mesh.hpp:988
bool Good() const
Definition: optparser.hpp:120