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