MFEM  v3.2 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  {
98  }
99  for (int j = 0; j < Nelem; j++)
100  {
101  int attribute = j + 1;
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  {
122  }
123  for (int j = 0; j < Nelem; j++)
124  {
125  int attribute = j + 1;
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  for (int l = 0; l < 5; l++)
154  {
155  mesh->RefineAtVertex(Vertex(0, 0, 1));
156  }
157  SnapNodes(*mesh);
158  }
159  else if (amr == 2)
160  {
161  for (int l = 0; l < 4; l++)
162  {
163  mesh->RandomRefinement(0.5); // 50% probability
164  }
165  SnapNodes(*mesh);
166  }
167
168  // 4. Define a finite element space on the mesh. Here we use isoparametric
169  // finite elements -- the same as the mesh nodes.
170  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, &fec);
171  cout << "Number of unknowns: " << fespace->GetTrueVSize() << endl;
172
173  // 5. Set up the linear form b(.) which corresponds to the right-hand side of
174  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
175  // the basis functions in the finite element fespace.
176  LinearForm *b = new LinearForm(fespace);
177  ConstantCoefficient one(1.0);
181  b->Assemble();
182
183  // 6. Define the solution vector x as a finite element grid function
184  // corresponding to fespace. Initialize x with initial guess of zero.
185  GridFunction x(fespace);
186  x = 0.0;
187
188  // 7. Set up the bilinear form a(.,.) on the finite element space
189  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
190  // and Mass domain integrators.
191  BilinearForm *a = new BilinearForm(fespace);
194
195  // 8. Assemble the linear system, apply conforming constraints, etc.
196  a->Assemble();
197  SparseMatrix A;
198  Vector B, X;
199  Array<int> empty_tdof_list;
200  a->FormLinearSystem(empty_tdof_list, x, *b, A, X, B);
201
202 #ifndef MFEM_USE_SUITESPARSE
203  // 9. Define a simple symmetric Gauss-Seidel preconditioner and use it to
204  // solve the system AX=B with PCG.
205  GSSmoother M(A);
206  PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
207 #else
208  // 9. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
209  UMFPackSolver umf_solver;
210  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
211  umf_solver.SetOperator(A);
212  umf_solver.Mult(B, X);
213 #endif
214
215  // 10. Recover the solution as a finite element grid function.
216  a->RecoverFEMSolution(X, *b, x);
217
218  // 11. Compute and print the L^2 norm of the error.
219  cout<<"\nL2 norm of error: " << x.ComputeL2Error(sol_coef) << endl;
220
221  // 12. Save the refined mesh and the solution. This output can be viewed
222  // later using GLVis: "glvis -m sphere_refined.mesh -g sol.gf".
223  {
224  ofstream mesh_ofs("sphere_refined.mesh");
225  mesh_ofs.precision(8);
226  mesh->Print(mesh_ofs);
227  ofstream sol_ofs("sol.gf");
228  sol_ofs.precision(8);
229  x.Save(sol_ofs);
230  }
231
232  // 13. Send the solution by socket to a GLVis server.
233  if (visualization)
234  {
235  char vishost[] = "localhost";
236  int visport = 19916;
237  socketstream sol_sock(vishost, visport);
238  sol_sock.precision(8);
239  sol_sock << "solution\n" << *mesh << x << flush;
240  }
241
242  // 14. Free the used memory.
243  delete a;
244  delete b;
245  delete fespace;
246  delete mesh;
247
248  return 0;
249 }
250
251 double analytic_solution(const Vector &x)
252 {
253  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
254  return x(0)*x(1)/l2;
255 }
256
257 double analytic_rhs(const Vector &x)
258 {
259  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
260  return 7*x(0)*x(1)/l2;
261 }
262
263 void SnapNodes(Mesh &mesh)
264 {
265  GridFunction &nodes = *mesh.GetNodes();
266  Vector node(mesh.SpaceDimension());
267  for (int i = 0; i < nodes.FESpace()->GetNDofs(); i++)
268  {
269  for (int d = 0; d < mesh.SpaceDimension(); d++)
270  {
271  node(d) = nodes(nodes.FESpace()->DofToVDof(i, d));
272  }
273
274  node /= node.Norml2();
275
276  for (int d = 0; d < mesh.SpaceDimension(); d++)
277  {
278  nodes(nodes.FESpace()->DofToVDof(i, d)) = node(d);
279  }
280  }
281 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:159
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:103
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:644
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:34
void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Data type for vertex.
Definition: vertex.hpp:21
int main(int argc, char *argv[])
void FormLinearSystem(Array< int > &ess_tdof_list, Vector &x, Vector &b, SparseMatrix &A, Vector &X, Vector &B, int copy_interior=0)
Data type for Gauss-Seidel smoother of sparse matrix.
void AddVertex(const double *)
Definition: mesh.cpp:842
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:333
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2116
void FinalizeTriMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Definition: mesh.cpp:989
double ComputeL2Error(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:193
Data type sparse matrix.
Definition: sparsemat.hpp:38
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6402
double analytic_solution(const Vector &x)
Definition: ex7.cpp:251
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:6085
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:3003
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:263
int Dimension() const
Definition: mesh.hpp:523
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
double analytic_rhs(const Vector &x)
Definition: ex7.cpp:257
int SpaceDimension() const
Definition: mesh.hpp:524
virtual void Print(std::ostream &out=std::cout) const
Print the mesh to the given stream using the default MFEM mesh format.
Definition: mesh.cpp:6736
Adds new Domain Integrator.
Definition: linearform.cpp:19
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:344
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
Definition: mesh.cpp:863
void RefineAtVertex(const Vertex &vert, double eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex. Uses GeneralRefinement.
Definition: mesh.cpp:6104
virtual int GetTrueVSize()
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:164
void AddTriangle(const int *vi, int attr=1)
Definition: mesh.cpp:858
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Definition: mesh.cpp:1020
Adds new Domain Integrator.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
void SnapNodes(Mesh &mesh)
Definition: ex7.cpp:263
class for C-function coefficient
Vector data type.
Definition: vector.hpp:33
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5114
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:77
Class for linear form - Vector with associated FE space and LFIntegrators.
Definition: linearform.hpp:23
virtual void Mult(const Vector &b, Vector &x) const
Operator application.
Definition: solvers.cpp:1741
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:1641
bool Good() const
Definition: optparser.hpp:120