MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
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
8 // ex17 -m ../data/beam-quad.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
21 // where sigma(u)=lambda*div(u)*I+mu*(grad*u+u*grad) is the stress
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 
48 using namespace std;
49 using namespace mfem;
50 
51 // Initial displacement, used for Dirichlet boundary conditions on boundary
52 // attributes 1 and 2.
53 void InitDisplacement(const Vector &x, Vector &u);
54 
55 // A Coefficient for computing the components of the stress.
56 class StressCoefficient : public Coefficient
57 {
58 protected:
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 
65 public:
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 double Eval(ElementTransformation &T, const IntegrationPoint &ip);
73 };
74 
75 // Simple GLVis visualization manager.
76 class VisMan : public iostream
77 {
78 protected:
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 
87 public:
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.
96 void new_window (VisMan &v) { v.NewWindow(); }
97 void position_window (VisMan &v) { v.PositionWindow(); }
98 void close_connection(VisMan &v) { v.CloseConnection(); }
99 ostream &operator<<(ostream &v, void (*f)(VisMan&));
100 
101 int 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  double alpha = -1.0;
108  double kappa = -1.0;
109  bool visualization = 1;
110 
111  OptionsParser args(argc, argv);
112  args.AddOption(&mesh_file, "-m", "--mesh",
113  "Mesh file to use.");
114  args.AddOption(&ref_levels, "-r", "--refine",
115  "Number of times to refine the mesh uniformly, -1 for auto.");
116  args.AddOption(&order, "-o", "--order",
117  "Finite element order (polynomial degree).");
118  args.AddOption(&alpha, "-a", "--alpha",
119  "One of the two DG penalty parameters, typically +1/-1."
120  " See the documentation of class DGElasticityIntegrator.");
121  args.AddOption(&kappa, "-k", "--kappa",
122  "One of the two DG penalty parameters, should be positive."
123  " Negative values are replaced with (order+1)^2.");
124  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
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.
167  DG_FECollection fec(order, dim, BasisType::GaussLobatto);
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.
201  PWConstCoefficient mu_c(mu);
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);
225  a.AddDomainIntegrator(new ElasticityIntegrator(lambda_c, mu_c));
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 double 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 
333 void InitDisplacement(const Vector &x, Vector &u)
334 {
335  u = 0.0;
336  u(u.Size()-1) = -0.2*x(0);
337 }
338 
339 
341  const IntegrationPoint &ip)
342 {
343  MFEM_ASSERT(u != NULL, "displacement field is not set");
344 
345  double L = lambda.Eval(T, ip);
346  double M = mu.Eval(T, ip);
347  u->GetVectorGradient(T, grad);
348  if (si == sj)
349  {
350  double div_u = grad.Trace();
351  return L*div_u + 2*M*grad(si,si);
352  }
353  else
354  {
355  return M*(grad(si,sj) + grad(sj,si));
356  }
357 }
358 
359 
360 VisMan::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 
373 void VisMan::NewWindow()
374 {
375  sock.Append(new socketstream(host, port));
376  sid = sock.Size()-1;
377  iostream::rdbuf(sock[sid]->rdbuf());
378 }
379 
380 void VisMan::CloseConnection()
381 {
382  if (sid < sock.Size())
383  {
384  delete sock[sid];
385  sock[sid] = NULL;
386  iostream::rdbuf(0);
387  }
388 }
389 
390 void 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 
398 VisMan::~VisMan()
399 {
400  for (int i = sock.Size()-1; i >= 0; i--)
401  {
402  delete sock[i];
403  }
404 }
405 
406 
407 ostream &operator<<(ostream &v, void (*f)(VisMan&))
408 {
409  VisMan *vp = dynamic_cast<VisMan*>(&v);
410  if (vp) { (*f)(*vp); }
411  return v;
412 }
std::ostream & operator<<(std::ostream &out, const Mesh &mesh)
Definition: mesh.cpp:10127
double Eval(ElementTransformation &T, const IntegrationPoint &ip, double t)
Evaluate the coefficient in the element described by T at the point ip at time t. ...
Definition: coefficient.hpp:55
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1188
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.
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:79
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().
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:693
void PrintInfo(std::ostream &out) const
Print various sparse matrix statistics.
Definition: sparsemat.cpp:2861
int main(int argc, char *argv[])
Definition: ex1.cpp:62
Data type for Gauss-Seidel smoother of sparse matrix.
void position_window(VisMan &v)
Definition: ex17.cpp:97
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:580
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2731
Data type sparse matrix.
Definition: sparsemat.hpp:40
void new_window(VisMan &v)
Definition: ex17.cpp:96
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7982
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3924
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
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:529
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:66
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:384
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3885
int Dimension() const
Definition: mesh.hpp:744
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void close_connection(VisMan &v)
Definition: ex17.cpp:98
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:189
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:591
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
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:76
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:191
class for piecewise constant coefficient
Definition: coefficient.hpp:82
Class for integration point with weight.
Definition: intrules.hpp:25
void InitDisplacement(const Vector &x, Vector &u)
Definition: ex17.cpp:333
int dim
Definition: ex24.cpp:43
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator. Assumes ownership of bfi.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:1685
const double alpha
Definition: ex15.cpp:336
double kappa
Definition: ex3.cpp:54
Vector data type.
Definition: vector.hpp:48
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6203
Class for linear form - Vector with associated FE space and LFIntegrators.
Definition: linearform.hpp:23
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:187
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:2382
int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, int m, double &tol, double atol, int printit)
GMRES method. (tolerances are squared)
Definition: solvers.cpp:942
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator. Assumes ownership of bfi.
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:143
void Neg()
(*this) = -(*this)
Definition: vector.cpp:230
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:2285
bool Good() const
Definition: optparser.hpp:125