MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex17.cpp
Go to the documentation of this file.
1// MFEM Example 17
2//
3// Compile with: make ex17
4//
5// Sample runs:
6//
7// ex17 -m ../data/beam-tri.mesh
9// ex17 -m ../data/beam-tet.mesh
10// ex17 -m ../data/beam-hex.mesh
11// ex17 -m ../data/beam-wedge.mesh
12// ex17 -m ../data/beam-quad.mesh -r 2 -o 3
13// ex17 -m ../data/beam-quad.mesh -r 2 -o 2 -a 1 -k 1
14// ex17 -m ../data/beam-hex.mesh -r 2 -o 2
15//
16// Description: This example code solves a simple linear elasticity problem
17// describing a multi-material cantilever beam using symmetric or
18// non-symmetric discontinuous Galerkin (DG) formulation.
19//
20// Specifically, we approximate the weak form of -div(sigma(u))=0
22// tensor corresponding to displacement field u, and lambda and mu
23// are the material Lame constants. The boundary conditions are
24// Dirichlet, u=u_D on the fixed part of the boundary, namely
25// boundary attributes 1 and 2; on the rest of the boundary we use
26// sigma(u).n=0 b.c. The geometry of the domain is assumed to be
27// as follows:
28//
29// +----------+----------+
30// boundary --->| material | material |<--- boundary
31// attribute 1 | 1 | 2 | attribute 2
32// (fixed) +----------+----------+ (fixed, nonzero)
33//
34// The example demonstrates the use of high-order DG vector finite
35// element spaces with the linear DG elasticity bilinear form,
36// meshes with curved elements, and the definition of piece-wise
37// constant and function vector-coefficient objects. The use of
38// non-homogeneous Dirichlet b.c. imposed weakly, is also
39// illustrated.
40//
41// We recommend viewing examples 2 and 14 before viewing this
42// example.
43
44#include "mfem.hpp"
45#include <fstream>
46#include <iostream>
47
48using namespace std;
49using namespace mfem;
50
51// Initial displacement, used for Dirichlet boundary conditions on boundary
52// attributes 1 and 2.
53void InitDisplacement(const Vector &x, Vector &u);
54
55// A Coefficient for computing the components of the stress.
56class StressCoefficient : public Coefficient
57{
58protected:
59 Coefficient &lambda, &mu;
60 GridFunction *u; // displacement
61 int si, sj; // component of the stress to evaluate, 0 <= si,sj < dim
62
63 DenseMatrix grad; // auxiliary matrix, used in Eval
64
65public:
66 StressCoefficient(Coefficient &lambda_, Coefficient &mu_)
67 : lambda(lambda_), mu(mu_), u(NULL), si(0), sj(0) { }
68
69 void SetDisplacement(GridFunction &u_) { u = &u_; }
70 void SetComponent(int i, int j) { si = i; sj = j; }
71
72 virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip);
73};
74
75// Simple GLVis visualization manager.
76class VisMan : public iostream
77{
78protected:
79 const char *host;
80 int port;
82 int sid; // active socket, index inside 'sock'.
83
84 int win_x, win_y, win_w, win_h;
85 int win_stride_x, win_stride_y, win_nx;
86
87public:
88 VisMan(const char *vishost, const int visport);
89 void NewWindow();
90 void CloseConnection();
91 void PositionWindow();
92 virtual ~VisMan();
93};
94
95// Manipulators for the GLVis visualization manager.
96void new_window (VisMan &v) { v.NewWindow(); }
97void position_window (VisMan &v) { v.PositionWindow(); }
98void close_connection(VisMan &v) { v.CloseConnection(); }
99ostream &operator<<(ostream &v, void (*f)(VisMan&));
100
101int main(int argc, char *argv[])
102{
103 // 1. Define and parse command-line options.
104 const char *mesh_file = "../data/beam-tri.mesh";
105 int ref_levels = -1;
106 int order = 1;
107 real_t alpha = -1.0;
108 real_t kappa = -1.0;
109 bool visualization = 1;
110
111 OptionsParser args(argc, argv);
113 "Mesh file to use.");
115 "Number of times to refine the mesh uniformly, -1 for auto.");
117 "Finite element order (polynomial degree).");
119 "One of the two DG penalty parameters, typically +1/-1."
120 " See the documentation of class DGElasticityIntegrator.");
122 "One of the two DG penalty parameters, should be positive."
123 " Negative values are replaced with (order+1)^2.");
125 "--no-visualization",
126 "Enable or disable GLVis visualization.");
127 args.Parse();
128 if (!args.Good())
129 {
130 args.PrintUsage(cout);
131 return 1;
132 }
133 if (kappa < 0)
134 {
135 kappa = (order+1)*(order+1);
136 }
137 args.PrintOptions(cout);
138
139 // 2. Read the mesh from the given mesh file.
140 Mesh mesh(mesh_file, 1, 1);
141 int dim = mesh.Dimension();
142
143 if (mesh.attributes.Max() < 2 || mesh.bdr_attributes.Max() < 2)
144 {
145 cerr << "\nInput mesh should have at least two materials and "
146 << "two boundary attributes! (See schematic in ex17.cpp)\n"
147 << endl;
148 return 3;
149 }
150
151 // 3. Refine the mesh to increase the resolution.
152 if (ref_levels < 0)
153 {
154 ref_levels = (int)floor(log(5000./mesh.GetNE())/log(2.)/dim);
155 }
156 for (int l = 0; l < ref_levels; l++)
157 {
158 mesh.UniformRefinement();
159 }
160 // Since NURBS meshes do not support DG integrators, we convert them to
161 // regular polynomial mesh of the specified (solution) order.
162 if (mesh.NURBSext) { mesh.SetCurvature(order); }
163
164 // 4. Define a DG vector finite element space on the mesh. Here, we use
165 // Gauss-Lobatto nodal basis because it gives rise to a sparser matrix
166 // compared to the default Gauss-Legendre nodal basis.
168 FiniteElementSpace fespace(&mesh, &fec, dim);
169
170 cout << "Number of finite element unknowns: " << fespace.GetTrueVSize()
171 << "\nAssembling: " << flush;
172
173 // 5. In this example, the Dirichlet boundary conditions are defined by
174 // marking boundary attributes 1 and 2 in the marker Array 'dir_bdr'.
175 // These b.c. are imposed weakly, by adding the appropriate boundary
176 // integrators over the marked 'dir_bdr' to the bilinear and linear forms.
177 // With this DG formulation, there are no essential boundary conditions.
178 Array<int> ess_tdof_list; // no essential b.c. (empty list)
179 Array<int> dir_bdr(mesh.bdr_attributes.Max());
180 dir_bdr = 0;
181 dir_bdr[0] = 1; // boundary attribute 1 is Dirichlet
182 dir_bdr[1] = 1; // boundary attribute 2 is Dirichlet
183
184 // 6. Define the DG solution vector 'x' as a finite element grid function
185 // corresponding to fespace. Initialize 'x' using the 'InitDisplacement'
186 // function.
187 GridFunction x(&fespace);
189 x.ProjectCoefficient(init_x);
190
191 // 7. Set up the Lame constants for the two materials. They are defined as
192 // piece-wise (with respect to the element attributes) constant
193 // coefficients, i.e. type PWConstCoefficient.
194 Vector lambda(mesh.attributes.Max());
195 lambda = 1.0; // Set lambda = 1 for all element attributes.
196 lambda(0) = 50.0; // Set lambda = 50 for element attribute 1.
197 PWConstCoefficient lambda_c(lambda);
198 Vector mu(mesh.attributes.Max());
199 mu = 1.0; // Set mu = 1 for all element attributes.
200 mu(0) = 50.0; // Set mu = 50 for element attribute 1.
202
203 // 8. Set up the linear form b(.) which corresponds to the right-hand side of
204 // the FEM linear system. In this example, the linear form b(.) consists
205 // only of the terms responsible for imposing weakly the Dirichlet
206 // boundary conditions, over the attributes marked in 'dir_bdr'. The
207 // values for the Dirichlet boundary condition are taken from the
208 // VectorFunctionCoefficient 'x_init' which in turn is based on the
209 // function 'InitDisplacement'.
210 LinearForm b(&fespace);
211 cout << "r.h.s. ... " << flush;
214 init_x, lambda_c, mu_c, alpha, kappa), dir_bdr);
215 b.Assemble();
216
217 // 9. Set up the bilinear form a(.,.) on the DG finite element space
218 // corresponding to the linear elasticity integrator with coefficients
219 // lambda and mu as defined above. The additional interior face integrator
220 // ensures the weak continuity of the displacement field. The additional
221 // boundary face integrator works together with the boundary integrator
222 // added to the linear form b(.) to impose weakly the Dirichlet boundary
223 // conditions.
224 BilinearForm a(&fespace);
227 new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa));
229 new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa), dir_bdr);
230
231 // 10. Assemble the bilinear form and the corresponding linear system.
232 cout << "matrix ... " << flush;
233 a.Assemble();
234
235 SparseMatrix A;
236 Vector B, X;
237 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
238 cout << "done." << endl;
239
240 // Print some information about the matrix of the linear system.
241 A.PrintInfo(cout);
242
243#ifndef MFEM_USE_SUITESPARSE
244 // 11. Define a simple symmetric Gauss-Seidel preconditioner and use it to
245 // solve the system Ax=b with PCG for the symmetric formulation, or GMRES
246 // for the non-symmetric.
247 GSSmoother M(A);
248 const real_t rtol = 1e-6;
249 if (alpha == -1.0)
250 {
251 PCG(A, M, B, X, 3, 5000, rtol*rtol, 0.0);
252 }
253 else
254 {
255 GMRES(A, M, B, X, 3, 5000, 100, rtol*rtol, 0.0);
256 }
257#else
258 // 11. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
259 UMFPackSolver umf_solver;
260 umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
261 umf_solver.SetOperator(A);
262 umf_solver.Mult(B, X);
263#endif
264
265 // 12. Recover the solution as a finite element grid function 'x'.
266 a.RecoverFEMSolution(X, b, x);
267
268 // 13. Use the DG solution space as the mesh nodal space. This allows us to
269 // save the displaced mesh as a curved DG mesh.
270 mesh.SetNodalFESpace(&fespace);
271
272 Vector reference_nodes;
273 if (visualization) { reference_nodes = *mesh.GetNodes(); }
274
275 // 14. Save the displaced mesh and minus the solution (which gives the
276 // backward displacements to the reference mesh). This output can be
277 // viewed later using GLVis: "glvis -m displaced.mesh -g sol.gf".
278 {
279 *mesh.GetNodes() += x;
280 x.Neg(); // x = -x
281 ofstream mesh_ofs("displaced.mesh");
282 mesh_ofs.precision(8);
283 mesh.Print(mesh_ofs);
284 ofstream sol_ofs("sol.gf");
285 sol_ofs.precision(8);
286 x.Save(sol_ofs);
287 }
288
289 // 15. Visualization: send data by socket to a GLVis server.
290 if (visualization)
291 {
292 char vishost[] = "localhost";
293 int visport = 19916;
294 VisMan vis(vishost, visport);
295 const char *glvis_keys = (dim < 3) ? "Rjlc" : "c";
296
297 // Visualize the deformed configuration.
298 vis << new_window << setprecision(8)
299 << "solution\n" << mesh << x << flush
300 << "keys " << glvis_keys << endl
301 << "window_title 'Deformed configuration'" << endl
302 << "plot_caption 'Backward displacement'" << endl
304
305 // Visualize the stress components.
306 const char *c = "xyz";
307 FiniteElementSpace scalar_dg_space(&mesh, &fec);
308 GridFunction stress(&scalar_dg_space);
309 StressCoefficient stress_c(lambda_c, mu_c);
310 *mesh.GetNodes() = reference_nodes;
311 x.Neg(); // x = -x
312 stress_c.SetDisplacement(x);
313 for (int si = 0; si < dim; si++)
314 {
315 for (int sj = si; sj < dim; sj++)
316 {
317 stress_c.SetComponent(si, sj);
318 stress.ProjectCoefficient(stress_c);
319
320 vis << new_window << setprecision(8)
321 << "solution\n" << mesh << stress << flush
322 << "keys " << glvis_keys << endl
323 << "window_title |Stress " << c[si] << c[sj] << '|' << endl
325 }
326 }
327 }
328
329 return 0;
330}
331
332
334{
335 u = 0.0;
336 u(u.Size()-1) = -0.2*x(0);
337}
338
339
340real_t StressCoefficient::Eval(ElementTransformation &T,
341 const IntegrationPoint &ip)
342{
343 MFEM_ASSERT(u != NULL, "displacement field is not set");
344
345 real_t L = lambda.Eval(T, ip);
346 real_t M = mu.Eval(T, ip);
348 if (si == sj)
349 {
352 }
353 else
354 {
356 }
357}
358
359
360VisMan::VisMan(const char *vishost, const int visport)
361 : iostream(0),
362 host(vishost), port(visport), sid(0)
363{
364 win_x = 0;
365 win_y = 0;
366 win_w = 400; // window width
367 win_h = 350; // window height
368 win_stride_x = win_w;
369 win_stride_y = win_h + 20;
370 win_nx = 4; // number of windows in a row
371}
372
373void VisMan::NewWindow()
374{
375 sock.Append(new socketstream(host, port));
376 sid = sock.Size()-1;
377 iostream::rdbuf(sock[sid]->rdbuf());
378}
379
380void VisMan::CloseConnection()
381{
382 if (sid < sock.Size())
383 {
384 delete sock[sid];
385 sock[sid] = NULL;
386 iostream::rdbuf(0);
387 }
388}
389
390void VisMan::PositionWindow()
391{
392 *this << "window_geometry "
393 << win_x + win_stride_x*(sid%win_nx) << ' '
394 << win_y + win_stride_y*(sid/win_nx) << ' '
395 << win_w << ' ' << win_h << endl;
396}
397
398VisMan::~VisMan()
399{
400 for (int i = sock.Size()-1; i >= 0; i--)
401 {
402 delete sock[i];
403 }
404}
405
406ostream &operator<<(ostream &v, void (*f)(VisMan&))
407{
408 VisMan *vp = dynamic_cast<VisMan*>(&v);
409 if (vp) { (*f)(*vp); }
410 return v;
411}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
@ GaussLobatto
Closed type.
Definition fe_base.hpp:32
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
real_t Trace() const
Trace of a square matrix.
Definition densemat.cpp:511
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
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
Compute the vector gradient with respect to the physical element variable.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Class for integration point with weight.
Definition intrules.hpp:35
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
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
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:290
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
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
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition mesh.cpp:6153
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition mesh.cpp:6211
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
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.
Data type sparse matrix.
Definition sparsemat.hpp:51
void PrintInfo(std::ostream &out) const
Print various sparse matrix statistics.
Direct sparse solver using UMFPACK.
Definition solvers.hpp:1096
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition solvers.cpp:3102
real_t Control[UMFPACK_CONTROL]
Definition solvers.hpp:1106
virtual void Mult(const Vector &b, Vector &x) const
Direct solution of the linear system using UMFPACK.
Definition solvers.cpp:3197
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
void Neg()
(*this) = -(*this)
Definition vector.cpp:300
const real_t alpha
Definition ex15.cpp:369
void close_connection(VisMan &v)
Definition ex17.cpp:98
void InitDisplacement(const Vector &x, Vector &u)
Definition ex17.cpp:333
void position_window(VisMan &v)
Definition ex17.cpp:97
void new_window(VisMan &v)
Definition ex17.cpp:96
real_t kappa
Definition ex24.cpp:54
int dim
Definition ex24.cpp:53
real_t mu
Definition ex25.cpp:140
int main()
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, int m, real_t &tol, real_t atol, int printit)
GMRES method. (tolerances are squared)
Definition solvers.cpp:1342
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition solvers.cpp:913
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[]