MFEM  v3.0
 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-quad-nurbs.mesh
10 // ex2 -m ../data/beam-hex-nurbs.mesh
11 //
12 // Description: This example code solves a simple linear elasticity problem
13 // describing a multi-material cantilever beam.
14 //
15 // Specifically, we approximate the weak form of -div(sigma(u))=0
16 // where sigma(u)=lambda*div(u)*I+mu*(grad*u+u*grad) is the stress
17 // tensor corresponding to displacement field u, and lambda and mu
18 // are the material Lame constants. The boundary conditions are
19 // u=0 on the fixed part of the boundary with attribute 1, and
20 // sigma(u).n=f on the remainder with f being a constant pull down
21 // vector on boundary elements with attribute 2, and zero
22 // otherwise. The geometry of the domain is assumed to be as
23 // follows:
24 //
25 // +----------+----------+
26 // boundary --->| material | material |<--- boundary
27 // attribute 1 | 1 | 2 | attribute 2
28 // (fixed) +----------+----------+ (pull down)
29 //
30 // The example demonstrates the use of high-order and NURBS vector
31 // finite element spaces with the linear elasticity bilinear form,
32 // meshes with curved elements, and the definition of piece-wise
33 // constant and vector coefficient objects.
34 //
35 // We recommend viewing Example 1 before viewing this example.
36 
37 #include "mfem.hpp"
38 #include <fstream>
39 #include <iostream>
40 
41 using namespace std;
42 using namespace mfem;
43 
44 int main(int argc, char *argv[])
45 {
46  // 1. Parse command-line options.
47  const char *mesh_file = "../data/beam-tri.mesh";
48  int order = 1;
49  bool visualization = 1;
50 
51  OptionsParser args(argc, argv);
52  args.AddOption(&mesh_file, "-m", "--mesh",
53  "Mesh file to use.");
54  args.AddOption(&order, "-o", "--order",
55  "Finite element order (polynomial degree).");
56  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
57  "--no-visualization",
58  "Enable or disable GLVis visualization.");
59  args.Parse();
60  if (!args.Good())
61  {
62  args.PrintUsage(cout);
63  return 1;
64  }
65  args.PrintOptions(cout);
66 
67  // 2. Read the mesh from the given mesh file. We can handle triangular,
68  // quadrilateral, tetrahedral or hexahedral elements with the same code.
69  Mesh *mesh;
70  ifstream imesh(mesh_file);
71  if (!imesh)
72  {
73  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
74  return 2;
75  }
76  mesh = new Mesh(imesh, 1, 1);
77  imesh.close();
78  int dim = mesh->Dimension();
79 
80  if (mesh->attributes.Max() < 2 || mesh->bdr_attributes.Max() < 2)
81  {
82  cerr << "\nInput mesh should have at least two materials and "
83  << "two boundary attributes! (See schematic in ex2.cpp)\n"
84  << endl;
85  return 3;
86  }
87 
88  // 3. Select the order of the finite element discretization space. For NURBS
89  // meshes, we increase the order by degree elevation.
90  if (mesh->NURBSext && order > mesh->NURBSext->GetOrder())
91  mesh->DegreeElevate(order - mesh->NURBSext->GetOrder());
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  mesh->UniformRefinement();
102  }
103 
104  // 5. Define a finite element space on the mesh. Here we use vector finite
105  // elements, i.e. dim copies of a scalar finite element space. The vector
106  // dimension is specified by the last argument of the FiniteElementSpace
107  // constructor. For NURBS meshes, we use the (degree elevated) NURBS space
108  // associated with the mesh nodes.
110  FiniteElementSpace *fespace;
111  if (mesh->NURBSext)
112  {
113  fec = NULL;
114  fespace = mesh->GetNodes()->FESpace();
115  }
116  else
117  {
118  fec = new H1_FECollection(order, dim);
119  fespace = new FiniteElementSpace(mesh, fec, dim);
120  }
121  cout << "Number of unknowns: " << fespace->GetVSize() << endl
122  << "Assembling: " << flush;
123 
124  // 6. Set up the linear form b(.) which corresponds to the right-hand side of
125  // the FEM linear system. In this case, b_i equals the boundary integral
126  // of f*phi_i where f represents a "pull down" force on the Neumann part
127  // of the boundary and phi_i are the basis functions in the finite element
128  // fespace. The force is defined by the VectorArrayCoefficient object f,
129  // which is a vector of Coefficient objects. The fact that f is non-zero
130  // on boundary attribute 2 is indicated by the use of piece-wise constants
131  // coefficient for its last component.
132  VectorArrayCoefficient f(dim);
133  for (int i = 0; i < dim-1; i++)
134  f.Set(i, new ConstantCoefficient(0.0));
135  {
136  Vector pull_force(mesh->bdr_attributes.Max());
137  pull_force = 0.0;
138  pull_force(1) = -1.0e-2;
139  f.Set(dim-1, new PWConstCoefficient(pull_force));
140  }
141 
142  LinearForm *b = new LinearForm(fespace);
144  cout << "r.h.s. ... " << flush;
145  b->Assemble();
146 
147  // 7. Define the solution vector x as a finite element grid function
148  // corresponding to fespace. Initialize x with initial guess of zero,
149  // which satisfies the boundary conditions.
150  GridFunction x(fespace);
151  x = 0.0;
152 
153  // 8. Set up the bilinear form a(.,.) on the finite element space
154  // corresponding to the linear elasticity integrator with piece-wise
155  // constants coefficient lambda and mu. The boundary conditions are
156  // implemented by marking only boundary attribute 1 as essential. After
157  // assembly and finalizing we extract the corresponding sparse matrix A.
158  Vector lambda(mesh->attributes.Max());
159  lambda = 1.0;
160  lambda(0) = lambda(1)*50;
161  PWConstCoefficient lambda_func(lambda);
162  Vector mu(mesh->attributes.Max());
163  mu = 1.0;
164  mu(0) = mu(1)*50;
165  PWConstCoefficient mu_func(mu);
166 
167  BilinearForm *a = new BilinearForm(fespace);
168  a->AddDomainIntegrator(new ElasticityIntegrator(lambda_func,mu_func));
169  cout << "matrix ... " << flush;
170  a->Assemble();
171  Array<int> ess_bdr(mesh->bdr_attributes.Max());
172  ess_bdr = 0;
173  ess_bdr[0] = 1;
174  a->EliminateEssentialBC(ess_bdr, x, *b);
175  a->Finalize();
176  cout << "done." << endl;
177  const SparseMatrix &A = a->SpMat();
178 
179 #ifndef MFEM_USE_SUITESPARSE
180  // 9. Define a simple symmetric Gauss-Seidel preconditioner and use it to
181  // solve the system Ax=b with PCG.
182  GSSmoother M(A);
183  PCG(A, M, *b, x, 1, 500, 1e-8, 0.0);
184 #else
185  // 9. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
186  UMFPackSolver umf_solver;
187  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
188  umf_solver.SetOperator(A);
189  umf_solver.Mult(*b, x);
190 #endif
191 
192  // 10. For non-NURBS meshes, make the mesh curved based on the finite element
193  // space. This means that we define the mesh elements through a fespace
194  // based transformation of the reference element. This allows us to save
195  // the displaced mesh as a curved mesh when using high-order finite
196  // element displacement field. We assume that the initial mesh (read from
197  // the file) is not higher order curved mesh compared to the chosen FE
198  // space.
199  if (!mesh->NURBSext)
200  mesh->SetNodalFESpace(fespace);
201 
202  // 11. Save the displaced mesh and the inverted solution (which gives the
203  // backward displacements to the original grid). This output can be
204  // viewed later using GLVis: "glvis -m displaced.mesh -g sol.gf".
205  {
206  GridFunction *nodes = mesh->GetNodes();
207  *nodes += x;
208  x *= -1;
209  ofstream mesh_ofs("displaced.mesh");
210  mesh_ofs.precision(8);
211  mesh->Print(mesh_ofs);
212  ofstream sol_ofs("sol.gf");
213  sol_ofs.precision(8);
214  x.Save(sol_ofs);
215  }
216 
217  // 12. Send the above data by socket to a GLVis server. Use the "n" and "b"
218  // keys in GLVis to visualize the displacements.
219  if (visualization)
220  {
221  char vishost[] = "localhost";
222  int visport = 19916;
223  socketstream sol_sock(vishost, visport);
224  sol_sock.precision(8);
225  sol_sock << "solution\n" << *mesh << x << flush;
226  }
227 
228  // 13. Free the used memory.
229  delete a;
230  delete b;
231  if (fec)
232  {
233  delete fespace;
234  delete fec;
235  }
236  delete mesh;
237 
238  return 0;
239 }
int GetVSize() const
Definition: fespace.hpp:151
Vector coefficient defined by an array of scalar coefficients.
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:26
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:2878
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:396
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:329
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:1797
Data type sparse matrix.
Definition: sparsemat.hpp:38
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6225
T Max() const
Definition: array.cpp:78
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:424
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3066
void AddBoundaryIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Integrator.
Definition: linearform.cpp:24
int Dimension() const
Definition: mesh.hpp:417
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:385
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:7136
Array< int > bdr_attributes
Definition: mesh.hpp:305
int main(int argc, char *argv[])
Definition: ex1.cpp:39
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:340
Abstract finite element space.
Definition: fespace.hpp:61
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:307
class for piecewise constant coefficient
Definition: coefficient.hpp:72
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:266
Vector data type.
Definition: vector.hpp:29
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:4948
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:52
Class for linear form - Vector with associated FE space and LFIntegrators.
Definition: linearform.hpp:23
const SparseMatrix & SpMat() const
Returns a reference to the sparse martix.
void Set(int i, Coefficient *c)
Sets coefficient in the vector.
void EliminateEssentialBC(Array< int > &bdr_attr_is_ess, Vector &sol, Vector &rhs, int d=0)
Array< int > attributes
Definition: mesh.hpp:304
virtual void Mult(const Vector &b, Vector &x) const
Operator application.
Definition: solvers.cpp:1597
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:1505
bool Good() const
Definition: optparser.hpp:120