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