MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
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
32using namespace std;
33using 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.
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 real_t vc[dimension];
49 vc[0] = 0.0; vc[1] = 0.0;
51 vc[0] = 1.0; vc[1] = 0.0;
53 vc[0] = offset; vc[1] = 1.0;
55 vc[0] = 1.0; vc[1] = 1.0;
57
58 // element
59 Array<int> vert(4);
60 vert[0] = 0; vert[1] = 1; vert[2] = 3; vert[3] = 2;
62
63 // boundary
64 Array<int> sv(2);
65 sv[0] = 0; sv[1] = 1;
67 sv[0] = 1; sv[1] = 3;
69 sv[0] = 2; sv[1] = 3;
71 sv[0] = 0; sv[1] = 2;
73
75
76 return mesh;
77}
78
79int main(int argc, char *argv[])
80{
81 // 1. Parse command-line options.
82 int order = 1;
83 bool visualization = 1;
84 real_t offset = 0.3;
85 bool visit = false;
86
87 OptionsParser args(argc, argv);
89 "Finite element order (polynomial degree).");
91 "--no-visualization",
92 "Enable or disable GLVis visualization.");
94 "How much to offset the trapezoid.");
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.
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);
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}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
A coefficient that is constant across space and time.
virtual void Mult(const Vector &f, Vector &x) const override
Solve for given .
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition fespace.cpp:588
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
Vector with associated FE space and LinearFormIntegrators.
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
Adds a quadrilateral to the mesh given by 4 vertices v1 through v4.
Definition mesh.cpp:1743
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
Definition mesh.cpp:1658
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
int AddBdrSegment(int v1, int v2, int attr=1)
Definition mesh.cpp:2035
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition mesh.cpp:2177
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition mesh.cpp:6153
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
A piecewise constant coefficient with the constants keyed off the element attribute numbers.
Solve constrained system by solving original mixed system; see ConstrainedSolver.
Data type sparse matrix.
Definition sparsemat.hpp:51
Vector coefficient defined by an array of scalar coefficients. Coefficients that are not set will eva...
Vector data type.
Definition vector.hpp:80
Vector & Set(const real_t a, const Vector &x)
(*this) = a * x
Definition vector.cpp:262
Data collection with VisIt I/O routines.
void SetLevelsOfDetail(int levels_of_detail)
Set VisIt parameter: default levels of detail for the MultiresControl.
virtual void Save() override
Save the collection and a VisIt root file.
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
int dim
Definition ex24.cpp:53
real_t mu
Definition ex25.cpp:140
Mesh * build_trapezoid_mesh(real_t offset)
Definition ex28.cpp:38
int main()
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
Definition hooke.cpp:45
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
SparseMatrix * BuildNormalConstraints(FiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts, bool parallel)
Build a matrix constraining normal components to zero.
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
const char vishost[]