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