MFEM  v3.4 Finite element discretization library
ex10.cpp
Go to the documentation of this file.
1 // MFEM Example 10
2 // SUNDIALS Modification
3 //
4 // Compile with: make ex10
5 //
6 // Sample runs:
7 // ex10 -m ../../data/beam-quad.mesh -r 2 -o 2 -s 5 -dt 0.15 -vs 10
8 // ex10 -m ../../data/beam-tri.mesh -r 2 -o 2 -s 7 -dt 0.3 -vs 5
9 // ex10 -m ../../data/beam-hex.mesh -r 1 -o 2 -s 5 -dt 0.2 -vs 5
10 // ex10 -m ../../data/beam-tri.mesh -r 2 -o 2 -s 2 -dt 3 -nls kinsol
11 // ex10 -m ../../data/beam-quad.mesh -r 2 -o 2 -s 2 -dt 3 -nls kinsol
12 // ex10 -m ../../data/beam-hex.mesh -r 1 -o 2 -s 2 -dt 3 -nls kinsol
13 // ex10 -m ../../data/beam-quad.mesh -r 2 -o 2 -s 15 -dt 5e-3 -vs 60
14 // ex10 -m ../../data/beam-tri.mesh -r 2 -o 2 -s 16 -dt 0.01 -vs 30
15 // ex10 -m ../../data/beam-hex.mesh -r 1 -o 2 -s 15 -dt 0.01 -vs 30
16 // ex10 -m ../../data/beam-quad-amr.mesh -r 2 -o 2 -s 5 -dt 0.15 -vs 10
17 //
18 // Description: This examples solves a time dependent nonlinear elasticity
19 // problem of the form dv/dt = H(x) + S v, dx/dt = v, where H is a
20 // hyperelastic model and S is a viscosity operator of Laplacian
21 // type. The geometry of the domain is assumed to be as follows:
22 //
23 // +---------------------+
24 // boundary --->| |
25 // attribute 1 | |
26 // (fixed) +---------------------+
27 //
28 // The example demonstrates the use of nonlinear operators (the
29 // class HyperelasticOperator defining H(x)), as well as their
30 // implicit time integration using a Newton method for solving an
31 // associated reduced backward-Euler type nonlinear equation
32 // (class ReducedSystemOperator). Each Newton step requires the
33 // inversion of a Jacobian matrix, which is done through a
34 // (preconditioned) inner solver. Note that implementing the
35 // method HyperelasticOperator::ImplicitSolve is the only
36 // requirement for high-order implicit (SDIRK) time integration.
37 //
38 // We recommend viewing examples 2 and 9 before viewing this
39 // example.
40
41 #include "mfem.hpp"
42 #include <memory>
43 #include <iostream>
44 #include <fstream>
45 #include <string>
46 #include <map>
47
48 #ifndef MFEM_USE_SUNDIALS
49 #error This example requires that MFEM is built with MFEM_USE_SUNDIALS=YES
50 #endif
51
52 using namespace std;
53 using namespace mfem;
54
55 class ReducedSystemOperator;
56 class SundialsJacSolver;
57
58 /** After spatial discretization, the hyperelastic model can be written as a
59  * system of ODEs:
60  * dv/dt = -M^{-1}*(H(x) + S*v)
61  * dx/dt = v,
62  * where x is the vector representing the deformation, v is the velocity field,
63  * M is the mass matrix, S is the viscosity matrix, and H(x) is the nonlinear
64  * hyperelastic operator.
65  *
66  * Class HyperelasticOperator represents the right-hand side of the above
67  * system of ODEs. */
68 class HyperelasticOperator : public TimeDependentOperator
69 {
70 protected:
71  FiniteElementSpace &fespace;
72
73  BilinearForm M, S;
74  NonlinearForm H;
75  double viscosity;
76  HyperelasticModel *model;
77
78  CGSolver M_solver; // Krylov solver for inverting the mass matrix M
79  DSmoother M_prec; // Preconditioner for the mass matrix M
80
81  /** Nonlinear operator defining the reduced backward Euler equation for the
82  velocity. Used in the implementation of method ImplicitSolve. */
83  ReducedSystemOperator *reduced_oper;
84
85  /// Newton solver for the reduced backward Euler equation
86  NewtonSolver *newton_solver;
87
88  /// Solver for the Jacobian solve in the Newton method
89  Solver *J_solver;
90  /// Preconditioner for the Jacobian solve in the Newton method
91  Solver *J_prec;
92
93  mutable Vector z; // auxiliary vector
94
95 public:
96  /// Solver type to use in the ImplicitSolve() method, used by SDIRK methods.
97  enum NonlinearSolverType
98  {
99  NEWTON = 0, ///< Use MFEM's plain NewtonSolver
100  KINSOL = 1 ///< Use SUNDIALS' KINSOL (through MFEM's class KinSolver)
101  };
102
103  HyperelasticOperator(FiniteElementSpace &f, Array<int> &ess_bdr,
104  double visc, double mu, double K,
105  NonlinearSolverType nls_type);
106
107  /// Compute the right-hand side of the ODE system.
108  virtual void Mult(const Vector &vx, Vector &dvx_dt) const;
109  /** Solve the Backward-Euler equation: k = f(x + dt*k, t), for the unknown k.
110  This is the only requirement for high-order SDIRK implicit integration.*/
111  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
112
113  /** Connect the Jacobian linear system solver (SundialsJacSolver) used by
114  SUNDIALS' CVODE and ARKODE time integrators to the internal objects
115  created by HyperelasticOperator. This method is called by the InitSystem
116  method of SundialsJacSolver. */
117  void InitSundialsJacSolver(SundialsJacSolver &sjsolv);
118
119  double ElasticEnergy(const Vector &x) const;
120  double KineticEnergy(const Vector &v) const;
121  void GetElasticEnergyDensity(const GridFunction &x, GridFunction &w) const;
122
123  virtual ~HyperelasticOperator();
124 };
125
126 /** Nonlinear operator of the form:
127  k --> (M + dt*S)*k + H(x + dt*v + dt^2*k) + S*v,
128  where M and S are given BilinearForms, H is a given NonlinearForm, v and x
129  are given vectors, and dt is a scalar. */
130 class ReducedSystemOperator : public Operator
131 {
132 private:
133  BilinearForm *M, *S;
134  NonlinearForm *H;
135  mutable SparseMatrix *Jacobian;
136  double dt;
137  const Vector *v, *x;
138  mutable Vector w, z;
139
140 public:
141  ReducedSystemOperator(BilinearForm *M_, BilinearForm *S_, NonlinearForm *H_);
142
143  /// Set current dt, v, x values - needed to compute action and Jacobian.
144  void SetParameters(double dt_, const Vector *v_, const Vector *x_);
145
146  /// Compute y = H(x + dt (v + dt k)) + M k + S (v + dt k).
147  virtual void Mult(const Vector &k, Vector &y) const;
148
149  /// Compute J = M + dt S + dt^2 grad_H(x + dt (v + dt k)).
150  virtual Operator &GetGradient(const Vector &k) const;
151
152  virtual ~ReducedSystemOperator();
153 };
154
155 /// Custom Jacobian system solver for the SUNDIALS time integrators.
156 /** For the ODE system represented by HyperelasticOperator
157
158  M dv/dt = -(H(x) + S*v)
159  dx/dt = v,
160
161  this class facilitates the solution of linear systems of the form
162
163  (M + γS) yv + γJ yx = M bv, J=(dH/dx)(x)
164  - γ yv + yx = bx
165
166  for given bv, bx, x, and γ = GetTimeStep(). */
167 class SundialsJacSolver : public SundialsODELinearSolver
168 {
169 private:
170  BilinearForm *M, *S;
171  NonlinearForm *H;
173  Solver *J_solver;
174
175 public:
176  SundialsJacSolver()
177  : M(), S(), H(), grad_H(), Jacobian(), J_solver() { }
178
179  /// Connect the solver to the objects created inside HyperelasticOperator.
180  void SetOperators(BilinearForm &M_, BilinearForm &S_,
181  NonlinearForm &H_, Solver &solver)
182  {
183  M = &M_; S = &S_; H = &H_; J_solver = &solver;
184  }
185
186  /** Linear solve applicable to the SUNDIALS format.
187  Solves (Mass - dt J) y = Mass b, where in our case:
188  Mass = | M 0 | J = | -S -grad_H | y = | v_hat | b = | b_v |
189  | 0 I | | I 0 | | x_hat | | b_x |
190  The result replaces the rhs b.
191  We substitute x_hat = b_x + dt v_hat and solve
192  (M + dt S + dt^2 grad_H) v_hat = M b_v - dt grad_H b_x. */
193  int InitSystem(void *sundials_mem);
194  int SetupSystem(void *sundials_mem, int conv_fail,
195  const Vector &y_pred, const Vector &f_pred, int &jac_cur,
196  Vector &v_temp1, Vector &v_temp2, Vector &v_temp3);
197  int SolveSystem(void *sundials_mem, Vector &b, const Vector &weight,
198  const Vector &y_cur, const Vector &f_cur);
199  int FreeSystem(void *sundials_mem);
200 };
201
202
203 /** Function representing the elastic energy density for the given hyperelastic
204  model+deformation. Used in HyperelasticOperator::GetElasticEnergyDensity. */
205 class ElasticEnergyCoefficient : public Coefficient
206 {
207 private:
208  HyperelasticModel &model;
209  const GridFunction &x;
210  DenseMatrix J;
211
212 public:
213  ElasticEnergyCoefficient(HyperelasticModel &m, const GridFunction &x_)
214  : model(m), x(x_) { }
215  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
216  virtual ~ElasticEnergyCoefficient() { }
217 };
218
219 void InitialDeformation(const Vector &x, Vector &y);
220
221 void InitialVelocity(const Vector &x, Vector &v);
222
223 void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
224  GridFunction *field, const char *field_name = NULL,
225  bool init_vis = false);
226
227
228 int main(int argc, char *argv[])
229 {
230  // 1. Parse command-line options.
231  const char *mesh_file = "../../data/beam-quad.mesh";
232  int ref_levels = 2;
233  int order = 2;
234  int ode_solver_type = 3;
235  double t_final = 300.0;
236  double dt = 3.0;
237  double visc = 1e-2;
238  double mu = 0.25;
239  double K = 5.0;
240  bool visualization = true;
241  const char *nls = "newton";
242  int vis_steps = 1;
243
244  // Relative and absolute tolerances for CVODE and ARKODE.
245  const double reltol = 1e-1, abstol = 1e-1;
246
247  OptionsParser args(argc, argv);
249  "Mesh file to use.");
251  "Number of times to refine the mesh uniformly.");
253  "Order (degree) of the finite elements.");
255  "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
256  " 4 - CVODE implicit, approximate Jacobian,\n\t"
257  " 5 - CVODE implicit, specified Jacobian,\n\t"
258  " 6 - ARKODE implicit, approximate Jacobian,\n\t"
259  " 7 - ARKODE implicit, specified Jacobian,\n\t"
260  " 11 - Forward Euler, 12 - RK2,\n\t"
261  " 13 - RK3 SSP, 14 - RK4,\n\t"
262  " 15 - CVODE (adaptive order) explicit,\n\t"
263  " 16 - ARKODE default (4th order) explicit.");
265  "Nonlinear systems solver: "
266  "\"newton\" (plain Newton) or \"kinsol\" (KINSOL).");
268  "Final time; start time is 0.");
270  "Time step.");
272  "Viscosity coefficient.");
274  "Shear modulus in the Neo-Hookean hyperelastic model.");
276  "Bulk modulus in the Neo-Hookean hyperelastic model.");
278  "--no-visualization",
279  "Enable or disable GLVis visualization.");
281  "Visualize every n-th timestep.");
282  args.Parse();
283  if (!args.Good())
284  {
285  args.PrintUsage(cout);
286  return 1;
287  }
288  args.PrintOptions(cout);
289
290  // 2. Read the mesh from the given mesh file. We can handle triangular,
291  // quadrilateral, tetrahedral and hexahedral meshes with the same code.
292  Mesh *mesh = new Mesh(mesh_file, 1, 1);
293  int dim = mesh->Dimension();
294
295  // 3. Define the ODE solver used for time integration. Several implicit
296  // singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as
297  // explicit Runge-Kutta methods are available.
298  ODESolver *ode_solver;
299  CVODESolver *cvode = NULL;
300  ARKODESolver *arkode = NULL;
301  SundialsJacSolver *sjsolver = NULL;
302  switch (ode_solver_type)
303  {
304  // Implicit L-stable methods
305  case 1: ode_solver = new BackwardEulerSolver; break;
306  case 2: ode_solver = new SDIRK23Solver(2); break;
307  case 3: ode_solver = new SDIRK33Solver; break;
308  case 4:
309  case 5:
310  cvode = new CVODESolver(CV_BDF, CV_NEWTON);
311  cvode->SetSStolerances(reltol, abstol);
312  cvode->SetMaxStep(dt);
313  if (ode_solver_type == 5)
314  {
315  sjsolver = new SundialsJacSolver;
316  cvode->SetLinearSolver(*sjsolver);
317  }
318  ode_solver = cvode; break;
319  case 6:
320  case 7:
321  arkode = new ARKODESolver(ARKODESolver::IMPLICIT);
322  arkode->SetSStolerances(reltol, abstol);
323  arkode->SetMaxStep(dt);
324  if (ode_solver_type == 7)
325  {
326  // Custom Jacobian inversion.
327  sjsolver = new SundialsJacSolver;
328  arkode->SetLinearSolver(*sjsolver);
329  }
330  ode_solver = arkode; break;
331  // Explicit methods
332  case 11: ode_solver = new ForwardEulerSolver; break;
333  case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
334  case 13: ode_solver = new RK3SSPSolver; break;
335  case 14: ode_solver = new RK4Solver; break;
336  case 15:
337  cvode = new CVODESolver(CV_ADAMS, CV_FUNCTIONAL);
338  cvode->SetSStolerances(reltol, abstol);
339  cvode->SetMaxStep(dt);
340  ode_solver = cvode; break;
341  case 16:
342  arkode = new ARKODESolver(ARKODESolver::IMPLICIT);
343  arkode->SetSStolerances(reltol, abstol);
344  arkode->SetMaxStep(dt);
345  ode_solver = arkode; break;
346  // Implicit A-stable methods (not L-stable)
347  case 22: ode_solver = new ImplicitMidpointSolver; break;
348  case 23: ode_solver = new SDIRK23Solver; break;
349  case 24: ode_solver = new SDIRK34Solver; break;
350  default:
351  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
352  delete mesh;
353  return 3;
354  }
355
356  map<string,HyperelasticOperator::NonlinearSolverType> nls_map;
357  nls_map["newton"] = HyperelasticOperator::NEWTON;
358  nls_map["kinsol"] = HyperelasticOperator::KINSOL;
359  if (nls_map.find(nls) == nls_map.end())
360  {
361  cout << "Unknown type of nonlinear solver: " << nls << endl;
362  return 4;
363  }
364
365  // 4. Refine the mesh to increase the resolution. In this example we do
366  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
367  // command-line parameter.
368  for (int lev = 0; lev < ref_levels; lev++)
369  {
370  mesh->UniformRefinement();
371  }
372
373  // 5. Define the vector finite element spaces representing the mesh
374  // deformation x, the velocity v, and the initial configuration, x_ref.
375  // Define also the elastic energy density, w, which is in a discontinuous
376  // higher-order space. Since x and v are integrated in time as a system,
377  // we group them together in block vector vx, with offsets given by the
378  // fe_offset array.
379  H1_FECollection fe_coll(order, dim);
380  FiniteElementSpace fespace(mesh, &fe_coll, dim);
381
382  int fe_size = fespace.GetTrueVSize();
383  cout << "Number of velocity/deformation unknowns: " << fe_size << endl;
384  Array<int> fe_offset(3);
385  fe_offset[0] = 0;
386  fe_offset[1] = fe_size;
387  fe_offset[2] = 2*fe_size;
388
389  BlockVector vx(fe_offset);
390  GridFunction v, x;
391  v.MakeTRef(&fespace, vx.GetBlock(0), 0);
392  x.MakeTRef(&fespace, vx.GetBlock(1), 0);
393
394  GridFunction x_ref(&fespace);
395  mesh->GetNodes(x_ref);
396
397  L2_FECollection w_fec(order + 1, dim);
398  FiniteElementSpace w_fespace(mesh, &w_fec);
399  GridFunction w(&w_fespace);
400
401  // 6. Set the initial conditions for v and x, and the boundary conditions on
402  // a beam-like mesh (see description above).
404  v.ProjectCoefficient(velo);
405  v.SetTrueVector();
407  x.ProjectCoefficient(deform);
408  x.SetTrueVector();
409
410  Array<int> ess_bdr(fespace.GetMesh()->bdr_attributes.Max());
411  ess_bdr = 0;
412  ess_bdr[0] = 1; // boundary attribute 1 (index 0) is fixed
413
414  // 7. Initialize the hyperelastic operator, the GLVis visualization and print
415  // the initial energies.
416  HyperelasticOperator oper(fespace, ess_bdr, visc, mu, K, nls_map[nls]);
417
418  socketstream vis_v, vis_w;
419  if (visualization)
420  {
421  char vishost[] = "localhost";
422  int visport = 19916;
423  vis_v.open(vishost, visport);
424  vis_v.precision(8);
426  visualize(vis_v, mesh, &x, &v, "Velocity", true);
427  vis_w.open(vishost, visport);
428  if (vis_w)
429  {
430  oper.GetElasticEnergyDensity(x, w);
431  vis_w.precision(8);
432  visualize(vis_w, mesh, &x, &w, "Elastic energy density", true);
433  }
434  }
435
436  double ee0 = oper.ElasticEnergy(x.GetTrueVector());
437  double ke0 = oper.KineticEnergy(v.GetTrueVector());
438  cout << "initial elastic energy (EE) = " << ee0 << endl;
439  cout << "initial kinetic energy (KE) = " << ke0 << endl;
440  cout << "initial total energy (TE) = " << (ee0 + ke0) << endl;
441
442  double t = 0.0;
443  oper.SetTime(t);
444  ode_solver->Init(oper);
445
446  // 8. Perform time-integration (looping over the time iterations, ti, with a
447  // time-step dt).
448  bool last_step = false;
449  for (int ti = 1; !last_step; ti++)
450  {
451  double dt_real = min(dt, t_final - t);
452
453  ode_solver->Step(vx, t, dt_real);
454
455  last_step = (t >= t_final - 1e-8*dt);
456
457  if (last_step || (ti % vis_steps) == 0)
458  {
459  double ee = oper.ElasticEnergy(x.GetTrueVector());
460  double ke = oper.KineticEnergy(v.GetTrueVector());
461
462  cout << "step " << ti << ", t = " << t << ", EE = " << ee << ", KE = "
463  << ke << ", ΔTE = " << (ee+ke)-(ee0+ke0) << endl;
464
465  if (cvode) { cvode->PrintInfo(); }
466  else if (arkode) { arkode->PrintInfo(); }
467
468  if (visualization)
469  {
471  visualize(vis_v, mesh, &x, &v);
472  if (vis_w)
473  {
474  oper.GetElasticEnergyDensity(x, w);
475  visualize(vis_w, mesh, &x, &w);
476  }
477  }
478  }
479  }
480
481  // 9. Save the displaced mesh, the velocity and elastic energy.
482  {
484  GridFunction *nodes = &x;
485  int owns_nodes = 0;
486  mesh->SwapNodes(nodes, owns_nodes);
487  ofstream mesh_ofs("deformed.mesh");
488  mesh_ofs.precision(8);
489  mesh->Print(mesh_ofs);
490  mesh->SwapNodes(nodes, owns_nodes);
491  ofstream velo_ofs("velocity.sol");
492  velo_ofs.precision(8);
493  v.Save(velo_ofs);
494  ofstream ee_ofs("elastic_energy.sol");
495  ee_ofs.precision(8);
496  oper.GetElasticEnergyDensity(x, w);
497  w.Save(ee_ofs);
498  }
499
500  // 10. Free the used memory.
501  delete ode_solver;
502  delete sjsolver;
503  delete mesh;
504
505  return 0;
506 }
507
508
509 void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
510  GridFunction *field, const char *field_name, bool init_vis)
511 {
512  if (!out)
513  {
514  return;
515  }
516
517  GridFunction *nodes = deformed_nodes;
518  int owns_nodes = 0;
519
520  mesh->SwapNodes(nodes, owns_nodes);
521
522  out << "solution\n" << *mesh << *field;
523
524  mesh->SwapNodes(nodes, owns_nodes);
525
526  if (init_vis)
527  {
528  out << "window_size 800 800\n";
529  out << "window_title '" << field_name << "'\n";
530  if (mesh->SpaceDimension() == 2)
531  {
532  out << "view 0 0\n"; // view from top
533  out << "keys jl\n"; // turn off perspective and light
534  }
535  out << "keys cm\n"; // show colorbar and mesh
536  out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
537  out << "pause\n";
538  }
539  out << flush;
540 }
541
542
543 ReducedSystemOperator::ReducedSystemOperator(
545  : Operator(M_->Height()), M(M_), S(S_), H(H_), Jacobian(NULL),
546  dt(0.0), v(NULL), x(NULL), w(height), z(height)
547 { }
548
549 void ReducedSystemOperator::SetParameters(double dt_, const Vector *v_,
550  const Vector *x_)
551 {
552  dt = dt_; v = v_; x = x_;
553 }
554
555 void ReducedSystemOperator::Mult(const Vector &k, Vector &y) const
556 {
557  // compute: y = H(x + dt*(v + dt*k)) + M*k + S*(v + dt*k)
560  H->Mult(z, y);
563 }
564
565 Operator &ReducedSystemOperator::GetGradient(const Vector &k) const
566 {
567  delete Jacobian;
568  Jacobian = Add(1.0, M->SpMat(), dt, S->SpMat());
573  return *Jacobian;
574 }
575
576 ReducedSystemOperator::~ReducedSystemOperator()
577 {
578  delete Jacobian;
579 }
580
581
582 int SundialsJacSolver::InitSystem(void *sundials_mem)
583 {
584  TimeDependentOperator *td_oper = GetTimeDependentOperator(sundials_mem);
585  HyperelasticOperator *he_oper;
586
587  // During development, we use dynamic_cast<> to ensure the setup is correct:
588  he_oper = dynamic_cast<HyperelasticOperator*>(td_oper);
589  MFEM_VERIFY(he_oper, "operator is not HyperelasticOperator");
590
591  // When the implementation is finalized, we can switch to static_cast<>:
592  // he_oper = static_cast<HyperelasticOperator*>(td_oper);
593
594  he_oper->InitSundialsJacSolver(*this);
595  return 0;
596 }
597
598 int SundialsJacSolver::SetupSystem(void *sundials_mem, int conv_fail,
599  const Vector &y_pred, const Vector &f_pred,
600  int &jac_cur, Vector &v_temp1,
601  Vector &v_temp2, Vector &v_temp3)
602 {
603  int sc = y_pred.Size() / 2;
604  const Vector x(y_pred.GetData() + sc, sc);
605  double dt = GetTimeStep(sundials_mem);
606
607  // J = M + dt*(S + dt*grad(H))
608  delete Jacobian;
609  Jacobian = Add(1.0, M->SpMat(), dt, S->SpMat());
612
613  J_solver->SetOperator(*Jacobian);
614
615  jac_cur = 1;
616  return 0;
617 }
618
619 int SundialsJacSolver::SolveSystem(void *sundials_mem, Vector &b,
620  const Vector &weight, const Vector &y_cur,
621  const Vector &f_cur)
622 {
623  int sc = b.Size() / 2;
624  // Vector x(y_cur.GetData() + sc, sc);
625  Vector b_v(b.GetData() + 0, sc);
626  Vector b_x(b.GetData() + sc, sc);
627  Vector rhs(sc);
628  double dt = GetTimeStep(sundials_mem);
629
630  // rhs = M b_v - dt*grad(H) b_x
632  rhs *= -dt;
634
635  J_solver->iterative_mode = false;
636  J_solver->Mult(rhs, b_v);
637
639
640  return 0;
641 }
642
643 int SundialsJacSolver::FreeSystem(void *sundials_mem)
644 {
645  delete Jacobian;
646  return 0;
647 }
648
649
650 HyperelasticOperator::HyperelasticOperator(FiniteElementSpace &f,
651  Array<int> &ess_bdr, double visc,
652  double mu, double K,
653  NonlinearSolverType nls_type)
654  : TimeDependentOperator(2*f.GetTrueVSize(), 0.0), fespace(f),
655  M(&fespace), S(&fespace), H(&fespace),
656  viscosity(visc), z(height/2)
657 {
658  const double rel_tol = 1e-8;
659  const int skip_zero_entries = 0;
660
661  const double ref_density = 1.0; // density in the reference configuration
662  ConstantCoefficient rho0(ref_density);
664  M.Assemble(skip_zero_entries);
665  Array<int> ess_tdof_list;
666  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
667  SparseMatrix tmp;
668  M.FormSystemMatrix(ess_tdof_list, tmp);
669
670  M_solver.iterative_mode = false;
671  M_solver.SetRelTol(rel_tol);
672  M_solver.SetAbsTol(0.0);
673  M_solver.SetMaxIter(30);
674  M_solver.SetPrintLevel(0);
675  M_solver.SetPreconditioner(M_prec);
676  M_solver.SetOperator(M.SpMat());
677
678  model = new NeoHookeanModel(mu, K);
680  H.SetEssentialTrueDofs(ess_tdof_list);
681
682  ConstantCoefficient visc_coeff(viscosity);
684  S.Assemble(skip_zero_entries);
685  S.FormSystemMatrix(ess_tdof_list, tmp);
686
687  reduced_oper = new ReducedSystemOperator(&M, &S, &H);
688
689 #ifndef MFEM_USE_SUITESPARSE
690  J_prec = new DSmoother(1);
691  MINRESSolver *J_minres = new MINRESSolver;
692  J_minres->SetRelTol(rel_tol);
693  J_minres->SetAbsTol(0.0);
694  J_minres->SetMaxIter(300);
695  J_minres->SetPrintLevel(-1);
696  J_minres->SetPreconditioner(*J_prec);
697  J_solver = J_minres;
698 #else
699  J_solver = new UMFPackSolver;
700  J_prec = NULL;
701 #endif
702
703  if (nls_type == KINSOL)
704  {
705  KinSolver *kinsolver = new KinSolver(KIN_NONE, true);
706  kinsolver->SetMaxSetupCalls(4);
707  newton_solver = kinsolver;
708  newton_solver->SetMaxIter(200);
709  newton_solver->SetRelTol(rel_tol);
710  newton_solver->SetPrintLevel(0);
711  }
712  else
713  {
714  newton_solver = new NewtonSolver();
715  newton_solver->SetMaxIter(10);
716  newton_solver->SetRelTol(rel_tol);
717  newton_solver->SetPrintLevel(-1);
718  }
719  newton_solver->SetSolver(*J_solver);
720  newton_solver->iterative_mode = false;
721  newton_solver->SetOperator(*reduced_oper);
722 }
723
724 void HyperelasticOperator::Mult(const Vector &vx, Vector &dvx_dt) const
725 {
726  // Create views to the sub-vectors v, x of vx, and dv_dt, dx_dt of dvx_dt
727  int sc = height/2;
728  Vector v(vx.GetData() + 0, sc);
729  Vector x(vx.GetData() + sc, sc);
730  Vector dv_dt(dvx_dt.GetData() + 0, sc);
731  Vector dx_dt(dvx_dt.GetData() + sc, sc);
732
733  H.Mult(x, z);
734  if (viscosity != 0.0)
735  {
737  }
738  z.Neg(); // z = -z
739  M_solver.Mult(z, dv_dt);
740
741  dx_dt = v;
742 }
743
744 void HyperelasticOperator::ImplicitSolve(const double dt,
745  const Vector &vx, Vector &dvx_dt)
746 {
747  int sc = height/2;
748  Vector v(vx.GetData() + 0, sc);
749  Vector x(vx.GetData() + sc, sc);
750  Vector dv_dt(dvx_dt.GetData() + 0, sc);
751  Vector dx_dt(dvx_dt.GetData() + sc, sc);
752
753  // By eliminating kx from the coupled system:
754  // kv = -M^{-1}*[H(x + dt*kx) + S*(v + dt*kv)]
755  // kx = v + dt*kv
756  // we reduce it to a nonlinear equation for kv, represented by the
757  // reduced_oper. This equation is solved with the newton_solver
758  // object (using J_solver and J_prec internally).
759  reduced_oper->SetParameters(dt, &v, &x);
760  Vector zero; // empty vector is interpreted as zero r.h.s. by NewtonSolver
761  newton_solver->Mult(zero, dv_dt);
762  MFEM_VERIFY(newton_solver->GetConverged(),
763  "Nonlinear solver did not converge.");
764 #ifdef MFEM_DEBUG
765  cout << " num nonlin sol iters = " << newton_solver->GetNumIterations()
766  << ", final norm = " << newton_solver->GetFinalNorm() << '\n';
767 #endif
769 }
770
771 void HyperelasticOperator::InitSundialsJacSolver(SundialsJacSolver &sjsolv)
772 {
773  sjsolv.SetOperators(M, S, H, *J_solver);
774 }
775
776 double HyperelasticOperator::ElasticEnergy(const Vector &x) const
777 {
778  return H.GetEnergy(x);
779 }
780
781 double HyperelasticOperator::KineticEnergy(const Vector &v) const
782 {
783  return 0.5*M.InnerProduct(v, v);
784 }
785
786 void HyperelasticOperator::GetElasticEnergyDensity(
787  const GridFunction &x, GridFunction &w) const
788 {
789  ElasticEnergyCoefficient w_coeff(*model, x);
790  w.ProjectCoefficient(w_coeff);
791 }
792
793 HyperelasticOperator::~HyperelasticOperator()
794 {
795  delete newton_solver;
796  delete J_solver;
797  delete J_prec;
798  delete reduced_oper;
799  delete model;
800 }
801
802
804  const IntegrationPoint &ip)
805 {
806  model.SetTransformation(T);
808  // return model.EvalW(J); // in reference configuration
809  return model.EvalW(J)/J.Det(); // in deformed configuration
810 }
811
812
813 void InitialDeformation(const Vector &x, Vector &y)
814 {
815  // set the initial configuration to be the same as the reference, stress
816  // free, configuration
817  y = x;
818 }
819
820 void InitialVelocity(const Vector &x, Vector &v)
821 {
822  const int dim = x.Size();
823  const double s = 0.1/64.;
824
825  v = 0.0;
826  v(dim-1) = s*x(0)*x(0)*(8.0-x(0));
827  v(0) = -s*x(0)*x(0);
828 }
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
Definition: ex10.cpp:376
double Eval(ElementTransformation &T, const IntegrationPoint &ip, double t)
Definition: coefficient.hpp:45
void Add(const int i, const int j, const double a)
Definition: sparsemat.cpp:2035
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
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
Data type for scaled Jacobi-type smoother of sparse matrix.
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:224
void SetMaxSetupCalls(int max_calls)
Set maximum number of nonlinear iterations without a Jacobian update.
Definition: sundials.cpp:980
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:122
Subclass constant coefficient.
Definition: coefficient.hpp:57
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.hpp:37
Base abstract class for time dependent operators.
Definition: operator.hpp:151
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
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:120
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
void SetLinearSolver(SundialsODELinearSolver &ls_spec)
Definition: sundials.cpp:233
MINRES method.
Definition: solvers.hpp:221
int main(int argc, char *argv[])
Definition: ex1.cpp:45
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:212
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:129
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:260
void SetSStolerances(double reltol, double abstol)
Specify the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:477
void InitialVelocity(const Vector &x, Vector &v)
Definition: ex10.cpp:595
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:344
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2349
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2942
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:116
int dim
Definition: ex3.cpp:47
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:234
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:72
Wrapper for SUNDIALS&#39; KINSOL library – Nonlinear solvers.
Definition: sundials.hpp:312
Data type sparse matrix.
Definition: sparsemat.hpp:38
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:224
Wrapper for SUNDIALS&#39; CVODE library – Multi-step time integration.
Definition: sundials.hpp:133
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6741
void MakeTRef(FiniteElementSpace *f, double *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
Definition: gridfunc.cpp:199
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
void SetMaxStep(double dt_max)
Set the maximum time step of the linear multistep method.
Definition: sundials.hpp:183
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:259
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:256
int Dimension() const
Definition: mesh.hpp:645
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
const Vector & GetTrueVector() const
Definition: gridfunc.hpp:105
Wrapper for SUNDIALS&#39; ARKODE library – Runge-Kutta time integration.
Definition: sundials.hpp:220
void SetLinearSolver(SundialsODELinearSolver &ls_spec)
Set a custom Jacobian system solver for implicit methods.
Definition: sundials.cpp:486
int SpaceDimension() const
Definition: mesh.hpp:646
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:147
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
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
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:134
Definition: gridfunc.cpp:1040
Abstract base class, wrapping the custom linear solvers interface in SUNDIALS&#39; CVODE and ARKODE solve...
Definition: sundials.hpp:47
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:225
Class for integration point with weight.
Definition: intrules.hpp:25
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:1377
Abstract class for hyperelastic models.
int open(const char hostname[], int port)
Vector data type.
Definition: vector.hpp:48
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5422
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:79
void PrintInfo() const
Print CVODE statistics.
Definition: sundials.cpp:400
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
The classical forward Euler method.
Definition: ode.hpp:101
void SetMaxStep(double dt_max)
Set the maximum time step of the Runge-Kutta method.
Definition: sundials.hpp:276
Abstract operator.
Definition: operator.hpp:21
void PrintInfo() const
Print ARKODE statistics.
Definition: sundials.cpp:691
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:77
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:128
void Neg()
(*this) = -(*this)
Definition: vector.cpp:252
bool Good() const
Definition: optparser.hpp:120