MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
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
27using namespace std;
28using namespace mfem;
29
30// Exact solution and r.h.s., see below for implementation.
32real_t analytic_rhs(const Vector &x);
33void SnapNodes(Mesh &mesh);
34
35int 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 real_t 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 real_t 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);
193 a->AddDomainIntegrator(new DiffusionIntegrator(one));
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
253{
254 real_t l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
255 return x(0)*x(1)/l2;
256}
257
259{
260 real_t l2 = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
261 return 7*x(0)*x(1)/l2;
262}
263
264void 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 if (mesh.Nonconforming())
283 {
284 // Snap hanging nodes to the master side.
285 Vector tnodes;
286 nodes.GetTrueDofs(tnodes);
287 nodes.SetFromTrueDofs(tnodes);
288 }
289}
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
A coefficient that is constant across space and time.
Class for domain integration .
Definition lininteg.hpp:109
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Definition fespace.cpp:249
A general function coefficient.
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
virtual real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition gridfunc.cpp:381
FiniteElementSpace * FESpace()
Definition gridfunc.hpp:696
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction.
Definition gridfunc.cpp:366
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Vector with associated FE space and LinearFormIntegrators.
Mesh data type.
Definition mesh.hpp:56
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Adds a quadrilateral to the mesh given by 4 vertices v1 through v4.
Definition mesh.cpp:1743
bool Nonconforming() const
Definition mesh.hpp:2229
int AddTriangle(int v1, int v2, int v3, int attr=1)
Adds a triangle to the mesh given by 3 vertices v1 through v3.
Definition mesh.cpp:1729
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
Definition mesh.cpp:1658
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void RandomRefinement(real_t prob, bool aniso=false, int nonconforming=-1, int nc_limit=0)
Refine each element with given probability. Uses GeneralRefinement.
Definition mesh.cpp:10650
void RefineAtVertex(const Vertex &vert, real_t eps=0.0, int nonconforming=-1)
Refine elements sharing the specified vertex. Uses GeneralRefinement.
Definition mesh.cpp:10669
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition mesh.cpp:2177
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
void FinalizeTriMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a triangular Mesh.
Definition mesh.cpp:2148
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition mesh.cpp:6153
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Data type sparse matrix.
Definition sparsemat.hpp:51
Direct sparse solver using UMFPACK.
Definition solvers.hpp:1096
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition solvers.cpp:3102
real_t Control[UMFPACK_CONTROL]
Definition solvers.hpp:1106
virtual void Mult(const Vector &b, Vector &x) const
Direct solution of the linear system using UMFPACK.
Definition solvers.cpp:3197
Vector data type.
Definition vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
Data type for vertex.
Definition vertex.hpp:23
void SnapNodes(Mesh &mesh)
Definition ex7.cpp:264
real_t analytic_rhs(const Vector &x)
Definition ex7.cpp:258
real_t analytic_solution(const Vector &x)
Definition ex7.cpp:252
int main()
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition solvers.cpp:913
float real_t
Definition config.hpp:43
const char vishost[]