MFEM  v3.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
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  mesh->RefineAtVertex(Vertex(0, 0, 1), 5); // 5 levels
154  SnapNodes(*mesh);
155  }
156  else if (amr == 2)
157  {
158  mesh->RandomRefinement(4, 2); // 4 levels, 1/2 of the elements
159  SnapNodes(*mesh);
160  }
161 
162  // 4. Define a finite element space on the mesh. Here we use isoparametric
163  // finite elements -- the same as the mesh nodes.
164  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, &fec);
165  cout << "Number of unknowns: " << fespace->GetTrueVSize() << endl;
166 
167  // 5. Set up the linear form b(.) which corresponds to the right-hand side of
168  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
169  // the basis functions in the finite element fespace.
170  LinearForm *b = new LinearForm(fespace);
171  ConstantCoefficient one(1.0);
174  b->AddDomainIntegrator(new DomainLFIntegrator(rhs_coef));
175  b->Assemble();
176 
177  // 6. Define the solution vector x as a finite element grid function
178  // corresponding to fespace. Initialize x with initial guess of zero.
179  GridFunction x(fespace);
180  x = 0.0;
181 
182  // 7. Set up the bilinear form a(.,.) on the finite element space
183  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
184  // and Mass domain integrators.
185  BilinearForm *a = new BilinearForm(fespace);
187  a->AddDomainIntegrator(new MassIntegrator(one));
188 
189  // 8. Assemble the linear system, apply conforming constraints, etc.
190  a->Assemble();
191  SparseMatrix A;
192  Vector B, X;
193  Array<int> empty_tdof_list;
194  a->FormLinearSystem(empty_tdof_list, x, *b, A, X, B);
195 
196 #ifndef MFEM_USE_SUITESPARSE
197  // 9. Define a simple symmetric Gauss-Seidel preconditioner and use it to
198  // solve the system AX=B with PCG.
199  GSSmoother M(A);
200  PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
201 #else
202  // 9. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
203  UMFPackSolver umf_solver;
204  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
205  umf_solver.SetOperator(A);
206  umf_solver.Mult(B, X);
207 #endif
208 
209  // 10. Recover the solution as a finite element grid function.
210  a->RecoverFEMSolution(X, *b, x);
211 
212  // 11. Compute and print the L^2 norm of the error.
213  cout<<"\nL2 norm of error: " << x.ComputeL2Error(sol_coef) << endl;
214 
215  // 12. Save the refined mesh and the solution. This output can be viewed
216  // later using GLVis: "glvis -m sphere_refined.mesh -g sol.gf".
217  {
218  ofstream mesh_ofs("sphere_refined.mesh");
219  mesh_ofs.precision(8);
220  mesh->Print(mesh_ofs);
221  ofstream sol_ofs("sol.gf");
222  sol_ofs.precision(8);
223  x.Save(sol_ofs);
224  }
225 
226  // 13. Send the solution by socket to a GLVis server.
227  if (visualization)
228  {
229  char vishost[] = "localhost";
230  int visport = 19916;
231  socketstream sol_sock(vishost, visport);
232  sol_sock.precision(8);
233  sol_sock << "solution\n" << *mesh << x << flush;
234  }
235 
236  // 14. Free the used memory.
237  delete a;
238  delete b;
239  delete fespace;
240  delete mesh;
241 
242  return 0;
243 }
244 
245 double analytic_solution(const Vector &x)
246 {
247  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
248  return x(0)*x(1)/l2;
249 }
250 
251 double analytic_rhs(const Vector &x)
252 {
253  double l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
254  return 7*x(0)*x(1)/l2;
255 }
256 
257 void SnapNodes(Mesh &mesh)
258 {
259  GridFunction &nodes = *mesh.GetNodes();
260  Vector node(mesh.SpaceDimension());
261  for (int i = 0; i < nodes.FESpace()->GetNDofs(); i++)
262  {
263  for (int d = 0; d < mesh.SpaceDimension(); d++)
264  {
265  node(d) = nodes(nodes.FESpace()->DofToVDof(i, d));
266  }
267 
268  node /= node.Norml2();
269 
270  for (int d = 0; d < mesh.SpaceDimension(); d++)
271  {
272  nodes(nodes.FESpace()->DofToVDof(i, d)) = node(d);
273  }
274  }
275 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:162
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:116
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
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 RefineAtVertex(const Vertex &vert, int levels, double eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex, &#39;levels&#39; times.
Definition: mesh.cpp:6961
void AddVertex(const double *)
Definition: mesh.cpp:845
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:332
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2138
void FinalizeTriMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Definition: mesh.cpp:992
double ComputeL2Error(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:190
Data type sparse matrix.
Definition: sparsemat.hpp:38
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7316
double analytic_solution(const Vector &x)
Definition: ex7.cpp:245
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:440
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3719
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:270
int Dimension() const
Definition: mesh.hpp:475
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
double analytic_rhs(const Vector &x)
Definition: ex7.cpp:251
int SpaceDimension() const
Definition: mesh.hpp:476
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:8332
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator.
Definition: linearform.cpp:19
int main(int argc, char *argv[])
Definition: ex1.cpp:45
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:343
Abstract finite element space.
Definition: fespace.hpp:62
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:866
virtual int GetTrueVSize()
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:167
void RandomRefinement(int levels, int frac=2, bool aniso=false, int nonconforming=-1, int nc_limit=-1, int seed=0)
Refine each element with 1/frac probability, repeat &#39;levels&#39; times.
Definition: mesh.cpp:6938
void AddTriangle(const int *vi, int attr=1)
Definition: mesh.cpp:861
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Definition: mesh.cpp:1020
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:257
class for C-function coefficient
Vector data type.
Definition: vector.hpp:33
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5844
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:54
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:1725
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:1625
bool Good() const
Definition: optparser.hpp:120