MFEM  v4.0 Finite element discretization library
ex19.cpp
Go to the documentation of this file.
1 // MFEM Example 19
2 //
3 // Compile with: make ex19
4 //
5 // Sample runs:
7 // ex19 -m ../data/beam-tri.mesh
8 // ex19 -m ../data/beam-hex.mesh
9 // ex19 -m ../data/beam-tet.mesh
10 // ex19 -m ../data/beam-wedge.mesh
11 //
12 // Description: This examples solves a quasi-static incompressible nonlinear
13 // elasticity problem of the form 0 = H(x), where H is an
14 // incompressible hyperelastic model and x is a block state vector
15 // containing displacement and pressure variables. The geometry of
16 // the domain is assumed to be as follows:
17 //
18 // +---------------------+
19 // boundary --->| |<--- boundary
20 // attribute 1 | | attribute 2
21 // (fixed) +---------------------+ (fixed, nonzero)
22 //
23 // The example demonstrates the use of block nonlinear operators
24 // (the class RubberOperator defining H(x)) as well as a nonlinear
25 // Newton solver for the quasi-static problem. Each Newton step
26 // requires the inversion of a Jacobian matrix, which is done
27 // through a (preconditioned) inner solver. The specialized block
28 // preconditioner is implemented as a user-defined solver.
29 //
30 // We recommend viewing examples 2, 5, and 10 before viewing this
31 // example.
32
33 #include "mfem.hpp"
34 #include <memory>
35 #include <iostream>
36 #include <fstream>
37
38 using namespace std;
39 using namespace mfem;
40
41 // Custom block preconditioner for the Jacobian of the incompressible nonlinear
42 // elasticity operator. It has the form
43 //
44 // P^-1 = [ K^-1 0 ][ I -B^T ][ I 0 ]
45 // [ 0 I ][ 0 I ][ 0 -\gamma S^-1 ]
46 //
47 // where the original Jacobian has the form
48 //
49 // J = [ K B^T ]
50 // [ B 0 ]
51 //
52 // and K^-1 is an approximation of the inverse of the displacement part of the
53 // Jacobian and S^-1 is an approximation of the inverse of the Schur
54 // complement S = B K^-1 B^T. The Schur complement is approximated using
55 // a mass matrix of the pressure variables.
56 class JacobianPreconditioner : public Solver
57 {
58 protected:
59  // Finite element spaces for setting up preconditioner blocks
61
62  // Offsets for extracting block vector segments
63  Array<int> &block_offsets;
64
65  // Jacobian for block access
66  BlockOperator *jacobian;
67
68  // Scaling factor for the pressure mass matrix in the block preconditioner
69  double gamma;
70
71  // Objects for the block preconditioner application
72  SparseMatrix *pressure_mass;
73  Solver *mass_pcg;
74  Solver *mass_prec;
75  Solver *stiff_pcg;
76  Solver *stiff_prec;
77
78 public:
79  JacobianPreconditioner(Array<FiniteElementSpace *> &fes,
80  SparseMatrix &mass, Array<int> &offsets);
81
82  virtual void Mult(const Vector &k, Vector &y) const;
83  virtual void SetOperator(const Operator &op);
84
85  virtual ~JacobianPreconditioner();
86 };
87
88 // After spatial discretization, the rubber model can be written as:
89 // 0 = H(x)
90 // where x is the block vector representing the deformation and pressure and
91 // H(x) is the nonlinear incompressible neo-Hookean operator.
92 class RubberOperator : public Operator
93 {
94 protected:
95  // Finite element spaces
97
98  // Block nonlinear form
99  BlockNonlinearForm *Hform;
100
101  // Pressure mass matrix for the preconditioner
102  SparseMatrix *pressure_mass;
103
104  // Newton solver for the hyperelastic operator
105  NewtonSolver newton_solver;
106
107  // Solver for the Jacobian solve in the Newton method
108  Solver *j_solver;
109
110  // Preconditioner for the Jacobian
111  Solver *j_prec;
112
113  // Shear modulus coefficient
114  Coefficient &mu;
115
116  // Block offsets for variable access
117  Array<int> &block_offsets;
118
119 public:
120  RubberOperator(Array<FiniteElementSpace *> &fes, Array<Array<int> *>&ess_bdr,
121  Array<int> &block_trueOffsets, double rel_tol, double abs_tol,
122  int iter, Coefficient &mu);
123
124  // Required to use the native newton solver
125  virtual Operator &GetGradient(const Vector &xp) const;
126  virtual void Mult(const Vector &k, Vector &y) const;
127
128  // Driver for the newton solver
129  void Solve(Vector &xp) const;
130
131  virtual ~RubberOperator();
132 };
133
134 // Visualization driver
135 void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
136  GridFunction *field, const char *field_name = NULL,
137  bool init_vis = false);
138
139 // Configuration definition functions
140 void ReferenceConfiguration(const Vector &x, Vector &y);
141 void InitialDeformation(const Vector &x, Vector &y);
142
143
144 int main(int argc, char *argv[])
145 {
146  // 1. Parse command-line options
147  const char *mesh_file = "../data/beam-tet.mesh";
148  int ref_levels = 0;
149  int order = 2;
150  bool visualization = true;
151  double newton_rel_tol = 1e-4;
152  double newton_abs_tol = 1e-6;
153  int newton_iter = 500;
154  double mu = 1.0;
155
156  OptionsParser args(argc, argv);
158  "Mesh file to use.");
160  "Number of times to refine the mesh uniformly.");
162  "Order (degree) of the finite elements.");
164  "--no-visualization",
165  "Enable or disable GLVis visualization.");
167  "Relative tolerance for the Newton solve.");
169  "Absolute tolerance for the Newton solve.");
171  "Maximum iterations for the Newton solve.");
173  "Shear modulus for the neo-Hookean material.");
174  args.Parse();
175  if (!args.Good())
176  {
177  args.PrintUsage(cout);
178  return 1;
179  }
180  args.PrintOptions(cout);
181
182  // 2. Read the mesh from the given mesh file. We can handle triangular,
183  // quadrilateral, tetrahedral and hexahedral meshes with the same code.
184  Mesh *mesh = new Mesh(mesh_file, 1, 1);
185  int dim = mesh->Dimension();
186
187  // 3. Refine the mesh to increase the resolution. In this example we do
188  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
189  // command-line parameter.
190  for (int lev = 0; lev < ref_levels; lev++)
191  {
192  mesh->UniformRefinement();
193  }
194
195  // 4. Define the shear modulus for the incompressible Neo-Hookean material
196  ConstantCoefficient c_mu(mu);
197
198  // 5. Define the finite element spaces for displacement and pressure
199  // (Taylor-Hood elements). By default, the displacement (u/x) is a second
200  // order vector field, while the pressure (p) is a linear scalar function.
202  H1_FECollection lin_coll(order-1, dim);
203
204  FiniteElementSpace R_space(mesh, &quad_coll, dim, Ordering::byVDIM);
205  FiniteElementSpace W_space(mesh, &lin_coll);
206
207  Array<FiniteElementSpace *> spaces(2);
208  spaces[0] = &R_space;
209  spaces[1] = &W_space;
210
211  int R_size = R_space.GetVSize();
212  int W_size = W_space.GetVSize();
213
214  // 6. Define the Dirichlet conditions (set to boundary attribute 1 and 2)
215  Array<Array<int> *> ess_bdr(2);
216
217  Array<int> ess_bdr_u(R_space.GetMesh()->bdr_attributes.Max());
218  Array<int> ess_bdr_p(W_space.GetMesh()->bdr_attributes.Max());
219
220  ess_bdr_p = 0;
221  ess_bdr_u = 0;
222  ess_bdr_u[0] = 1;
223  ess_bdr_u[1] = 1;
224
225  ess_bdr[0] = &ess_bdr_u;
226  ess_bdr[1] = &ess_bdr_p;
227
228  // 7. Print the mesh statistics
229  std::cout << "***********************************************************\n";
230  std::cout << "dim(u) = " << R_size << "\n";
231  std::cout << "dim(p) = " << W_size << "\n";
232  std::cout << "dim(u+p) = " << R_size + W_size << "\n";
233  std::cout << "***********************************************************\n";
234
235  // 8. Define the block structure of the solution vector (u then p)
236  Array<int> block_offsets(3);
237  block_offsets[0] = 0;
238  block_offsets[1] = R_space.GetVSize();
239  block_offsets[2] = W_space.GetVSize();
240  block_offsets.PartialSum();
241
242  BlockVector xp(block_offsets);
243
244  // 9. Define grid functions for the current configuration, reference
245  // configuration, final deformation, and pressure
246  GridFunction x_gf(&R_space);
247  GridFunction x_ref(&R_space);
248  GridFunction x_def(&R_space);
249  GridFunction p_gf(&W_space);
250
251  x_gf.MakeRef(&R_space, xp.GetBlock(0), 0);
252  p_gf.MakeRef(&W_space, xp.GetBlock(1), 0);
253
256
257  x_gf.ProjectCoefficient(deform);
258  x_ref.ProjectCoefficient(refconfig);
259  p_gf = 0.0;
260
261  // 10. Initialize the incompressible neo-Hookean operator
262  RubberOperator oper(spaces, ess_bdr, block_offsets,
263  newton_rel_tol, newton_abs_tol, newton_iter, c_mu);
264
265  // 11. Solve the Newton system
266  oper.Solve(xp);
267
268  // 12. Compute the final deformation
269  subtract(x_gf, x_ref, x_def);
270
271  // 13. Visualize the results if requested
272  socketstream vis_u, vis_p;
273  if (visualization)
274  {
275  char vishost[] = "localhost";
276  int visport = 19916;
277  vis_u.open(vishost, visport);
278  vis_u.precision(8);
279  visualize(vis_u, mesh, &x_gf, &x_def, "Deformation", true);
280  vis_p.open(vishost, visport);
281  vis_p.precision(8);
282  visualize(vis_p, mesh, &x_gf, &p_gf, "Pressure", true);
283  }
284
285  // 14. Save the displaced mesh, the final deformation, and the pressure
286  {
287  GridFunction *nodes = &x_gf;
288  int owns_nodes = 0;
289  mesh->SwapNodes(nodes, owns_nodes);
290
291  ofstream mesh_ofs("deformed.mesh");
292  mesh_ofs.precision(8);
293  mesh->Print(mesh_ofs);
294
295  ofstream pressure_ofs("pressure.sol");
296  pressure_ofs.precision(8);
297  p_gf.Save(pressure_ofs);
298
299  ofstream deformation_ofs("deformation.sol");
300  deformation_ofs.precision(8);
301  x_def.Save(deformation_ofs);
302  }
303
304  // 15. Free the used memory
305  delete mesh;
306
307  return 0;
308 }
309
310
311 JacobianPreconditioner::JacobianPreconditioner(Array<FiniteElementSpace *> &fes,
312  SparseMatrix &mass,
313  Array<int> &offsets)
314  : Solver(offsets[2]), block_offsets(offsets), pressure_mass(&mass)
315 {
316  fes.Copy(spaces);
317
318  gamma = 0.00001;
319
320  // The mass matrix and preconditioner do not change every Newton cycle, so we
321  // only need to define them once
322  GSSmoother *mass_prec_gs = new GSSmoother(*pressure_mass);
323
324  mass_prec = mass_prec_gs;
325
326  CGSolver *mass_pcg_iter = new CGSolver();
327  mass_pcg_iter->SetRelTol(1e-12);
328  mass_pcg_iter->SetAbsTol(1e-12);
329  mass_pcg_iter->SetMaxIter(200);
330  mass_pcg_iter->SetPrintLevel(0);
331  mass_pcg_iter->SetPreconditioner(*mass_prec);
332  mass_pcg_iter->SetOperator(*pressure_mass);
333  mass_pcg_iter->iterative_mode = false;
334
335  mass_pcg = mass_pcg_iter;
336
337  // The stiffness matrix does change every Newton cycle, so we will define it
338  // during SetOperator
339  stiff_pcg = NULL;
340  stiff_prec = NULL;
341 }
342
343 void JacobianPreconditioner::Mult(const Vector &k, Vector &y) const
344 {
345  // Extract the blocks from the input and output vectors
346  Vector disp_in(k.GetData() + block_offsets[0],
347  block_offsets[1]-block_offsets[0]);
348  Vector pres_in(k.GetData() + block_offsets[1],
349  block_offsets[2]-block_offsets[1]);
350
351  Vector disp_out(y.GetData() + block_offsets[0],
352  block_offsets[1]-block_offsets[0]);
353  Vector pres_out(y.GetData() + block_offsets[1],
354  block_offsets[2]-block_offsets[1]);
355
356  Vector temp(block_offsets[1]-block_offsets[0]);
357  Vector temp2(block_offsets[1]-block_offsets[0]);
358
359  // Perform the block elimination for the preconditioner
360  mass_pcg->Mult(pres_in, pres_out);
361  pres_out *= -gamma;
362
363  jacobian->GetBlock(0,1).Mult(pres_out, temp);
364  subtract(disp_in, temp, temp2);
365
366  stiff_pcg->Mult(temp2, disp_out);
367 }
368
369 void JacobianPreconditioner::SetOperator(const Operator &op)
370 {
371  jacobian = (BlockOperator *) &op;
372
373  // Initialize the stiffness preconditioner and solver
374  if (stiff_prec == NULL)
375  {
376  GSSmoother *stiff_prec_gs = new GSSmoother();
377
378  stiff_prec = stiff_prec_gs;
379
380  GMRESSolver *stiff_pcg_iter = new GMRESSolver();
381  stiff_pcg_iter->SetRelTol(1e-8);
382  stiff_pcg_iter->SetAbsTol(1e-8);
383  stiff_pcg_iter->SetMaxIter(200);
384  stiff_pcg_iter->SetPrintLevel(0);
385  stiff_pcg_iter->SetPreconditioner(*stiff_prec);
386  stiff_pcg_iter->iterative_mode = false;
387
388  stiff_pcg = stiff_pcg_iter;
389  }
390
391  // At each Newton cycle, compute the new stiffness preconditioner by updating
392  // the iterative solver which, in turn, updates its preconditioner
393  stiff_pcg->SetOperator(jacobian->GetBlock(0,0));
394 }
395
396 JacobianPreconditioner::~JacobianPreconditioner()
397 {
398  delete mass_pcg;
399  delete mass_prec;
400  delete stiff_prec;
401  delete stiff_pcg;
402 }
403
404
405 RubberOperator::RubberOperator(Array<FiniteElementSpace *> &fes,
406  Array<Array<int> *> &ess_bdr,
407  Array<int> &offsets,
408  double rel_tol,
409  double abs_tol,
410  int iter,
411  Coefficient &c_mu)
412  : Operator(fes[0]->GetVSize() + fes[1]->GetVSize()),
413  newton_solver(), mu(c_mu), block_offsets(offsets)
414 {
415  Array<Vector *> rhs(2);
416  rhs = NULL; // Set all entries in the array
417
418  fes.Copy(spaces);
419
420  // Define the block nonlinear form
421  Hform = new BlockNonlinearForm(spaces);
422
423  // Add the incompressible neo-Hookean integrator
425
426  // Set the essential boundary conditions
427  Hform->SetEssentialBC(ess_bdr, rhs);
428
429  // Compute the pressure mass stiffness matrix
430  BilinearForm *a = new BilinearForm(spaces[1]);
431  ConstantCoefficient one(1.0);
433  a->Assemble();
434  a->Finalize();
435  pressure_mass = a->LoseMat();
436  delete a;
437
438  // Initialize the Jacobian preconditioner
439  JacobianPreconditioner *jac_prec =
440  new JacobianPreconditioner(fes, *pressure_mass, block_offsets);
441  j_prec = jac_prec;
442
443  // Set up the Jacobian solver
444  GMRESSolver *j_gmres = new GMRESSolver();
445  j_gmres->iterative_mode = false;
446  j_gmres->SetRelTol(1e-12);
447  j_gmres->SetAbsTol(1e-12);
448  j_gmres->SetMaxIter(300);
449  j_gmres->SetPrintLevel(0);
450  j_gmres->SetPreconditioner(*j_prec);
451  j_solver = j_gmres;
452
453  // Set the newton solve parameters
454  newton_solver.iterative_mode = true;
455  newton_solver.SetSolver(*j_solver);
456  newton_solver.SetOperator(*this);
457  newton_solver.SetPrintLevel(1);
458  newton_solver.SetRelTol(rel_tol);
459  newton_solver.SetAbsTol(abs_tol);
460  newton_solver.SetMaxIter(iter);
461 }
462
463 // Solve the Newton system
464 void RubberOperator::Solve(Vector &xp) const
465 {
466  Vector zero;
467  newton_solver.Mult(zero, xp);
468  MFEM_VERIFY(newton_solver.GetConverged(),
469  "Newton Solver did not converge.");
470 }
471
472 // compute: y = H(x,p)
473 void RubberOperator::Mult(const Vector &k, Vector &y) const
474 {
475  Hform->Mult(k, y);
476 }
477
478 // Compute the Jacobian from the nonlinear form
479 Operator &RubberOperator::GetGradient(const Vector &xp) const
480 {
482 }
483
484 RubberOperator::~RubberOperator()
485 {
486  delete Hform;
487  delete pressure_mass;
488  delete j_solver;
489  delete j_prec;
490 }
491
492
493 // Inline visualization
494 void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
495  GridFunction *field, const char *field_name, bool init_vis)
496 {
497  if (!out)
498  {
499  return;
500  }
501
502  GridFunction *nodes = deformed_nodes;
503  int owns_nodes = 0;
504
505  mesh->SwapNodes(nodes, owns_nodes);
506
507  out << "solution\n" << *mesh << *field;
508
509  mesh->SwapNodes(nodes, owns_nodes);
510
511  if (init_vis)
512  {
513  out << "window_size 800 800\n";
514  out << "window_title '" << field_name << "'\n";
515  if (mesh->SpaceDimension() == 2)
516  {
517  out << "view 0 0\n"; // view from top
518  out << "keys jlA\n"; // turn off perspective and light, +anti-aliasing
519  }
520  out << "keys cmA\n"; // show colorbar and mesh, +anti-aliasing
521  out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
522  }
523  out << flush;
524 }
525
527 {
528  // Set the reference, stress free, configuration
529  y = x;
530 }
531
532 void InitialDeformation(const Vector &x, Vector &y)
533 {
534  // Set the initial configuration. Having this different from the reference
535  // configuration can help convergence
536  y = x;
537  y[1] = x[1] + 0.25*x[0];
538 }
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
Definition: ex10.cpp:377
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:347
void InitialDeformation(const Vector &x, Vector &y)
Definition: ex10.cpp:589
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1156
Definition: solvers.hpp:111
void ReferenceConfiguration(const Vector &x, Vector &y)
Definition: ex19.cpp:526
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
Subclass constant coefficient.
Definition: coefficient.hpp:67
virtual Operator & GetGradient(const Vector &x) const
Evaluate the gradient operator at the point x. The default behavior in class Operator is to generate ...
Definition: operator.hpp:68
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Definition: mesh.cpp:5875
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
void Copy(Array &copy) const
Create a copy of the current array.
Definition: array.hpp:795
A class representing a general block nonlinear operator defined on the Cartesian product of multiple ...
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:283
int main(int argc, char *argv[])
Definition: ex1.cpp:58
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:159
Data type for Gauss-Seidel smoother of sparse matrix.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2627
int dim
Definition: ex3.cpp:48
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:67
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:94
Data type sparse matrix.
Definition: sparsemat.hpp:40
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:272
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7591
void SetMaxIter(int max_it)
Definition: solvers.hpp:63
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:259
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition: gridfunc.cpp:188
int Dimension() const
Definition: mesh.hpp:713
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
int SpaceDimension() const
Definition: mesh.hpp:714
void SetAbsTol(double atol)
Definition: solvers.hpp:62
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:179
void SetRelTol(double rtol)
Definition: solvers.hpp:61
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:85
SparseMatrix * LoseMat()
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:385
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
void PartialSum()
Partial Sum.
Definition: array.cpp:103
GMRES method.
Definition: solvers.hpp:143
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
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:1582
int open(const char hostname[], int port)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:125
Vector data type.
Definition: vector.hpp:48
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:88
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
Base class for solvers.
Definition: operator.hpp:279
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:64
Abstract operator.
Definition: operator.hpp:21
A class to handle Block systems in a matrix-free implementation.
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:84
bool Good() const
Definition: optparser.hpp:122