MFEM  v4.2.0 Finite element discretization library
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 &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
189  ParGridFunction *field, const char *field_name = NULL,
190  bool init_vis = false);
191
192 // Configuration definition functions
193 void ReferenceConfiguration(const Vector &x, Vector &y);
194 void InitialDeformation(const Vector &x, Vector &y);
195
196
197 int main(int argc, char *argv[])
198 {
199  // 1. Initialize MPI
200  MPI_Session mpi;
201  const int myid = mpi.WorldRank();
202
203  // 2. Parse command-line options
204  const char *mesh_file = "../data/beam-tet.mesh";
205  int ser_ref_levels = 0;
206  int par_ref_levels = 0;
207  int order = 2;
208  bool visualization = true;
209  double newton_rel_tol = 1e-4;
210  double newton_abs_tol = 1e-6;
211  int newton_iter = 500;
212  double mu = 1.0;
213
214  OptionsParser args(argc, argv);
216  "Mesh file to use.");
218  "Number of times to refine the mesh uniformly in serial.");
220  "Number of times to refine the mesh uniformly in parallel.");
222  "Order (degree) of the finite elements.");
224  "--no-visualization",
225  "Enable or disable GLVis visualization.");
227  "Relative tolerance for the Newton solve.");
229  "Absolute tolerance for the Newton solve.");
231  "Maximum iterations for the Newton solve.");
233  "Shear modulus for the neo-Hookean material.");
234  args.Parse();
235  if (!args.Good())
236  {
237  if (myid == 0)
238  {
239  args.PrintUsage(cout);
240  }
241  return 1;
242  }
243  if (myid == 0)
244  {
245  args.PrintOptions(cout);
246  }
247
248  // 3. Read the (serial) mesh from the given mesh file on all processors. We
249  // can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
250  // with the same code.
251  Mesh *mesh = new Mesh(mesh_file, 1, 1);
252  int dim = mesh->Dimension();
253
254  // 4. Refine the mesh in serial to increase the resolution. In this example
255  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
256  // a command-line parameter.
257  for (int lev = 0; lev < ser_ref_levels; lev++)
258  {
259  mesh->UniformRefinement();
260  }
261
262  // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
263  // this mesh further in parallel to increase the resolution. Once the
264  // parallel mesh is defined, the serial mesh can be deleted.
265  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
266  delete mesh;
267  for (int lev = 0; lev < par_ref_levels; lev++)
268  {
269  pmesh->UniformRefinement();
270  }
271
272  // 6. Define the shear modulus for the incompressible Neo-Hookean material
273  ConstantCoefficient c_mu(mu);
274
275  // 7. Define the finite element spaces for displacement and pressure
276  // (Taylor-Hood elements). By default, the displacement (u/x) is a second
277  // order vector field, while the pressure (p) is a linear scalar function.
279  H1_FECollection lin_coll(order-1, dim);
280
281  ParFiniteElementSpace R_space(pmesh, &quad_coll, dim, Ordering::byVDIM);
282  ParFiniteElementSpace W_space(pmesh, &lin_coll);
283
285  spaces[0] = &R_space;
286  spaces[1] = &W_space;
287
288  HYPRE_Int glob_R_size = R_space.GlobalTrueVSize();
289  HYPRE_Int glob_W_size = W_space.GlobalTrueVSize();
290
291  // 8. Define the Dirichlet conditions (set to boundary attribute 1 and 2)
292  Array<Array<int> *> ess_bdr(2);
293
294  Array<int> ess_bdr_u(R_space.GetMesh()->bdr_attributes.Max());
295  Array<int> ess_bdr_p(W_space.GetMesh()->bdr_attributes.Max());
296
297  ess_bdr_p = 0;
298  ess_bdr_u = 0;
299  ess_bdr_u[0] = 1;
300  ess_bdr_u[1] = 1;
301
302  ess_bdr[0] = &ess_bdr_u;
303  ess_bdr[1] = &ess_bdr_p;
304
305  // 9. Print the mesh statistics
306  if (myid == 0)
307  {
308  std::cout << "***********************************************************\n";
309  std::cout << "dim(u) = " << glob_R_size << "\n";
310  std::cout << "dim(p) = " << glob_W_size << "\n";
311  std::cout << "dim(u+p) = " << glob_R_size + glob_W_size << "\n";
312  std::cout << "***********************************************************\n";
313  }
314
315  // 10. Define the block structure of the solution vector (u then p)
316  Array<int> block_trueOffsets(3);
317  block_trueOffsets[0] = 0;
318  block_trueOffsets[1] = R_space.TrueVSize();
319  block_trueOffsets[2] = W_space.TrueVSize();
320  block_trueOffsets.PartialSum();
321
322  BlockVector xp(block_trueOffsets);
323
324  // 11. Define grid functions for the current configuration, reference
325  // configuration, final deformation, and pressure
326  ParGridFunction x_gf(&R_space);
327  ParGridFunction x_ref(&R_space);
328  ParGridFunction x_def(&R_space);
329  ParGridFunction p_gf(&W_space);
330
333
334  x_gf.ProjectCoefficient(deform);
335  x_ref.ProjectCoefficient(refconfig);
336  p_gf = 0.0;
337
338  // 12. Set up the block solution vectors
339  x_gf.GetTrueDofs(xp.GetBlock(0));
340  p_gf.GetTrueDofs(xp.GetBlock(1));
341
342  // 13. Initialize the incompressible neo-Hookean operator
343  RubberOperator oper(spaces, ess_bdr, block_trueOffsets,
344  newton_rel_tol, newton_abs_tol, newton_iter, c_mu);
345
346  // 14. Solve the Newton system
347  oper.Solve(xp);
348
349  // 15. Distribute the shared degrees of freedom
350  x_gf.Distribute(xp.GetBlock(0));
351  p_gf.Distribute(xp.GetBlock(1));
352
353  // 16. Compute the final deformation
354  subtract(x_gf, x_ref, x_def);
355
356  // 17. Visualize the results if requested
357  socketstream vis_u, vis_p;
358  if (visualization)
359  {
360  char vishost[] = "localhost";
361  int visport = 19916;
362  vis_u.open(vishost, visport);
363  vis_u.precision(8);
364  visualize(vis_u, pmesh, &x_gf, &x_def, "Deformation", true);
365  // Make sure all ranks have sent their 'u' solution before initiating
366  // another set of GLVis connections (one from each rank):
367  MPI_Barrier(pmesh->GetComm());
368  vis_p.open(vishost, visport);
369  vis_p.precision(8);
370  visualize(vis_p, pmesh, &x_gf, &p_gf, "Pressure", true);
371  }
372
373  // 18. Save the displaced mesh, the final deformation, and the pressure
374  {
375  GridFunction *nodes = &x_gf;
376  int owns_nodes = 0;
377  pmesh->SwapNodes(nodes, owns_nodes);
378
379  ostringstream mesh_name, pressure_name, deformation_name;
380  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
381  pressure_name << "pressure." << setfill('0') << setw(6) << myid;
382  deformation_name << "deformation." << setfill('0') << setw(6) << myid;
383
384  ofstream mesh_ofs(mesh_name.str().c_str());
385  mesh_ofs.precision(8);
386  pmesh->Print(mesh_ofs);
387
388  ofstream pressure_ofs(pressure_name.str().c_str());
389  pressure_ofs.precision(8);
390  p_gf.Save(pressure_ofs);
391
392  ofstream deformation_ofs(deformation_name.str().c_str());
393  deformation_ofs.precision(8);
394  x_def.Save(deformation_ofs);
395  }
396
397  // 19. Free the used memory
398  delete pmesh;
399
400  return 0;
401 }
402
403
404 JacobianPreconditioner::JacobianPreconditioner(Array<ParFiniteElementSpace *>
405  &fes,
406  Operator &mass,
407  Array<int> &offsets)
408  : Solver(offsets[2]), block_trueOffsets(offsets), pressure_mass(&mass)
409 {
410  fes.Copy(spaces);
411
412  gamma = 0.00001;
413
414  // The mass matrix and preconditioner do not change every Newton cycle, so
415  // we only need to define them once
416  HypreBoomerAMG *mass_prec_amg = new HypreBoomerAMG();
417  mass_prec_amg->SetPrintLevel(0);
418
419  mass_prec = mass_prec_amg;
420
421  CGSolver *mass_pcg_iter = new CGSolver(spaces[0]->GetComm());
422  mass_pcg_iter->SetRelTol(1e-12);
423  mass_pcg_iter->SetAbsTol(1e-12);
424  mass_pcg_iter->SetMaxIter(200);
425  mass_pcg_iter->SetPrintLevel(0);
426  mass_pcg_iter->SetPreconditioner(*mass_prec);
427  mass_pcg_iter->SetOperator(*pressure_mass);
428  mass_pcg_iter->iterative_mode = false;
429
430  mass_pcg = mass_pcg_iter;
431
432  // The stiffness matrix does change every Newton cycle, so we will define it
433  // during SetOperator
434  stiff_pcg = NULL;
435  stiff_prec = NULL;
436 }
437
438 void JacobianPreconditioner::Mult(const Vector &k, Vector &y) const
439 {
440  // Extract the blocks from the input and output vectors
441  Vector disp_in(k.GetData() + block_trueOffsets[0],
442  block_trueOffsets[1]-block_trueOffsets[0]);
443  Vector pres_in(k.GetData() + block_trueOffsets[1],
444  block_trueOffsets[2]-block_trueOffsets[1]);
445
446  Vector disp_out(y.GetData() + block_trueOffsets[0],
447  block_trueOffsets[1]-block_trueOffsets[0]);
448  Vector pres_out(y.GetData() + block_trueOffsets[1],
449  block_trueOffsets[2]-block_trueOffsets[1]);
450
451  Vector temp(block_trueOffsets[1]-block_trueOffsets[0]);
452  Vector temp2(block_trueOffsets[1]-block_trueOffsets[0]);
453
454  // Perform the block elimination for the preconditioner
455  mass_pcg->Mult(pres_in, pres_out);
456  pres_out *= -gamma;
457
458  jacobian->GetBlock(0,1).Mult(pres_out, temp);
459  subtract(disp_in, temp, temp2);
460
461  stiff_pcg->Mult(temp2, disp_out);
462 }
463
464 void JacobianPreconditioner::SetOperator(const Operator &op)
465 {
466  jacobian = (BlockOperator *) &op;
467
468  // Initialize the stiffness preconditioner and solver
469  if (stiff_prec == NULL)
470  {
471  HypreBoomerAMG *stiff_prec_amg = new HypreBoomerAMG();
472  stiff_prec_amg->SetPrintLevel(0);
473
474  if (!spaces[0]->GetParMesh()->Nonconforming())
475  {
476  stiff_prec_amg->SetElasticityOptions(spaces[0]);
477  }
478
479  stiff_prec = stiff_prec_amg;
480
481  GMRESSolver *stiff_pcg_iter = new GMRESSolver(spaces[0]->GetComm());
482  stiff_pcg_iter->SetRelTol(1e-8);
483  stiff_pcg_iter->SetAbsTol(1e-8);
484  stiff_pcg_iter->SetMaxIter(200);
485  stiff_pcg_iter->SetPrintLevel(0);
486  stiff_pcg_iter->SetPreconditioner(*stiff_prec);
487  stiff_pcg_iter->iterative_mode = false;
488
489  stiff_pcg = stiff_pcg_iter;
490  }
491
492  // At each Newton cycle, compute the new stiffness AMG preconditioner by
493  // updating the iterative solver which, in turn, updates its preconditioner
494  stiff_pcg->SetOperator(jacobian->GetBlock(0,0));
495 }
496
497 JacobianPreconditioner::~JacobianPreconditioner()
498 {
499  delete mass_pcg;
500  delete mass_prec;
501  delete stiff_prec;
502  delete stiff_pcg;
503 }
504
505
506 RubberOperator::RubberOperator(Array<ParFiniteElementSpace *> &fes,
507  Array<Array<int> *> &ess_bdr,
508  Array<int> &trueOffsets,
509  double rel_tol,
510  double abs_tol,
511  int iter,
512  Coefficient &c_mu)
513  : Operator(fes[0]->TrueVSize() + fes[1]->TrueVSize()),
514  newton_solver(fes[0]->GetComm()),
515  newton_monitor(fes[0]->GetComm(), "Newton", 1),
516  j_monitor(fes[0]->GetComm(), " GMRES", 3),
517  mu(c_mu), block_trueOffsets(trueOffsets)
518 {
519  Array<Vector *> rhs(2);
520  rhs = NULL; // Set all entries in the array
521
522  fes.Copy(spaces);
523
524  // Define the block nonlinear form
525  Hform = new ParBlockNonlinearForm(spaces);
526
527  // Add the incompressible neo-Hookean integrator
529
530  // Set the essential boundary conditions
531  Hform->SetEssentialBC(ess_bdr, rhs);
532
533  // Compute the pressure mass stiffness matrix
534  ParBilinearForm *a = new ParBilinearForm(spaces[1]);
535  ConstantCoefficient one(1.0);
536  OperatorHandle mass(Operator::Hypre_ParCSR);
538  a->Assemble();
539  a->Finalize();
540  a->ParallelAssemble(mass);
541  delete a;
542
543  mass.SetOperatorOwner(false);
544  pressure_mass = mass.Ptr();
545
546  // Initialize the Jacobian preconditioner
547  JacobianPreconditioner *jac_prec =
548  new JacobianPreconditioner(fes, *pressure_mass, block_trueOffsets);
549  j_prec = jac_prec;
550
551  // Set up the Jacobian solver
552  GMRESSolver *j_gmres = new GMRESSolver(spaces[0]->GetComm());
553  j_gmres->iterative_mode = false;
554  j_gmres->SetRelTol(1e-12);
555  j_gmres->SetAbsTol(1e-12);
556  j_gmres->SetMaxIter(300);
557  j_gmres->SetPrintLevel(-1);
558  j_gmres->SetMonitor(j_monitor);
559  j_gmres->SetPreconditioner(*j_prec);
560  j_solver = j_gmres;
561
562  // Set the newton solve parameters
563  newton_solver.iterative_mode = true;
564  newton_solver.SetSolver(*j_solver);
565  newton_solver.SetOperator(*this);
566  newton_solver.SetPrintLevel(-1);
567  newton_solver.SetMonitor(newton_monitor);
568  newton_solver.SetRelTol(rel_tol);
569  newton_solver.SetAbsTol(abs_tol);
570  newton_solver.SetMaxIter(iter);
571 }
572
573 // Solve the Newton system
574 void RubberOperator::Solve(Vector &xp) const
575 {
576  Vector zero;
577  newton_solver.Mult(zero, xp);
578  MFEM_VERIFY(newton_solver.GetConverged(),
579  "Newton Solver did not converge.");
580 }
581
582 // compute: y = H(x,p)
583 void RubberOperator::Mult(const Vector &k, Vector &y) const
584 {
585  Hform->Mult(k, y);
586 }
587
588 // Compute the Jacobian from the nonlinear form
589 Operator &RubberOperator::GetGradient(const Vector &xp) const
590 {
592 }
593
594 RubberOperator::~RubberOperator()
595 {
596  delete Hform;
597  delete pressure_mass;
598  delete j_solver;
599  delete j_prec;
600 }
601
602
603 // Inline visualization
604 void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
605  ParGridFunction *field, const char *field_name, bool init_vis)
606 {
607  if (!out)
608  {
609  return;
610  }
611
612  GridFunction *nodes = deformed_nodes;
613  int owns_nodes = 0;
614
615  mesh->SwapNodes(nodes, owns_nodes);
616
617  out << "parallel " << mesh->GetNRanks() << " " << mesh->GetMyRank() << "\n";
618  out << "solution\n" << *mesh << *field;
619
620  mesh->SwapNodes(nodes, owns_nodes);
621
622  if (init_vis)
623  {
624  out << "window_size 800 800\n";
625  out << "window_title '" << field_name << "'\n";
626  if (mesh->SpaceDimension() == 2)
627  {
628  out << "view 0 0\n"; // view from top
629  out << "keys jlA\n"; // turn off perspective and light, +anti-aliasing
630  }
631  out << "keys cmA\n"; // show colorbar and mesh, +anti-aliasing
632  out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
633  }
634  out << flush;
635 }
636
637 void ReferenceConfiguration(const Vector &x, Vector &y)
638 {
639  // Set the reference, stress free, configuration
640  y = x;
641 }
642
643 void InitialDeformation(const Vector &x, Vector &y)
644 {
645  // Set the initial configuration. Having this different from the reference
646  // configuration can help convergence
647  y = x;
648  y[1] = x[1] + 0.25*x[0];
649 }
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
Definition: ex10.cpp:379
void InitialDeformation(const Vector &x, Vector &y)
Definition: ex10.cpp:591
Definition: solvers.hpp:258
void ReferenceConfiguration(const Vector &x, Vector &y)
Definition: ex19.cpp:577
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:78
int TrueVSize() const
Obsolete, kept for backward compatibility.
Definition: pfespace.hpp:382
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:97
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Definition: mesh.cpp:6641
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
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:831
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:819
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:474
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:638
int main(int argc, char *argv[])
Definition: ex1.cpp:66
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:169
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
int GetNRanks() const
Definition: pmesh.hpp:237
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1079
A simple convenience class that calls MPI_Init() at construction and MPI_Finalize() at destruction...
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:70
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:97
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
constexpr char vishost[]
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:314
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:8382
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1119
constexpr int visport
void SetMaxIter(int max_it)
Definition: solvers.hpp:98
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
HYPRE_Int GlobalTrueVSize() const
Definition: pfespace.hpp:251
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:406
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4169
void SetElasticityOptions(ParFiniteElementSpace *fespace)
Definition: hypre.cpp:3561
int Dimension() const
Definition: mesh.hpp:788
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:434
A general vector function coefficient.
int SpaceDimension() const
Definition: mesh.hpp:789
Abstract base class for an iterative solver monitor.
Definition: solvers.hpp:34
void SetAbsTol(double atol)
Definition: solvers.hpp:97
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:201
int GetMyRank() const
Definition: pmesh.hpp:238
void SetRelTol(double rtol)
Definition: solvers.hpp:96
int WorldRank() const
Return MPI_COMM_WORLD&#39;s rank.
MPI_Comm GetComm() const
Definition: pmesh.hpp:236
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:408
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
GMRES method.
Definition: solvers.hpp:290
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:112
double a
Definition: lissajous.cpp:41
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
double mu
Definition: ex25.cpp:135
int dim
Definition: ex24.cpp:53
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:304
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:272
Vector data type.
Definition: vector.hpp:51
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:91
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:159
Base class for solvers.
Definition: operator.hpp:634
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
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:87
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:145