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