MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
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
6 // ex2 -m ../data/beam-quad.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
11 // ex2 -m ../data/beam-quad-nurbs.mesh
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
18 // where sigma(u)=lambda*div(u)*I+mu*(grad*u+u*grad) is the stress
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 
44 using namespace std;
45 using namespace mfem;
46 
47 int 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);
56  args.AddOption(&mesh_file, "-m", "--mesh",
57  "Mesh file to use.");
58  args.AddOption(&order, "-o", "--order",
59  "Finite element order (polynomial degree).");
60  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
61  "--no-static-condensation", "Enable static condensation.");
62  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
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);
179  a->AddDomainIntegrator(new ElasticityIntegrator(lambda_func,mu_func));
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 }
Vector coefficient defined by an array of scalar coefficients. Coefficients that are not set will eva...
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:168
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition: fespace.cpp:564
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:1025
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3676
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
Data type sparse matrix.
Definition: sparsemat.hpp:50
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
constexpr char vishost[]
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9816
constexpr int visport
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void Set(int i, Coefficient *c, bool own=true)
Sets coefficient in the vector.
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:904
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:590
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:5285
void AddBoundaryIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:72
int Dimension() const
Definition: mesh.hpp:1006
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:1036
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
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 &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
A piecewise constant coefficient with the constants keyed off the element attribute numbers...
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1671
double mu
Definition: ex25.cpp:139
void DegreeElevate(int rel_degree, int degree=16)
Definition: mesh.cpp:5107
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7908
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
int main()
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:268
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3189
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...
double f(const Vector &p)
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:3094
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150