MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex19p.cpp
Go to the documentation of this file.
1 // MFEM Example 19 - Parallel Version
2 //
3 // Compile with: make ex19p
4 //
5 // Sample runs:
6 // mpirun -np 2 ex19p -m ../data/beam-quad.mesh
7 // mpirun -np 2 ex19p -m ../data/beam-tri.mesh
8 // mpirun -np 2 ex19p -m ../data/beam-hex.mesh
9 // mpirun -np 2 ex19p -m ../data/beam-tet.mesh
10 // mpirun -np 2 ex19p -m ../data/beam-wedge.mesh
11 // mpirun -np 2 ex19p -m ../data/beam-quad-amr.mesh
12 //
13 // Description: This examples solves a quasi-static incompressible nonlinear
14 // elasticity problem of the form 0 = H(x), where H is an
15 // incompressible hyperelastic model and x is a block state vector
16 // containing displacement and pressure variables. The geometry of
17 // the domain is assumed to be as follows:
18 //
19 // +---------------------+
20 // boundary --->| |<--- boundary
21 // attribute 1 | | attribute 2
22 // (fixed) +---------------------+ (fixed, nonzero)
23 //
24 // The example demonstrates the use of block nonlinear operators
25 // (the class RubberOperator defining H(x)) as well as a nonlinear
26 // Newton solver for the quasi-static problem. Each Newton step
27 // requires the inversion of a Jacobian matrix, which is done
28 // through a (preconditioned) inner solver. The specialized block
29 // preconditioner is implemented as a user-defined solver.
30 //
31 // We recommend viewing examples 2, 5, and 10 before viewing this
32 // example.
33 
34 #include "mfem.hpp"
35 #include <memory>
36 #include <iostream>
37 #include <fstream>
38 
39 using namespace std;
40 using namespace mfem;
41 
42 class GeneralResidualMonitor : public IterativeSolverMonitor
43 {
44 public:
45  GeneralResidualMonitor(MPI_Comm comm, const std::string& prefix_,
46  int print_lvl)
47  : prefix(prefix_)
48  {
49 #ifndef MFEM_USE_MPI
50  print_level = print_lvl;
51 #else
52  int rank;
53  MPI_Comm_rank(comm, &rank);
54  if (rank == 0)
55  {
56  print_level = print_lvl;
57  }
58  else
59  {
60  print_level = -1;
61  }
62 #endif
63  }
64 
65  virtual void MonitorResidual(int it, double norm, const Vector &r, bool final);
66 
67 private:
68  const std::string prefix;
69  int print_level;
70  mutable double norm0;
71 };
72 
73 void GeneralResidualMonitor::MonitorResidual(int it, double norm,
74  const Vector &r, bool final)
75 {
76  if (print_level == 1 || (print_level == 3 && (final || it == 0)))
77  {
78  mfem::out << prefix << " iteration " << setw(2) << it
79  << " : ||r|| = " << norm;
80  if (it > 0)
81  {
82  mfem::out << ", ||r||/||r_0|| = " << norm/norm0;
83  }
84  else
85  {
86  norm0 = norm;
87  }
88  mfem::out << '\n';
89  }
90 }
91 
92 // Custom block preconditioner for the Jacobian of the incompressible nonlinear
93 // elasticity operator. It has the form
94 //
95 // P^-1 = [ K^-1 0 ][ I -B^T ][ I 0 ]
96 // [ 0 I ][ 0 I ][ 0 -\gamma S^-1 ]
97 //
98 // where the original Jacobian has the form
99 //
100 // J = [ K B^T ]
101 // [ B 0 ]
102 //
103 // and K^-1 is an approximation of the inverse of the displacement part of the
104 // Jacobian and S^-1 is an approximation of the inverse of the Schur
105 // complement S = B K^-1 B^T. The Schur complement is approximated using
106 // a mass matrix of the pressure variables.
107 class JacobianPreconditioner : public Solver
108 {
109 protected:
110  // Finite element spaces for setting up preconditioner blocks
112 
113  // Offsets for extracting block vector segments
114  Array<int> &block_trueOffsets;
115 
116  // Jacobian for block access
117  BlockOperator *jacobian;
118 
119  // Scaling factor for the pressure mass matrix in the block preconditioner
120  double gamma;
121 
122  // Objects for the block preconditioner application
123  Operator *pressure_mass;
124  Solver *mass_pcg;
125  Solver *mass_prec;
126  Solver *stiff_pcg;
127  Solver *stiff_prec;
128 
129 public:
130  JacobianPreconditioner(Array<ParFiniteElementSpace *> &fes,
131  Operator &mass, Array<int> &offsets);
132 
133  virtual void Mult(const Vector &k, Vector &y) const;
134  virtual void SetOperator(const Operator &op);
135 
136  virtual ~JacobianPreconditioner();
137 };
138 
139 // After spatial discretization, the rubber model can be written as:
140 // 0 = H(x)
141 // where x is the block vector representing the deformation and pressure and
142 // H(x) is the nonlinear incompressible neo-Hookean operator.
143 class RubberOperator : public Operator
144 {
145 protected:
146  // Finite element spaces
148 
149  // Block nonlinear form
150  ParBlockNonlinearForm *Hform;
151 
152  // Pressure mass matrix for the preconditioner
153  Operator *pressure_mass;
154 
155  // Newton solver for the hyperelastic operator
156  NewtonSolver newton_solver;
157  GeneralResidualMonitor newton_monitor;
158 
159  // Solver for the Jacobian solve in the Newton method
160  Solver *j_solver;
161  GeneralResidualMonitor j_monitor;
162 
163  // Preconditioner for the Jacobian
164  Solver *j_prec;
165 
166  // Shear modulus coefficient
167  Coefficient &mu;
168 
169  // Block offsets for variable access
170  Array<int> &block_trueOffsets;
171 
172 public:
173  RubberOperator(Array<ParFiniteElementSpace *> &fes, Array<Array<int> *>&ess_bdr,
174  Array<int> &block_trueOffsets, double rel_tol, double abs_tol,
175  int iter, Coefficient &mu);
176 
177  // Required to use the native newton solver
178  virtual Operator &GetGradient(const Vector &xp) const;
179  virtual void Mult(const Vector &k, Vector &y) const;
180 
181  // Driver for the newton solver
182  void Solve(Vector &xp) const;
183 
184  virtual ~RubberOperator();
185 };
186 
187 // Visualization driver
188 void visualize(ostream &os, ParMesh *mesh,
189  ParGridFunction *deformed_nodes,
190  ParGridFunction *field, const char *field_name = NULL,
191  bool init_vis = false);
192 
193 // Configuration definition functions
194 void ReferenceConfiguration(const Vector &x, Vector &y);
195 void InitialDeformation(const Vector &x, Vector &y);
196 
197 
198 int main(int argc, char *argv[])
199 {
200 #ifdef HYPRE_USING_GPU
201  cout << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this example\n"
202  << "is NOT supported with the GPU version of hypre.\n\n";
203  return 242;
204 #endif
205 
206  // 1. Initialize MPI and HYPRE.
207  Mpi::Init();
208  const int myid = Mpi::WorldRank();
209  Hypre::Init();
210 
211  // 2. Parse command-line options
212  const char *mesh_file = "../data/beam-tet.mesh";
213  int ser_ref_levels = 0;
214  int par_ref_levels = 0;
215  int order = 2;
216  bool visualization = true;
217  double newton_rel_tol = 1e-4;
218  double newton_abs_tol = 1e-6;
219  int newton_iter = 500;
220  double mu = 1.0;
221 
222  OptionsParser args(argc, argv);
223  args.AddOption(&mesh_file, "-m", "--mesh",
224  "Mesh file to use.");
225  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
226  "Number of times to refine the mesh uniformly in serial.");
227  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
228  "Number of times to refine the mesh uniformly in parallel.");
229  args.AddOption(&order, "-o", "--order",
230  "Order (degree) of the finite elements.");
231  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
232  "--no-visualization",
233  "Enable or disable GLVis visualization.");
234  args.AddOption(&newton_rel_tol, "-rel", "--relative-tolerance",
235  "Relative tolerance for the Newton solve.");
236  args.AddOption(&newton_abs_tol, "-abs", "--absolute-tolerance",
237  "Absolute tolerance for the Newton solve.");
238  args.AddOption(&newton_iter, "-it", "--newton-iterations",
239  "Maximum iterations for the Newton solve.");
240  args.AddOption(&mu, "-mu", "--shear-modulus",
241  "Shear modulus for the neo-Hookean material.");
242  args.Parse();
243  if (!args.Good())
244  {
245  if (myid == 0)
246  {
247  args.PrintUsage(cout);
248  }
249  return 1;
250  }
251  if (myid == 0)
252  {
253  args.PrintOptions(cout);
254  }
255 
256  // 3. Read the (serial) mesh from the given mesh file on all processors. We
257  // can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
258  // with the same code.
259  Mesh *mesh = new Mesh(mesh_file, 1, 1);
260  int dim = mesh->Dimension();
261 
262  // 4. Refine the mesh in serial to increase the resolution. In this example
263  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
264  // a command-line parameter.
265  for (int lev = 0; lev < ser_ref_levels; lev++)
266  {
267  mesh->UniformRefinement();
268  }
269 
270  // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
271  // this mesh further in parallel to increase the resolution. Once the
272  // parallel mesh is defined, the serial mesh can be deleted.
273  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
274  delete mesh;
275  for (int lev = 0; lev < par_ref_levels; lev++)
276  {
277  pmesh->UniformRefinement();
278  }
279 
280  // 6. Define the shear modulus for the incompressible Neo-Hookean material
281  ConstantCoefficient c_mu(mu);
282 
283  // 7. Define the finite element spaces for displacement and pressure
284  // (Taylor-Hood elements). By default, the displacement (u/x) is a second
285  // order vector field, while the pressure (p) is a linear scalar function.
286  H1_FECollection quad_coll(order, dim);
287  H1_FECollection lin_coll(order-1, dim);
288 
289  ParFiniteElementSpace R_space(pmesh, &quad_coll, dim, Ordering::byVDIM);
290  ParFiniteElementSpace W_space(pmesh, &lin_coll);
291 
293  spaces[0] = &R_space;
294  spaces[1] = &W_space;
295 
296  HYPRE_BigInt glob_R_size = R_space.GlobalTrueVSize();
297  HYPRE_BigInt glob_W_size = W_space.GlobalTrueVSize();
298 
299  // 8. Define the Dirichlet conditions (set to boundary attribute 1 and 2)
300  Array<Array<int> *> ess_bdr(2);
301 
302  Array<int> ess_bdr_u(R_space.GetMesh()->bdr_attributes.Max());
303  Array<int> ess_bdr_p(W_space.GetMesh()->bdr_attributes.Max());
304 
305  ess_bdr_p = 0;
306  ess_bdr_u = 0;
307  ess_bdr_u[0] = 1;
308  ess_bdr_u[1] = 1;
309 
310  ess_bdr[0] = &ess_bdr_u;
311  ess_bdr[1] = &ess_bdr_p;
312 
313  // 9. Print the mesh statistics
314  if (myid == 0)
315  {
316  std::cout << "***********************************************************\n";
317  std::cout << "dim(u) = " << glob_R_size << "\n";
318  std::cout << "dim(p) = " << glob_W_size << "\n";
319  std::cout << "dim(u+p) = " << glob_R_size + glob_W_size << "\n";
320  std::cout << "***********************************************************\n";
321  }
322 
323  // 10. Define the block structure of the solution vector (u then p)
324  Array<int> block_trueOffsets(3);
325  block_trueOffsets[0] = 0;
326  block_trueOffsets[1] = R_space.TrueVSize();
327  block_trueOffsets[2] = W_space.TrueVSize();
328  block_trueOffsets.PartialSum();
329 
330  BlockVector xp(block_trueOffsets);
331 
332  // 11. Define grid functions for the current configuration, reference
333  // configuration, final deformation, and pressure
334  ParGridFunction x_gf(&R_space);
335  ParGridFunction x_ref(&R_space);
336  ParGridFunction x_def(&R_space);
337  ParGridFunction p_gf(&W_space);
338 
341 
342  x_gf.ProjectCoefficient(deform);
343  x_ref.ProjectCoefficient(refconfig);
344  p_gf = 0.0;
345 
346  // 12. Set up the block solution vectors
347  x_gf.GetTrueDofs(xp.GetBlock(0));
348  p_gf.GetTrueDofs(xp.GetBlock(1));
349 
350  // 13. Initialize the incompressible neo-Hookean operator
351  RubberOperator oper(spaces, ess_bdr, block_trueOffsets,
352  newton_rel_tol, newton_abs_tol, newton_iter, c_mu);
353 
354  // 14. Solve the Newton system
355  oper.Solve(xp);
356 
357  // 15. Distribute the shared degrees of freedom
358  x_gf.Distribute(xp.GetBlock(0));
359  p_gf.Distribute(xp.GetBlock(1));
360 
361  // 16. Compute the final deformation
362  subtract(x_gf, x_ref, x_def);
363 
364  // 17. Visualize the results if requested
365  socketstream vis_u, vis_p;
366  if (visualization)
367  {
368  char vishost[] = "localhost";
369  int visport = 19916;
370  vis_u.open(vishost, visport);
371  vis_u.precision(8);
372  visualize(vis_u, pmesh, &x_gf, &x_def, "Deformation", true);
373  // Make sure all ranks have sent their 'u' solution before initiating
374  // another set of GLVis connections (one from each rank):
375  MPI_Barrier(pmesh->GetComm());
376  vis_p.open(vishost, visport);
377  vis_p.precision(8);
378  visualize(vis_p, pmesh, &x_gf, &p_gf, "Pressure", true);
379  }
380 
381  // 18. Save the displaced mesh, the final deformation, and the pressure
382  {
383  GridFunction *nodes = &x_gf;
384  int owns_nodes = 0;
385  pmesh->SwapNodes(nodes, owns_nodes);
386 
387  ostringstream mesh_name, pressure_name, deformation_name;
388  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
389  pressure_name << "pressure." << setfill('0') << setw(6) << myid;
390  deformation_name << "deformation." << setfill('0') << setw(6) << myid;
391 
392  ofstream mesh_ofs(mesh_name.str().c_str());
393  mesh_ofs.precision(8);
394  pmesh->Print(mesh_ofs);
395 
396  ofstream pressure_ofs(pressure_name.str().c_str());
397  pressure_ofs.precision(8);
398  p_gf.Save(pressure_ofs);
399 
400  ofstream deformation_ofs(deformation_name.str().c_str());
401  deformation_ofs.precision(8);
402  x_def.Save(deformation_ofs);
403  }
404 
405  // 19. Free the used memory
406  delete pmesh;
407 
408  return 0;
409 }
410 
411 
412 JacobianPreconditioner::JacobianPreconditioner(Array<ParFiniteElementSpace *>
413  &fes,
414  Operator &mass,
415  Array<int> &offsets)
416  : Solver(offsets[2]), block_trueOffsets(offsets), pressure_mass(&mass)
417 {
418  fes.Copy(spaces);
419 
420  gamma = 0.00001;
421 
422  // The mass matrix and preconditioner do not change every Newton cycle, so
423  // we only need to define them once
424  HypreBoomerAMG *mass_prec_amg = new HypreBoomerAMG();
425  mass_prec_amg->SetPrintLevel(0);
426 
427  mass_prec = mass_prec_amg;
428 
429  CGSolver *mass_pcg_iter = new CGSolver(spaces[0]->GetComm());
430  mass_pcg_iter->SetRelTol(1e-12);
431  mass_pcg_iter->SetAbsTol(1e-12);
432  mass_pcg_iter->SetMaxIter(200);
433  mass_pcg_iter->SetPrintLevel(0);
434  mass_pcg_iter->SetPreconditioner(*mass_prec);
435  mass_pcg_iter->SetOperator(*pressure_mass);
436  mass_pcg_iter->iterative_mode = false;
437 
438  mass_pcg = mass_pcg_iter;
439 
440  // The stiffness matrix does change every Newton cycle, so we will define it
441  // during SetOperator
442  stiff_pcg = NULL;
443  stiff_prec = NULL;
444 }
445 
446 void JacobianPreconditioner::Mult(const Vector &k, Vector &y) const
447 {
448  // Extract the blocks from the input and output vectors
449  Vector disp_in;
450  disp_in.MakeRef(const_cast<Vector&>(k), block_trueOffsets[0],
451  block_trueOffsets[1]-block_trueOffsets[0]);
452  Vector pres_in;
453  pres_in.MakeRef(const_cast<Vector&>(k), block_trueOffsets[1],
454  block_trueOffsets[2]-block_trueOffsets[1]);
455 
456  Vector disp_out;
457  disp_out.MakeRef(y, block_trueOffsets[0],
458  block_trueOffsets[1]-block_trueOffsets[0]);
459  Vector pres_out;
460  pres_out.MakeRef(y, block_trueOffsets[1],
461  block_trueOffsets[2]-block_trueOffsets[1]);
462 
463  Vector temp(block_trueOffsets[1]-block_trueOffsets[0]);
464  Vector temp2(block_trueOffsets[1]-block_trueOffsets[0]);
465 
466  // Perform the block elimination for the preconditioner
467  mass_pcg->Mult(pres_in, pres_out);
468  pres_out *= -gamma;
469 
470  jacobian->GetBlock(0,1).Mult(pres_out, temp);
471  subtract(disp_in, temp, temp2);
472 
473  stiff_pcg->Mult(temp2, disp_out);
474 
475  disp_out.SyncAliasMemory(y);
476  pres_out.SyncAliasMemory(y);
477 }
478 
479 void JacobianPreconditioner::SetOperator(const Operator &op)
480 {
481  jacobian = (BlockOperator *) &op;
482 
483  // Initialize the stiffness preconditioner and solver
484  if (stiff_prec == NULL)
485  {
486  HypreBoomerAMG *stiff_prec_amg = new HypreBoomerAMG();
487  stiff_prec_amg->SetPrintLevel(0);
488 
489  if (!spaces[0]->GetParMesh()->Nonconforming())
490  {
491 #if !defined(HYPRE_USING_GPU)
492  // Not available yet when hypre is built with GPU support
493  stiff_prec_amg->SetElasticityOptions(spaces[0]);
494 #endif
495  }
496 
497  stiff_prec = stiff_prec_amg;
498 
499  GMRESSolver *stiff_pcg_iter = new GMRESSolver(spaces[0]->GetComm());
500  stiff_pcg_iter->SetRelTol(1e-8);
501  stiff_pcg_iter->SetAbsTol(1e-8);
502  stiff_pcg_iter->SetMaxIter(200);
503  stiff_pcg_iter->SetPrintLevel(0);
504  stiff_pcg_iter->SetPreconditioner(*stiff_prec);
505  stiff_pcg_iter->iterative_mode = false;
506 
507  stiff_pcg = stiff_pcg_iter;
508  }
509 
510  // At each Newton cycle, compute the new stiffness AMG preconditioner by
511  // updating the iterative solver which, in turn, updates its preconditioner
512  stiff_pcg->SetOperator(jacobian->GetBlock(0,0));
513 }
514 
515 JacobianPreconditioner::~JacobianPreconditioner()
516 {
517  delete mass_pcg;
518  delete mass_prec;
519  delete stiff_prec;
520  delete stiff_pcg;
521 }
522 
523 
524 RubberOperator::RubberOperator(Array<ParFiniteElementSpace *> &fes,
525  Array<Array<int> *> &ess_bdr,
526  Array<int> &trueOffsets,
527  double rel_tol,
528  double abs_tol,
529  int iter,
530  Coefficient &c_mu)
531  : Operator(fes[0]->TrueVSize() + fes[1]->TrueVSize()),
532  newton_solver(fes[0]->GetComm()),
533  newton_monitor(fes[0]->GetComm(), "Newton", 1),
534  j_monitor(fes[0]->GetComm(), " GMRES", 3),
535  mu(c_mu), block_trueOffsets(trueOffsets)
536 {
537  Array<Vector *> rhs(2);
538  rhs = NULL; // Set all entries in the array
539 
540  fes.Copy(spaces);
541 
542  // Define the block nonlinear form
543  Hform = new ParBlockNonlinearForm(spaces);
544 
545  // Add the incompressible neo-Hookean integrator
546  Hform->AddDomainIntegrator(new IncompressibleNeoHookeanIntegrator(mu));
547 
548  // Set the essential boundary conditions
549  Hform->SetEssentialBC(ess_bdr, rhs);
550 
551  // Compute the pressure mass stiffness matrix
552  ParBilinearForm *a = new ParBilinearForm(spaces[1]);
553  ConstantCoefficient one(1.0);
554  OperatorHandle mass(Operator::Hypre_ParCSR);
555  a->AddDomainIntegrator(new MassIntegrator(one));
556  a->Assemble();
557  a->Finalize();
558  a->ParallelAssemble(mass);
559  delete a;
560 
561  mass.SetOperatorOwner(false);
562  pressure_mass = mass.Ptr();
563 
564  // Initialize the Jacobian preconditioner
565  JacobianPreconditioner *jac_prec =
566  new JacobianPreconditioner(fes, *pressure_mass, block_trueOffsets);
567  j_prec = jac_prec;
568 
569  // Set up the Jacobian solver
570  GMRESSolver *j_gmres = new GMRESSolver(spaces[0]->GetComm());
571  j_gmres->iterative_mode = false;
572  j_gmres->SetRelTol(1e-12);
573  j_gmres->SetAbsTol(1e-12);
574  j_gmres->SetMaxIter(300);
575  j_gmres->SetPrintLevel(-1);
576  j_gmres->SetMonitor(j_monitor);
577  j_gmres->SetPreconditioner(*j_prec);
578  j_solver = j_gmres;
579 
580  // Set the newton solve parameters
581  newton_solver.iterative_mode = true;
582  newton_solver.SetSolver(*j_solver);
583  newton_solver.SetOperator(*this);
584  newton_solver.SetPrintLevel(-1);
585  newton_solver.SetMonitor(newton_monitor);
586  newton_solver.SetRelTol(rel_tol);
587  newton_solver.SetAbsTol(abs_tol);
588  newton_solver.SetMaxIter(iter);
589 }
590 
591 // Solve the Newton system
592 void RubberOperator::Solve(Vector &xp) const
593 {
594  Vector zero;
595  newton_solver.Mult(zero, xp);
596  MFEM_VERIFY(newton_solver.GetConverged(),
597  "Newton Solver did not converge.");
598 }
599 
600 // compute: y = H(x,p)
601 void RubberOperator::Mult(const Vector &k, Vector &y) const
602 {
603  Hform->Mult(k, y);
604 }
605 
606 // Compute the Jacobian from the nonlinear form
607 Operator &RubberOperator::GetGradient(const Vector &xp) const
608 {
609  return Hform->GetGradient(xp);
610 }
611 
612 RubberOperator::~RubberOperator()
613 {
614  delete Hform;
615  delete pressure_mass;
616  delete j_solver;
617  delete j_prec;
618 }
619 
620 
621 // Inline visualization
622 void visualize(ostream &os, ParMesh *mesh,
623  ParGridFunction *deformed_nodes,
624  ParGridFunction *field, const char *field_name, bool init_vis)
625 {
626  if (!os)
627  {
628  return;
629  }
630 
631  GridFunction *nodes = deformed_nodes;
632  int owns_nodes = 0;
633 
634  mesh->SwapNodes(nodes, owns_nodes);
635 
636  os << "parallel " << mesh->GetNRanks() << " " << mesh->GetMyRank() <<
637  "\n";
638  os << "solution\n" << *mesh << *field;
639 
640  mesh->SwapNodes(nodes, owns_nodes);
641 
642  if (init_vis)
643  {
644  os << "window_size 800 800\n";
645  os << "window_title '" << field_name << "'\n";
646  if (mesh->SpaceDimension() == 2)
647  {
648  os << "view 0 0\n"; // view from top
649  // turn off perspective and light, +anti-aliasing
650  os << "keys jlA\n";
651  }
652  os << "keys cmA\n"; // show colorbar and mesh, +anti-aliasing
653  // update value-range; keep mesh-extents fixed
654  os << "autoscale value\n";
655  }
656  os << flush;
657 }
658 
659 void ReferenceConfiguration(const Vector &x, Vector &y)
660 {
661  // Set the reference, stress free, configuration
662  y = x;
663 }
664 
665 void InitialDeformation(const Vector &x, Vector &y)
666 {
667  // Set the initial configuration. Having this different from the reference
668  // configuration can help convergence
669  y = x;
670  y[1] = x[1] + 0.25*x[0];
671 }
void InitialDeformation(const Vector &x, Vector &y)
Definition: ex10.cpp:594
Conjugate gradient method.
Definition: solvers.hpp:465
void ReferenceConfiguration(const Vector &x, Vector &y)
Definition: ex19.cpp:579
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
int TrueVSize() const
Obsolete, kept for backward compatibility.
Definition: pfespace.hpp:434
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:99
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Definition: mesh.cpp:7957
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:856
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:873
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:659
void SyncAliasMemory(const Vector &v) const
Update the alias memory location of the vector to match v.
Definition: vector.hpp:245
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
int GetNRanks() const
Definition: pmesh.hpp:352
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:179
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1579
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
constexpr char vishost[]
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9816
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1659
constexpr int visport
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void visualize(ostream &os, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
Definition: ex10.cpp:381
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:613
void SetElasticityOptions(ParFiniteElementSpace *fespace)
Definition: hypre.cpp:5083
int Dimension() const
Definition: mesh.hpp:1006
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
A general vector function coefficient.
int SpaceDimension() const
Definition: mesh.hpp:1007
Abstract base class for an iterative solver monitor.
Definition: solvers.hpp:36
void SetAbsTol(double atol)
Definition: solvers.hpp:199
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
int GetMyRank() const
Definition: pmesh.hpp:353
void SetRelTol(double rtol)
Definition: solvers.hpp:198
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:460
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)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition: array.cpp:103
HYPRE_Int HYPRE_BigInt
GMRES method.
Definition: solvers.hpp:497
void Distribute(const Vector *tv)
Definition: pgridfunc.cpp:132
A class representing a general parallel block nonlinear operator defined on the Cartesian product of ...
void SetMonitor(IterativeSolverMonitor &m)
Set the iterative solver monitor.
Definition: solvers.hpp:258
double a
Definition: lissajous.cpp:41
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
double mu
Definition: ex25.cpp:139
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:143
Class for parallel bilinear form.
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:479
Vector data type.
Definition: vector.hpp:60
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:577
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4770
Base class for solvers.
Definition: operator.hpp:655
Class for parallel grid function.
Definition: pgridfunc.hpp:32
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:66
Abstract operator.
Definition: operator.hpp:24
A class to handle Block systems in a matrix-free implementation.
Class for parallel meshes.
Definition: pmesh.hpp:32
int main()
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:90
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150