MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex2.cpp
Go to the documentation of this file.
1// MFEM Example 2
2//
3// Compile with: make ex2
4//
5// Sample runs: ex2 -m ../data/beam-tri.mesh
7// ex2 -m ../data/beam-tet.mesh
8// ex2 -m ../data/beam-hex.mesh
9// ex2 -m ../data/beam-wedge.mesh
10// ex2 -m ../data/beam-quad.mesh -o 3 -sc
12// ex2 -m ../data/beam-hex-nurbs.mesh
13//
14// Description: This example code solves a simple linear elasticity problem
15// describing a multi-material cantilever beam.
16//
17// Specifically, we approximate the weak form of -div(sigma(u))=0
19// tensor corresponding to displacement field u, and lambda and mu
20// are the material Lame constants. The boundary conditions are
21// u=0 on the fixed part of the boundary with attribute 1, and
22// sigma(u).n=f on the remainder with f being a constant pull down
23// vector on boundary elements with attribute 2, and zero
24// otherwise. The geometry of the domain is assumed to be as
25// follows:
26//
27// +----------+----------+
28// boundary --->| material | material |<--- boundary
29// attribute 1 | 1 | 2 | attribute 2
30// (fixed) +----------+----------+ (pull down)
31//
32// The example demonstrates the use of high-order and NURBS vector
33// finite element spaces with the linear elasticity bilinear form,
34// meshes with curved elements, and the definition of piece-wise
35// constant and vector coefficient objects. Static condensation is
36// also illustrated.
37//
38// We recommend viewing Example 1 before viewing this example.
39
40#include "mfem.hpp"
41#include <fstream>
42#include <iostream>
43
44using namespace std;
45using namespace mfem;
46
47int main(int argc, char *argv[])
48{
49 // 1. Parse command-line options.
50 const char *mesh_file = "../data/beam-tri.mesh";
51 int order = 1;
52 bool static_cond = false;
53 bool visualization = 1;
54
55 OptionsParser args(argc, argv);
57 "Mesh file to use.");
59 "Finite element order (polynomial degree).");
61 "--no-static-condensation", "Enable static condensation.");
63 "--no-visualization",
64 "Enable or disable GLVis visualization.");
65 args.Parse();
66 if (!args.Good())
67 {
68 args.PrintUsage(cout);
69 return 1;
70 }
71 args.PrintOptions(cout);
72
73 // 2. Read the mesh from the given mesh file. We can handle triangular,
74 // quadrilateral, tetrahedral or hexahedral elements with the same code.
75 Mesh *mesh = new Mesh(mesh_file, 1, 1);
76 int dim = mesh->Dimension();
77
78 if (mesh->attributes.Max() < 2 || mesh->bdr_attributes.Max() < 2)
79 {
80 cerr << "\nInput mesh should have at least two materials and "
81 << "two boundary attributes! (See schematic in ex2.cpp)\n"
82 << endl;
83 return 3;
84 }
85
86 // 3. Select the order of the finite element discretization space. For NURBS
87 // meshes, we increase the order by degree elevation.
88 if (mesh->NURBSext)
89 {
90 mesh->DegreeElevate(order, order);
91 }
92
93 // 4. Refine the mesh to increase the resolution. In this example we do
94 // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
95 // largest number that gives a final mesh with no more than 5,000
96 // elements.
97 {
98 int ref_levels =
99 (int)floor(log(5000./mesh->GetNE())/log(2.)/dim);
100 for (int l = 0; l < ref_levels; l++)
101 {
102 mesh->UniformRefinement();
103 }
104 }
105
106 // 5. Define a finite element space on the mesh. Here we use vector finite
107 // elements, i.e. dim copies of a scalar finite element space. The vector
108 // dimension is specified by the last argument of the FiniteElementSpace
109 // constructor. For NURBS meshes, we use the (degree elevated) NURBS space
110 // associated with the mesh nodes.
112 FiniteElementSpace *fespace;
113 if (mesh->NURBSext)
114 {
115 fec = NULL;
116 fespace = mesh->GetNodes()->FESpace();
117 }
118 else
119 {
120 fec = new H1_FECollection(order, dim);
121 fespace = new FiniteElementSpace(mesh, fec, dim);
122 }
123 cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
124 << endl << "Assembling: " << flush;
125
126 // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
127 // In this example, the boundary conditions are defined by marking only
128 // boundary attribute 1 from the mesh as essential and converting it to a
129 // list of true dofs.
130 Array<int> ess_tdof_list, ess_bdr(mesh->bdr_attributes.Max());
131 ess_bdr = 0;
132 ess_bdr[0] = 1;
133 fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
134
135 // 7. Set up the linear form b(.) which corresponds to the right-hand side of
136 // the FEM linear system. In this case, b_i equals the boundary integral
137 // of f*phi_i where f represents a "pull down" force on the Neumann part
138 // of the boundary and phi_i are the basis functions in the finite element
139 // fespace. The force is defined by the VectorArrayCoefficient object f,
140 // which is a vector of Coefficient objects. The fact that f is non-zero
141 // on boundary attribute 2 is indicated by the use of piece-wise constants
142 // coefficient for its last component.
144 for (int i = 0; i < dim-1; i++)
145 {
146 f.Set(i, new ConstantCoefficient(0.0));
147 }
148 {
149 Vector pull_force(mesh->bdr_attributes.Max());
150 pull_force = 0.0;
151 pull_force(1) = -1.0e-2;
152 f.Set(dim-1, new PWConstCoefficient(pull_force));
153 }
154
155 LinearForm *b = new LinearForm(fespace);
157 cout << "r.h.s. ... " << flush;
158 b->Assemble();
159
160 // 8. Define the solution vector x as a finite element grid function
161 // corresponding to fespace. Initialize x with initial guess of zero,
162 // which satisfies the boundary conditions.
163 GridFunction x(fespace);
164 x = 0.0;
165
166 // 9. Set up the bilinear form a(.,.) on the finite element space
167 // corresponding to the linear elasticity integrator with piece-wise
168 // constants coefficient lambda and mu.
169 Vector lambda(mesh->attributes.Max());
170 lambda = 1.0;
171 lambda(0) = lambda(1)*50;
172 PWConstCoefficient lambda_func(lambda);
173 Vector mu(mesh->attributes.Max());
174 mu = 1.0;
175 mu(0) = mu(1)*50;
176 PWConstCoefficient mu_func(mu);
177
178 BilinearForm *a = new BilinearForm(fespace);
180
181 // 10. Assemble the bilinear form and the corresponding linear system,
182 // applying any necessary transformations such as: eliminating boundary
183 // conditions, applying conforming constraints for non-conforming AMR,
184 // static condensation, etc.
185 cout << "matrix ... " << flush;
186 if (static_cond) { a->EnableStaticCondensation(); }
187 a->Assemble();
188
189 SparseMatrix A;
190 Vector B, X;
191 a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
192 cout << "done." << endl;
193
194 cout << "Size of linear system: " << A.Height() << endl;
195
196#ifndef MFEM_USE_SUITESPARSE
197 // 11. 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, 500, 1e-8, 0.0);
201#else
202 // 11. 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 // 12. Recover the solution as a finite element grid function.
210 a->RecoverFEMSolution(X, *b, x);
211
212 // 13. For non-NURBS meshes, make the mesh curved based on the finite element
213 // space. This means that we define the mesh elements through a fespace
214 // based transformation of the reference element. This allows us to save
215 // the displaced mesh as a curved mesh when using high-order finite
216 // element displacement field. We assume that the initial mesh (read from
217 // the file) is not higher order curved mesh compared to the chosen FE
218 // space.
219 if (!mesh->NURBSext)
220 {
221 mesh->SetNodalFESpace(fespace);
222 }
223
224 // 14. Save the displaced mesh and the inverted solution (which gives the
225 // backward displacements to the original grid). This output can be
226 // viewed later using GLVis: "glvis -m displaced.mesh -g sol.gf".
227 {
228 GridFunction *nodes = mesh->GetNodes();
229 *nodes += x;
230 x *= -1;
231 ofstream mesh_ofs("displaced.mesh");
232 mesh_ofs.precision(8);
233 mesh->Print(mesh_ofs);
234 ofstream sol_ofs("sol.gf");
235 sol_ofs.precision(8);
236 x.Save(sol_ofs);
237 }
238
239 // 15. Send the above data by socket to a GLVis server. Use the "n" and "b"
240 // keys in GLVis to visualize the displacements.
241 if (visualization)
242 {
243 char vishost[] = "localhost";
244 int visport = 19916;
245 socketstream sol_sock(vishost, visport);
246 sol_sock.precision(8);
247 sol_sock << "solution\n" << *mesh << x << flush;
248 }
249
250 // 16. Free the used memory.
251 delete a;
252 delete b;
253 if (fec)
254 {
255 delete fespace;
256 delete fec;
257 }
258 delete mesh;
259
260 return 0;
261}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
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.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
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
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition fespace.cpp:588
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.
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
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:290
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
void DegreeElevate(int rel_degree, int degree=16)
Definition mesh.cpp:5779
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition mesh.cpp:6153
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
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.
A piecewise constant coefficient with the constants keyed off the element attribute numbers.
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 coefficient defined by an array of scalar coefficients. Coefficients that are not set will eva...
Vector data type.
Definition vector.hpp:80
Vector & Set(const real_t a, const Vector &x)
(*this) = a * x
Definition vector.cpp:262
int dim
Definition ex24.cpp:53
real_t mu
Definition ex25.cpp:140
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
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
const char vishost[]