MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex28.cpp
Go to the documentation of this file.
1 // MFEM Example 28
2 //
3 // Compile with: make ex28
4 //
5 // Sample runs: ex28
6 // ex28 --visit-datafiles
7 // ex28 --order 2
8 //
9 // Description: Demonstrates a sliding boundary condition in an elasticity
10 // problem. A trapezoid, roughly as pictured below, is pushed
11 // from the right into a rigid notch. Normal displacement is
12 // restricted, but tangential movement is allowed, so the
13 // trapezoid compresses into the notch.
14 //
15 // /-------+
16 // normal constrained --->/ | <--- boundary force (2)
17 // boundary (4) /---------+
18 // ^
19 // |
20 // normal constrained boundary (1)
21 //
22 // This example demonstrates the use of the ConstrainedSolver
23 // framework.
24 //
25 // We recommend viewing Example 2 before viewing this example.
26 
27 #include "mfem.hpp"
28 #include <fstream>
29 #include <iostream>
30 #include <set>
31 
32 using namespace std;
33 using namespace mfem;
34 
35 // Return a mesh with a single element with vertices (0, 0), (1, 0), (1, 1),
36 // (offset, 1) to demonstrate boundary conditions on a surface that is not
37 // axis-aligned.
38 Mesh * build_trapezoid_mesh(double offset)
39 {
40  MFEM_VERIFY(offset < 0.9, "offset is too large!");
41 
42  const int dimension = 2;
43  const int nvt = 4; // vertices
44  const int nbe = 4; // num boundary elements
45  Mesh * mesh = new Mesh(dimension, nvt, 1, nbe);
46 
47  // vertices
48  double vc[dimension];
49  vc[0] = 0.0; vc[1] = 0.0;
50  mesh->AddVertex(vc);
51  vc[0] = 1.0; vc[1] = 0.0;
52  mesh->AddVertex(vc);
53  vc[0] = offset; vc[1] = 1.0;
54  mesh->AddVertex(vc);
55  vc[0] = 1.0; vc[1] = 1.0;
56  mesh->AddVertex(vc);
57 
58  // element
59  Array<int> vert(4);
60  vert[0] = 0; vert[1] = 1; vert[2] = 3; vert[3] = 2;
61  mesh->AddQuad(vert, 1);
62 
63  // boundary
64  Array<int> sv(2);
65  sv[0] = 0; sv[1] = 1;
66  mesh->AddBdrSegment(sv, 1);
67  sv[0] = 1; sv[1] = 3;
68  mesh->AddBdrSegment(sv, 2);
69  sv[0] = 2; sv[1] = 3;
70  mesh->AddBdrSegment(sv, 3);
71  sv[0] = 0; sv[1] = 2;
72  mesh->AddBdrSegment(sv, 4);
73 
74  mesh->FinalizeQuadMesh(1, 0, true);
75 
76  return mesh;
77 }
78 
79 int main(int argc, char *argv[])
80 {
81  // 1. Parse command-line options.
82  int order = 1;
83  bool visualization = 1;
84  double offset = 0.3;
85  bool visit = false;
86 
87  OptionsParser args(argc, argv);
88  args.AddOption(&order, "-o", "--order",
89  "Finite element order (polynomial degree).");
90  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
91  "--no-visualization",
92  "Enable or disable GLVis visualization.");
93  args.AddOption(&offset, "--offset", "--offset",
94  "How much to offset the trapezoid.");
95  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
96  "--no-visit-datafiles",
97  "Save data files for VisIt (visit.llnl.gov) visualization.");
98  args.Parse();
99  if (!args.Good())
100  {
101  args.PrintUsage(cout);
102  return 1;
103  }
104  args.PrintOptions(cout);
105 
106  // 2. Build a trapezoidal mesh with a single quadrilateral element, where
107  // 'offset' determines how far off it is from a rectangle.
108  Mesh *mesh = build_trapezoid_mesh(offset);
109  int dim = mesh->Dimension();
110 
111  // 3. Refine the mesh to increase the resolution. In this example we do
112  // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
113  // largest number that gives a final mesh with no more than 1,000
114  // elements.
115  {
116  int ref_levels =
117  (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
118  for (int l = 0; l < ref_levels; l++)
119  {
120  mesh->UniformRefinement();
121  }
122  }
123 
124  // 4. Define a finite element space on the mesh. Here we use vector finite
125  // elements, i.e. dim copies of a scalar finite element space. The vector
126  // dimension is specified by the last argument of the FiniteElementSpace
127  // constructor.
128  FiniteElementCollection *fec = new H1_FECollection(order, dim);
129  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec, dim);
130  cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
131  << endl;
132  cout << "Assembling matrix and r.h.s... " << flush;
133 
134  // 5. Determine the list of true (i.e. parallel conforming) essential
135  // boundary dofs. In this example, there are no essential boundary
136  // conditions in the usual sense, but we leave the machinery here for
137  // users to modify if they wish.
138  Array<int> ess_tdof_list, ess_bdr(mesh->bdr_attributes.Max());
139  ess_bdr = 0;
140  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
141 
142  // 6. 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 "push" force on the right side of the
145  // trapezoid.
147  for (int i = 0; i < dim-1; i++)
148  {
149  f.Set(i, new ConstantCoefficient(0.0));
150  }
151  {
152  Vector push_force(mesh->bdr_attributes.Max());
153  push_force = 0.0;
154  push_force(1) = -5.0e-2; // index 1 attribute 2
155  f.Set(0, new PWConstCoefficient(push_force));
156  }
157  LinearForm *b = new LinearForm(fespace);
159  b->Assemble();
160 
161  // 7. Define the solution vector x as a finite element grid function
162  // corresponding to fespace.
163  GridFunction x(fespace);
164  x = 0.0;
165 
166  // 8. 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. We use constant coefficients,
169  // but see ex2 for how to set up piecewise constant coefficients based
170  // on attribute.
171  Vector lambda(mesh->attributes.Max());
172  lambda = 1.0;
173  PWConstCoefficient lambda_func(lambda);
174  Vector mu(mesh->attributes.Max());
175  mu = 1.0;
176  PWConstCoefficient mu_func(mu);
177 
178  BilinearForm *a = new BilinearForm(fespace);
179  a->AddDomainIntegrator(new ElasticityIntegrator(lambda_func, mu_func));
180 
181  // 9. 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  a->Assemble();
186 
187  SparseMatrix A;
188  Vector B, X;
189  a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
190  cout << "done." << endl;
191  cout << "Size of linear system: " << A.Height() << endl;
192 
193  // 10. Set up constraint matrix to constrain normal displacement (but
194  // allow tangential displacement) on specified boundaries.
195  Array<int> constraint_atts(2);
196  constraint_atts[0] = 1; // attribute 1 bottom
197  constraint_atts[1] = 4; // attribute 4 left side
198  Array<int> lagrange_rowstarts;
199  SparseMatrix* local_constraints =
200  BuildNormalConstraints(*fespace, constraint_atts, lagrange_rowstarts);
201 
202  // 11. Define and apply an iterative solver for the constrained system
203  // in saddle-point form with a Gauss-Seidel smoother for the
204  // displacement block.
205  GSSmoother M(A);
206  SchurConstrainedSolver * solver =
207  new SchurConstrainedSolver(A, *local_constraints, M);
208  solver->SetRelTol(1e-5);
209  solver->SetMaxIter(2000);
210  solver->SetPrintLevel(1);
211  solver->Mult(B, X);
212 
213  // 12. Recover the solution as a finite element grid function. Move the
214  // mesh to reflect the displacement of the elastic body being
215  // simulated, for purposes of output.
216  a->RecoverFEMSolution(X, *b, x);
217  mesh->SetNodalFESpace(fespace);
218  GridFunction *nodes = mesh->GetNodes();
219  *nodes += x;
220 
221  // 13. Save the refined mesh and the solution in VisIt format.
222  if (visit)
223  {
224  VisItDataCollection visit_dc("ex28", mesh);
225  visit_dc.SetLevelsOfDetail(4);
226  visit_dc.RegisterField("displacement", &x);
227  visit_dc.Save();
228  }
229 
230  // 14. Save the displaced mesh and the inverted solution (which gives the
231  // backward displacements to the original grid). This output can be
232  // viewed later using GLVis: "glvis -m displaced.mesh -g sol.gf".
233  {
234  x *= -1; // sign convention for GLVis displacements
235  ofstream mesh_ofs("displaced.mesh");
236  mesh_ofs.precision(8);
237  mesh->Print(mesh_ofs);
238  ofstream sol_ofs("sol.gf");
239  sol_ofs.precision(8);
240  x.Save(sol_ofs);
241  }
242 
243  // 15. Send the above data by socket to a GLVis server. Use the "n" and "b"
244  // keys in GLVis to visualize the displacements.
245  if (visualization)
246  {
247  char vishost[] = "localhost";
248  int visport = 19916;
249  socketstream sol_sock(vishost, visport);
250  sol_sock.precision(8);
251  sol_sock << "solution\n" << *mesh << x << flush;
252  }
253 
254  // 16. Free the used memory.
255  delete local_constraints;
256  delete solver;
257  delete a;
258  delete b;
259  if (fec)
260  {
261  delete fespace;
262  delete fec;
263  }
264  delete mesh;
265 
266  return 0;
267 }
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
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Definition: mesh.cpp:1674
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
Solve constrained system by solving original mixed system; see ConstrainedSolver. ...
SparseMatrix * BuildNormalConstraints(FiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts, bool parallel)
Build a matrix constraining normal components to zero.
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.
virtual void Mult(const Vector &f, Vector &x) const override
Solve for given .
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3676
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1613
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
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
Data collection with VisIt I/O routines.
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
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.
int AddBdrSegment(int v1, int v2, int attr=1)
Definition: mesh.cpp:1823
void SetLevelsOfDetail(int levels_of_detail)
Set VisIt parameter: default levels of detail for the MultiresControl.
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
void SetRelTol(double rtol)
Definition: solvers.hpp:198
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
virtual void Save() override
Save the collection and a VisIt root file.
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
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
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition: mesh.cpp:1968
double mu
Definition: ex25.cpp:139
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
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
Definition: hooke.cpp:45
Mesh * build_trapezoid_mesh(double offset)
Definition: ex28.cpp:38
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
double f(const Vector &p)
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150