MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
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
39using namespace std;
40using namespace mfem;
41
42class GeneralResidualMonitor : public IterativeSolverMonitor
43{
44public:
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, real_t norm, const Vector &r, bool final);
66
67private:
68 const std::string prefix;
69 int print_level;
70 mutable real_t norm0;
71};
72
73void GeneralResidualMonitor::MonitorResidual(int it, real_t 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.
107class JacobianPreconditioner : public Solver
108{
109protected:
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 real_t 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
129public:
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.
143class RubberOperator : public Operator
144{
145protected:
146 // Finite element spaces
148
149 // Block nonlinear form
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
172public:
173 RubberOperator(Array<ParFiniteElementSpace *> &fes, Array<Array<int> *>&ess_bdr,
174 Array<int> &block_trueOffsets, real_t rel_tol, real_t 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
188void 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
194void ReferenceConfiguration(const Vector &x, Vector &y);
195void InitialDeformation(const Vector &x, Vector &y);
196
197
198int 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 MFEM_SKIP_RETURN_VALUE;
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 real_t newton_rel_tol = 1e-4;
218 real_t newton_abs_tol = 1e-6;
219 int newton_iter = 500;
220 real_t mu = 1.0;
221
222 OptionsParser args(argc, argv);
224 "Mesh file to use.");
226 "Number of times to refine the mesh uniformly in serial.");
228 "Number of times to refine the mesh uniformly in parallel.");
230 "Order (degree) of the finite elements.");
232 "--no-visualization",
233 "Enable or disable GLVis visualization.");
235 "Relative tolerance for the Newton solve.");
237 "Absolute tolerance for the Newton solve.");
239 "Maximum iterations for the Newton solve.");
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
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.
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
412JacobianPreconditioner::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
446void 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
479void 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
515JacobianPreconditioner::~JacobianPreconditioner()
516{
517 delete mass_pcg;
518 delete mass_prec;
519 delete stiff_prec;
520 delete stiff_pcg;
521}
522
523
524RubberOperator::RubberOperator(Array<ParFiniteElementSpace *> &fes,
525 Array<Array<int> *> &ess_bdr,
526 Array<int> &trueOffsets,
527 real_t rel_tol,
528 real_t 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
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);
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
592void 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)
601void RubberOperator::Mult(const Vector &k, Vector &y) const
602{
603 Hform->Mult(k, y);
604}
605
606// Compute the Jacobian from the nonlinear form
608{
610}
611
612RubberOperator::~RubberOperator()
613{
614 delete Hform;
615 delete pressure_mass;
616 delete j_solver;
617 delete j_prec;
618}
619
620
621// Inline visualization
622void 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
660{
661 // Set the reference, stress free, configuration
662 y = x;
663}
664
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}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition array.cpp:103
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition array.hpp:874
virtual Operator & GetGradient(const Vector &x) const
virtual void Mult(const Vector &x, Vector &y) const
A class to handle Block systems in a matrix-free implementation.
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
A class to handle Vectors in a block fashion.
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition solvers.hpp:513
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:526
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
A coefficient that is constant across space and time.
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
GMRES method.
Definition solvers.hpp:547
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetElasticityOptions(ParFiniteElementSpace *fespace, bool interp_refine=true)
Definition hypre.cpp:5238
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Abstract base class for an iterative solver monitor.
Definition solvers.hpp:37
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
void SetMonitor(IterativeSolverMonitor &m)
Set the iterative solver monitor.
Definition solvers.hpp:296
bool GetConverged() const
Returns true if the last call to Mult() converged successfully.
Definition solvers.hpp:262
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Swap the internal node GridFunction pointer and ownership flag members with the given ones.
Definition mesh.cpp:9022
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
Newton's method for solving F(x)=b for a given operator F.
Definition solvers.hpp:667
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition solvers.cpp:1821
Pointer to an Operator of a specified type.
Definition handle.hpp:34
Abstract operator.
Definition operator.hpp:25
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
@ Hypre_ParCSR
ID for class HypreParMatrix.
Definition operator.hpp:287
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Class for parallel bilinear form.
A class representing a general parallel block nonlinear operator defined on the Cartesian product of ...
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
int TrueVSize() const
Obsolete, kept for backward compatibility.
Definition pfespace.hpp:436
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Save(std::ostream &out) const override
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void Distribute(const Vector *tv)
Class for parallel meshes.
Definition pmesh.hpp:34
MPI_Comm GetComm() const
Definition pmesh.hpp:402
int GetMyRank() const
Definition pmesh.hpp:404
int GetNRanks() const
Definition pmesh.hpp:403
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Base class for solvers.
Definition operator.hpp:683
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition operator.hpp:686
virtual void SetOperator(const Operator &op)=0
Set/update the solver for the given operator.
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
void SyncAliasMemory(const Vector &v) const
Update the alias memory location of the vector to match v.
Definition vector.hpp:259
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:602
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
void ReferenceConfiguration(const Vector &x, Vector &y)
Definition ex19p.cpp:659
void visualize(ostream &os, ParMesh *mesh, ParGridFunction *deformed_nodes, ParGridFunction *field, const char *field_name=NULL, bool init_vis=false)
Definition ex19p.cpp:622
void InitialDeformation(const Vector &x, Vector &y)
Definition ex19p.cpp:665
int dim
Definition ex24.cpp:53
real_t mu
Definition ex25.cpp:140
int main()
HYPRE_Int HYPRE_BigInt
real_t a
Definition lissajous.cpp:41
const int visport
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
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition vector.cpp:472
float real_t
Definition config.hpp:43
const char vishost[]