MFEM  v3.1 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
7 // ex2 -m ../data/beam-tet.mesh
8 // ex2 -m ../data/beam-hex.mesh
9 // ex2 -m ../data/beam-quad.mesh -o 3 -sc
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
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);
56  "Mesh file to use.");
58  "Finite element order (polynomial degree).");
60  "--no-static-condensation", "Enable static condensation.");
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;
75  ifstream imesh(mesh_file);
76  if (!imesh)
77  {
78  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
79  return 2;
80  }
81  mesh = new Mesh(imesh, 1, 1);
82  imesh.close();
83  int dim = mesh->Dimension();
84
85  if (mesh->attributes.Max() < 2 || mesh->bdr_attributes.Max() < 2)
86  {
87  cerr << "\nInput mesh should have at least two materials and "
88  << "two boundary attributes! (See schematic in ex2.cpp)\n"
89  << endl;
90  return 3;
91  }
92
93  // 3. Select the order of the finite element discretization space. For NURBS
94  // meshes, we increase the order by degree elevation.
95  if (mesh->NURBSext && order > mesh->NURBSext->GetOrder())
96  {
97  mesh->DegreeElevate(order - mesh->NURBSext->GetOrder());
98  }
99
100  // 4. Refine the mesh to increase the resolution. In this example we do
101  // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
102  // largest number that gives a final mesh with no more than 5,000
103  // elements.
104  {
105  int ref_levels =
106  (int)floor(log(5000./mesh->GetNE())/log(2.)/dim);
107  for (int l = 0; l < ref_levels; l++)
108  {
109  mesh->UniformRefinement();
110  }
111  }
112
113  // 5. Define a finite element space on the mesh. Here we use vector finite
114  // elements, i.e. dim copies of a scalar finite element space. The vector
115  // dimension is specified by the last argument of the FiniteElementSpace
116  // constructor. For NURBS meshes, we use the (degree elevated) NURBS space
117  // associated with the mesh nodes.
119  FiniteElementSpace *fespace;
120  if (mesh->NURBSext)
121  {
122  fec = NULL;
123  fespace = mesh->GetNodes()->FESpace();
124  }
125  else
126  {
127  fec = new H1_FECollection(order, dim);
128  fespace = new FiniteElementSpace(mesh, fec, dim);
129  }
130  cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
131  << endl << "Assembling: " << flush;
132
133  // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
134  // In this example, the boundary conditions are defined by marking only
135  // boundary attribute 1 from the mesh as essential and converting it to a
136  // list of true dofs.
137  Array<int> ess_tdof_list, ess_bdr(mesh->bdr_attributes.Max());
138  ess_bdr = 0;
139  ess_bdr[0] = 1;
140  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
141
142  // 7. Set up the linear form b(.) which corresponds to the right-hand side of
143  // the FEM linear system. In this case, b_i equals the boundary integral
144  // of f*phi_i where f represents a "pull down" force on the Neumann part
145  // of the boundary and phi_i are the basis functions in the finite element
146  // fespace. The force is defined by the VectorArrayCoefficient object f,
147  // which is a vector of Coefficient objects. The fact that f is non-zero
148  // on boundary attribute 2 is indicated by the use of piece-wise constants
149  // coefficient for its last component.
150  VectorArrayCoefficient f(dim);
151  for (int i = 0; i < dim-1; i++)
152  {
153  f.Set(i, new ConstantCoefficient(0.0));
154  }
155  {
156  Vector pull_force(mesh->bdr_attributes.Max());
157  pull_force = 0.0;
158  pull_force(1) = -1.0e-2;
159  f.Set(dim-1, new PWConstCoefficient(pull_force));
160  }
161
162  LinearForm *b = new LinearForm(fespace);
164  cout << "r.h.s. ... " << flush;
165  b->Assemble();
166
167  // 8. Define the solution vector x as a finite element grid function
168  // corresponding to fespace. Initialize x with initial guess of zero,
169  // which satisfies the boundary conditions.
170  GridFunction x(fespace);
171  x = 0.0;
172
173  // 9. Set up the bilinear form a(.,.) on the finite element space
174  // corresponding to the linear elasticity integrator with piece-wise
175  // constants coefficient lambda and mu.
176  Vector lambda(mesh->attributes.Max());
177  lambda = 1.0;
178  lambda(0) = lambda(1)*50;
179  PWConstCoefficient lambda_func(lambda);
180  Vector mu(mesh->attributes.Max());
181  mu = 1.0;
182  mu(0) = mu(1)*50;
183  PWConstCoefficient mu_func(mu);
184
185  BilinearForm *a = new BilinearForm(fespace);
187
188  // 10. Assemble the bilinear form and the corresponding linear system,
189  // applying any necessary transformations such as: eliminating boundary
190  // conditions, applying conforming constraints for non-conforming AMR,
191  // static condensation, etc.
192  cout << "matrix ... " << flush;
193  if (static_cond) { a->EnableStaticCondensation(); }
194  a->Assemble();
195
196  SparseMatrix A;
197  Vector B, X;
198  a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
199  cout << "done." << endl;
200
201  cout << "Size of linear system: " << A.Height() << endl;
202
203 #ifndef MFEM_USE_SUITESPARSE
204  // 11. 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, 500, 1e-8, 0.0);
208 #else
209  // 11. 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  // 12. Recover the solution as a finite element grid function.
217  a->RecoverFEMSolution(X, *b, x);
218
219  // 13. For non-NURBS meshes, make the mesh curved based on the finite element
220  // space. This means that we define the mesh elements through a fespace
221  // based transformation of the reference element. This allows us to save
222  // the displaced mesh as a curved mesh when using high-order finite
223  // element displacement field. We assume that the initial mesh (read from
224  // the file) is not higher order curved mesh compared to the chosen FE
225  // space.
226  if (!mesh->NURBSext)
227  {
228  mesh->SetNodalFESpace(fespace);
229  }
230
231  // 14. Save the displaced mesh and the inverted solution (which gives the
232  // backward displacements to the original grid). This output can be
233  // viewed later using GLVis: "glvis -m displaced.mesh -g sol.gf".
234  {
235  GridFunction *nodes = mesh->GetNodes();
236  *nodes += x;
237  x *= -1;
238  ofstream mesh_ofs("displaced.mesh");
239  mesh_ofs.precision(8);
240  mesh->Print(mesh_ofs);
241  ofstream sol_ofs("sol.gf");
242  sol_ofs.precision(8);
243  x.Save(sol_ofs);
244  }
245
246  // 15. Send the above data by socket to a GLVis server. Use the "n" and "b"
247  // keys in GLVis to visualize the displacements.
248  if (visualization)
249  {
250  char vishost[] = "localhost";
251  int visport = 19916;
252  socketstream sol_sock(vishost, visport);
253  sol_sock.precision(8);
254  sol_sock << "solution\n" << *mesh << x << flush;
255  }
256
257  // 16. Free the used memory.
258  delete a;
259  delete b;
260  if (fec)
261  {
262  delete fespace;
263  delete fec;
264  }
265  delete mesh;
266
267  return 0;
268 }
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:34
void DegreeElevate(int t)
Definition: mesh.cpp:3497
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list)
Definition: fespace.cpp:469
void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:454
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.
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
int dim
Definition: ex3.cpp:48
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:35
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7316
T Max() const
Definition: array.cpp:90
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
Definition: linearform.cpp:24
int Dimension() const
Definition: mesh.hpp:475
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
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
Array< int > bdr_attributes
Definition: mesh.hpp:140
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
NURBSExtension * NURBSext
Definition: mesh.hpp:142
class for piecewise constant coefficient
Definition: coefficient.hpp:72
virtual int GetTrueVSize()
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:167
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
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
void Set(int i, Coefficient *c)
Sets coefficient in the vector.
Array< int > attributes
Definition: mesh.hpp:139
virtual void Mult(const Vector &b, Vector &x) const
Operator application.
Definition: solvers.cpp:1725
void EnableStaticCondensation()
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