MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex10p.cpp
Go to the documentation of this file.
1// MFEM Example 10 - Parallel Version
2// PETSc Modification
3//
4// Compile with: make ex10p
5//
6// Sample runs:
7// mpirun -np 4 ex10p -m ../../data/beam-quad.mesh --petscopts rc_ex10p -s 3 -rs 2 -dt 3
8// mpirun -np 4 ex10p -m ../../data/beam-quad-amr.mesh --petscopts rc_ex10p -s 3 -rs 2 -dt 3
9//
10// Description: This examples solves a time dependent nonlinear elasticity
11// problem of the form dv/dt = H(x) + S v, dx/dt = v, where H is a
12// hyperelastic model and S is a viscosity operator of Laplacian
13// type. The geometry of the domain is assumed to be as follows:
14//
15// +---------------------+
16// boundary --->| |
17// attribute 1 | |
18// (fixed) +---------------------+
19//
20// The example demonstrates the use of nonlinear operators (the
21// class HyperelasticOperator defining H(x)), as well as their
22// implicit time integration using a Newton method for solving an
23// associated reduced backward-Euler type nonlinear equation
24// (class ReducedSystemOperator). Each Newton step requires the
25// inversion of a Jacobian matrix, which is done through a
26// (preconditioned) inner solver. Note that implementing the
27// method HyperelasticOperator::ImplicitSolve is the only
28// requirement for high-order implicit (SDIRK) time integration.
29// If using PETSc to solve the nonlinear problem, use the option
30// files provided (see rc_ex10p, rc_ex10p_mf, rc_ex10p_mfop) that
31// customize the Newton-Krylov method.
32// When option --jfnk is used, PETSc will use a Jacobian-free
33// Newton-Krylov method, using a user-defined preconditioner
34// constructed with the PetscPreconditionerFactory class.
35//
36// We recommend viewing examples 2 and 9 before viewing this
37// example.
38
39#include "mfem.hpp"
40#include <memory>
41#include <iostream>
42#include <fstream>
43
44#ifndef MFEM_USE_PETSC
45#error This example requires that MFEM is built with MFEM_USE_PETSC=YES
46#endif
47
48using namespace std;
49using namespace mfem;
50
51class ReducedSystemOperator;
52
53/** After spatial discretization, the hyperelastic model can be written as a
54 * system of ODEs:
55 * dv/dt = -M^{-1}*(H(x) + S*v)
56 * dx/dt = v,
57 * where x is the vector representing the deformation, v is the velocity field,
58 * M is the mass matrix, S is the viscosity matrix, and H(x) is the nonlinear
59 * hyperelastic operator.
60 *
61 * Class HyperelasticOperator represents the right-hand side of the above
62 * system of ODEs. */
63class HyperelasticOperator : public TimeDependentOperator
64{
65protected:
66 ParFiniteElementSpace &fespace;
67 Array<int> ess_tdof_list;
68
69 ParBilinearForm M, S;
71 real_t viscosity;
72 HyperelasticModel *model;
73
74 HypreParMatrix *Mmat; // Mass matrix from ParallelAssemble()
75 CGSolver M_solver; // Krylov solver for inverting the mass matrix M
76 HypreSmoother M_prec; // Preconditioner for the mass matrix M
77
78 /** Nonlinear operator defining the reduced backward Euler equation for the
79 velocity. Used in the implementation of method ImplicitSolve. */
80 ReducedSystemOperator *reduced_oper;
81
82 /// Newton solver for the reduced backward Euler equation
83 NewtonSolver newton_solver;
84
85 /// Newton solver for the reduced backward Euler equation (PETSc SNES)
86 PetscNonlinearSolver* pnewton_solver;
87
88 /// Solver for the Jacobian solve in the Newton method
89 Solver *J_solver;
90 /// Preconditioner for the Jacobian solve in the Newton method
91 Solver *J_prec;
92 /// Preconditioner factory for JFNK
94
95 mutable Vector z; // auxiliary vector
96
97public:
98 HyperelasticOperator(ParFiniteElementSpace &f, Array<int> &ess_bdr,
99 real_t visc, real_t mu, real_t K,
100 bool use_petsc, bool petsc_use_jfnk);
101
102 /// Compute the right-hand side of the ODE system.
103 virtual void Mult(const Vector &vx, Vector &dvx_dt) const;
104 /** Solve the Backward-Euler equation: k = f(x + dt*k, t), for the unknown k.
105 This is the only requirement for high-order SDIRK implicit integration.*/
106 virtual void ImplicitSolve(const real_t dt, const Vector &x, Vector &k);
107
108 real_t ElasticEnergy(const ParGridFunction &x) const;
109 real_t KineticEnergy(const ParGridFunction &v) const;
110 void GetElasticEnergyDensity(const ParGridFunction &x,
111 ParGridFunction &w) const;
112
113 virtual ~HyperelasticOperator();
114};
115
116/** Nonlinear operator of the form:
117 k --> (M + dt*S)*k + H(x + dt*v + dt^2*k) + S*v,
118 where M and S are given BilinearForms, H is a given NonlinearForm, v and x
119 are given vectors, and dt is a scalar. */
120class ReducedSystemOperator : public Operator
121{
122private:
123 ParBilinearForm *M, *S;
125 mutable HypreParMatrix *Jacobian;
126 real_t dt;
127 const Vector *v, *x;
128 mutable Vector w, z;
129 const Array<int> &ess_tdof_list;
130
131public:
132 ReducedSystemOperator(ParBilinearForm *M_, ParBilinearForm *S_,
133 ParNonlinearForm *H_, const Array<int> &ess_tdof_list);
134
135 /// Set current dt, v, x values - needed to compute action and Jacobian.
136 void SetParameters(real_t dt_, const Vector *v_, const Vector *x_);
137
138 /// Compute y = H(x + dt (v + dt k)) + M k + S (v + dt k).
139 virtual void Mult(const Vector &k, Vector &y) const;
140
141 /// Compute J = M + dt S + dt^2 grad_H(x + dt (v + dt k)).
142 virtual Operator &GetGradient(const Vector &k) const;
143
144 virtual ~ReducedSystemOperator();
145
146};
147
148/** Auxiliary class to provide preconditioners for matrix-free methods */
149class PreconditionerFactory : public PetscPreconditionerFactory
150{
151private:
152 // const ReducedSystemOperator& op; // unused for now (generates warning)
153
154public:
155 PreconditionerFactory(const ReducedSystemOperator& op_, const string& name_)
156 : PetscPreconditionerFactory(name_) /* , op(op_) */ {}
157 virtual mfem::Solver* NewPreconditioner(const mfem::OperatorHandle&);
158 virtual ~PreconditionerFactory() {}
159};
160
161/** Function representing the elastic energy density for the given hyperelastic
162 model+deformation. Used in HyperelasticOperator::GetElasticEnergyDensity. */
163class ElasticEnergyCoefficient : public Coefficient
164{
165private:
166 HyperelasticModel &model;
167 const ParGridFunction &x;
168 DenseMatrix J;
169
170public:
171 ElasticEnergyCoefficient(HyperelasticModel &m, const ParGridFunction &x_)
172 : model(m), x(x_) { }
173 virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip);
174 virtual ~ElasticEnergyCoefficient() { }
175};
176
177void InitialDeformation(const Vector &x, Vector &y);
178
179void InitialVelocity(const Vector &x, Vector &v);
180
181void visualize(ostream &os, ParMesh *mesh, ParGridFunction *deformed_nodes,
182 ParGridFunction *field, const char *field_name = NULL,
183 bool init_vis = false);
184
185
186int main(int argc, char *argv[])
187{
188 // 1. Initialize MPI and HYPRE.
189 Mpi::Init(argc, argv);
190 int myid = Mpi::WorldRank();
191 Hypre::Init();
192
193 // 2. Parse command-line options.
194 const char *mesh_file = "../../data/beam-quad.mesh";
195 int ser_ref_levels = 2;
196 int par_ref_levels = 0;
197 int order = 2;
198 int ode_solver_type = 3;
199 real_t t_final = 300.0;
200 real_t dt = 3.0;
201 real_t visc = 1e-2;
202 real_t mu = 0.25;
203 real_t K = 5.0;
204 bool visualization = true;
205 int vis_steps = 1;
206 bool use_petsc = true;
207 const char *petscrc_file = "";
208 bool petsc_use_jfnk = false;
209
210 OptionsParser args(argc, argv);
212 "Mesh file to use.");
214 "Number of times to refine the mesh uniformly in serial.");
216 "Number of times to refine the mesh uniformly in parallel.");
218 "Order (degree) of the finite elements.");
220 "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
221 " 11 - Forward Euler, 12 - RK2,\n\t"
222 " 13 - RK3 SSP, 14 - RK4.");
224 "Final time; start time is 0.");
226 "Time step.");
228 "Viscosity coefficient.");
230 "Shear modulus in the Neo-Hookean hyperelastic model.");
232 "Bulk modulus in the Neo-Hookean hyperelastic model.");
234 "--no-visualization",
235 "Enable or disable GLVis visualization.");
237 "Visualize every n-th timestep.");
239 "--no-petsc",
240 "Use or not PETSc to solve the nonlinear system.");
242 "PetscOptions file to use.");
244 "--no-jfnk",
245 "Use JFNK with user-defined preconditioner factory.");
246 args.Parse();
247 if (!args.Good())
248 {
249 if (myid == 0)
250 {
251 args.PrintUsage(cout);
252 }
253 return 1;
254 }
255 if (myid == 0)
256 {
257 args.PrintOptions(cout);
258 }
259
260 // 2b. We initialize PETSc
261 if (use_petsc)
262 {
263 MFEMInitializePetsc(NULL,NULL,petscrc_file,NULL);
264 }
265
266 // 3. Read the serial mesh from the given mesh file on all processors. We can
267 // handle triangular, quadrilateral, tetrahedral and hexahedral meshes
268 // with the same code.
269 Mesh *mesh = new Mesh(mesh_file, 1, 1);
270 int dim = mesh->Dimension();
271
272 // 4. Define the ODE solver used for time integration. Several implicit
273 // singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as
274 // explicit Runge-Kutta methods are available.
275 ODESolver *ode_solver;
276 switch (ode_solver_type)
277 {
278 // Implicit L-stable methods
279 case 1: ode_solver = new BackwardEulerSolver; break;
280 case 2: ode_solver = new SDIRK23Solver(2); break;
281 case 3: ode_solver = new SDIRK33Solver; break;
282 // Explicit methods
283 case 11: ode_solver = new ForwardEulerSolver; break;
284 case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
285 case 13: ode_solver = new RK3SSPSolver; break;
286 case 14: ode_solver = new RK4Solver; break;
287 // Implicit A-stable methods (not L-stable)
288 case 22: ode_solver = new ImplicitMidpointSolver; break;
289 case 23: ode_solver = new SDIRK23Solver; break;
290 case 24: ode_solver = new SDIRK34Solver; break;
291 default:
292 if (myid == 0)
293 {
294 cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
295 }
296 return 3;
297 }
298
299 // 5. Refine the mesh in serial to increase the resolution. In this example
300 // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
301 // a command-line parameter.
302 for (int lev = 0; lev < ser_ref_levels; lev++)
303 {
304 mesh->UniformRefinement();
305 }
306
307 // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
308 // this mesh further in parallel to increase the resolution. Once the
309 // parallel mesh is defined, the serial mesh can be deleted.
310 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
311 delete mesh;
312 for (int lev = 0; lev < par_ref_levels; lev++)
313 {
314 pmesh->UniformRefinement();
315 }
316
317 // 7. Define the parallel vector finite element spaces representing the mesh
318 // deformation x_gf, the velocity v_gf, and the initial configuration,
319 // x_ref. Define also the elastic energy density, w_gf, which is in a
320 // discontinuous higher-order space. Since x and v are integrated in time
321 // as a system, we group them together in block vector vx, on the unique
322 // parallel degrees of freedom, with offsets given by array true_offset.
323 H1_FECollection fe_coll(order, dim);
324 ParFiniteElementSpace fespace(pmesh, &fe_coll, dim);
325
326 HYPRE_BigInt glob_size = fespace.GlobalTrueVSize();
327 if (myid == 0)
328 {
329 cout << "Number of velocity/deformation unknowns: " << glob_size << endl;
330 }
331 int true_size = fespace.TrueVSize();
332 Array<int> true_offset(3);
333 true_offset[0] = 0;
334 true_offset[1] = true_size;
335 true_offset[2] = 2*true_size;
336
337 BlockVector vx(true_offset);
338 ParGridFunction v_gf, x_gf;
339 v_gf.MakeTRef(&fespace, vx, true_offset[0]);
340 x_gf.MakeTRef(&fespace, vx, true_offset[1]);
341
342 ParGridFunction x_ref(&fespace);
343 pmesh->GetNodes(x_ref);
344
345 L2_FECollection w_fec(order + 1, dim);
346 ParFiniteElementSpace w_fespace(pmesh, &w_fec);
347 ParGridFunction w_gf(&w_fespace);
348
349 // 8. Set the initial conditions for v_gf, x_gf and vx, and define the
350 // boundary conditions on a beam-like mesh (see description above).
352 v_gf.ProjectCoefficient(velo);
353 v_gf.SetTrueVector();
355 x_gf.ProjectCoefficient(deform);
356 x_gf.SetTrueVector();
357
359
360 Array<int> ess_bdr(fespace.GetMesh()->bdr_attributes.Max());
361 ess_bdr = 0;
362 ess_bdr[0] = 1; // boundary attribute 1 (index 0) is fixed
363
364 // 9. Initialize the hyperelastic operator, the GLVis visualization and print
365 // the initial energies.
366 HyperelasticOperator *oper = new HyperelasticOperator(fespace, ess_bdr, visc,
367 mu, K, use_petsc,
368 petsc_use_jfnk);
369
370 socketstream vis_v, vis_w;
371 if (visualization)
372 {
373 char vishost[] = "localhost";
374 int visport = 19916;
375 vis_v.open(vishost, visport);
376 vis_v.precision(8);
377 visualize(vis_v, pmesh, &x_gf, &v_gf, "Velocity", true);
378 // Make sure all ranks have sent their 'v' solution before initiating
379 // another set of GLVis connections (one from each rank):
380 MPI_Barrier(pmesh->GetComm());
381 vis_w.open(vishost, visport);
382 if (vis_w)
383 {
384 oper->GetElasticEnergyDensity(x_gf, w_gf);
385 vis_w.precision(8);
386 visualize(vis_w, pmesh, &x_gf, &w_gf, "Elastic energy density", true);
387 }
388 }
389
390 real_t ee0 = oper->ElasticEnergy(x_gf);
391 real_t ke0 = oper->KineticEnergy(v_gf);
392 if (myid == 0)
393 {
394 cout << "initial elastic energy (EE) = " << ee0 << endl;
395 cout << "initial kinetic energy (KE) = " << ke0 << endl;
396 cout << "initial total energy (TE) = " << (ee0 + ke0) << endl;
397 }
398
399 real_t t = 0.0;
400 oper->SetTime(t);
401 ode_solver->Init(*oper);
402
403 // 10. Perform time-integration
404 // (looping over the time iterations, ti, with a time-step dt).
405 bool last_step = false;
406 for (int ti = 1; !last_step; ti++)
407 {
408 real_t dt_real = min(dt, t_final - t);
409
410 ode_solver->Step(vx, t, dt_real);
411
412 last_step = (t >= t_final - 1e-8*dt);
413
414 if (last_step || (ti % vis_steps) == 0)
415 {
417
418 real_t ee = oper->ElasticEnergy(x_gf);
419 real_t ke = oper->KineticEnergy(v_gf);
420
421 if (myid == 0)
422 {
423 cout << "step " << ti << ", t = " << t << ", EE = " << ee
424 << ", KE = " << ke << ", Î”TE = " << (ee+ke)-(ee0+ke0) << endl;
425 }
426
427 if (visualization)
428 {
429 visualize(vis_v, pmesh, &x_gf, &v_gf);
430 if (vis_w)
431 {
432 oper->GetElasticEnergyDensity(x_gf, w_gf);
433 visualize(vis_w, pmesh, &x_gf, &w_gf);
434 }
435 }
436 }
437 }
438
439 // 11. Save the displaced mesh, the velocity and elastic energy.
440 {
442 GridFunction *nodes = &x_gf;
443 int owns_nodes = 0;
444 pmesh->SwapNodes(nodes, owns_nodes);
445
446 ostringstream mesh_name, velo_name, ee_name;
447 mesh_name << "deformed." << setfill('0') << setw(6) << myid;
448 velo_name << "velocity." << setfill('0') << setw(6) << myid;
449 ee_name << "elastic_energy." << setfill('0') << setw(6) << myid;
450
451 ofstream mesh_ofs(mesh_name.str().c_str());
452 mesh_ofs.precision(8);
453 pmesh->Print(mesh_ofs);
454 pmesh->SwapNodes(nodes, owns_nodes);
455 ofstream velo_ofs(velo_name.str().c_str());
456 velo_ofs.precision(8);
457 v_gf.Save(velo_ofs);
458 ofstream ee_ofs(ee_name.str().c_str());
459 ee_ofs.precision(8);
460 oper->GetElasticEnergyDensity(x_gf, w_gf);
461 w_gf.Save(ee_ofs);
462 }
463
464 // 12. Free the used memory.
465 delete ode_solver;
466 delete pmesh;
467 delete oper;
468
469 // We finalize PETSc
470 if (use_petsc) { MFEMFinalizePetsc(); }
471
472 return 0;
473}
474
475void visualize(ostream &os, ParMesh *mesh, ParGridFunction *deformed_nodes,
476 ParGridFunction *field, const char *field_name, bool init_vis)
477{
478 if (!os)
479 {
480 return;
481 }
482
483 GridFunction *nodes = deformed_nodes;
484 int owns_nodes = 0;
485
486 mesh->SwapNodes(nodes, owns_nodes);
487
488 os << "parallel " << mesh->GetNRanks() << " " << mesh->GetMyRank() << "\n";
489 os << "solution\n" << *mesh << *field;
490
491 mesh->SwapNodes(nodes, owns_nodes);
492
493 if (init_vis)
494 {
495 os << "window_size 800 800\n";
496 os << "window_title '" << field_name << "'\n";
497 if (mesh->SpaceDimension() == 2)
498 {
499 os << "view 0 0\n"; // view from top
500 os << "keys jl\n"; // turn off perspective and light
501 }
502 os << "keys cm\n"; // show colorbar and mesh
503 os << "autoscale value\n"; // update value-range; keep mesh-extents fixed
504 os << "pause\n";
505 }
506 os << flush;
507}
508
509
510ReducedSystemOperator::ReducedSystemOperator(
512 const Array<int> &ess_tdof_list_)
513 : Operator(M_->ParFESpace()->TrueVSize()), M(M_), S(S_), H(H_),
514 Jacobian(NULL), dt(0.0), v(NULL), x(NULL), w(height), z(height),
515 ess_tdof_list(ess_tdof_list_)
516{ }
517
518void ReducedSystemOperator::SetParameters(real_t dt_, const Vector *v_,
519 const Vector *x_)
520{
521 dt = dt_; v = v_; x = x_;
522}
523
524void ReducedSystemOperator::Mult(const Vector &k, Vector &y) const
525{
526 // compute: y = H(x + dt*(v + dt*k)) + M*k + S*(v + dt*k)
529 H->Mult(z, y);
532 y.SetSubVector(ess_tdof_list, 0.0);
533}
534
536{
537 delete Jacobian;
538 SparseMatrix *localJ = Add(1.0, M->SpMat(), dt, S->SpMat());
542 // if we are using PETSc, the HypreParCSR Jacobian will be converted to
543 // PETSc's AIJ on the fly
544 Jacobian = M->ParallelAssemble(localJ);
545 delete localJ;
546 HypreParMatrix *Je = Jacobian->EliminateRowsCols(ess_tdof_list);
547 delete Je;
548 return *Jacobian;
549}
550
551ReducedSystemOperator::~ReducedSystemOperator()
552{
553 delete Jacobian;
554}
555
556
557HyperelasticOperator::HyperelasticOperator(ParFiniteElementSpace &f,
558 Array<int> &ess_bdr, real_t visc,
559 real_t mu, real_t K, bool use_petsc,
560 bool use_petsc_factory)
561 : TimeDependentOperator(2*f.TrueVSize(), static_cast<real_t>(0.0)), fespace(f),
562 M(&fespace), S(&fespace), H(&fespace),
563 viscosity(visc), M_solver(f.GetComm()),
564 newton_solver(f.GetComm()), pnewton_solver(NULL), z(height/2)
565{
566#if defined(MFEM_USE_DOUBLE)
567 const real_t rel_tol = 1e-8;
568 const real_t newton_abs_tol = 0.0;
569#elif defined(MFEM_USE_SINGLE)
570 const real_t rel_tol = 1e-3;
571 const real_t newton_abs_tol = 1e-4;
572#endif
573 const int skip_zero_entries = 0;
574
575 const real_t ref_density = 1.0; // density in the reference configuration
576 ConstantCoefficient rho0(ref_density);
578 M.Assemble(skip_zero_entries);
579 M.Finalize(skip_zero_entries);
580 Mmat = M.ParallelAssemble();
581 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
582 HypreParMatrix *Me = Mmat->EliminateRowsCols(ess_tdof_list);
583 delete Me;
584
585 M_solver.iterative_mode = false;
586 M_solver.SetRelTol(rel_tol);
587 M_solver.SetAbsTol(0.0);
588 M_solver.SetMaxIter(30);
589 M_solver.SetPrintLevel(0);
590 M_prec.SetType(HypreSmoother::Jacobi);
591 M_solver.SetPreconditioner(M_prec);
592 M_solver.SetOperator(*Mmat);
593
594 model = new NeoHookeanModel(mu, K);
596 H.SetEssentialTrueDofs(ess_tdof_list);
597
598 ConstantCoefficient visc_coeff(viscosity);
600 S.Assemble(skip_zero_entries);
601 S.Finalize(skip_zero_entries);
602
603 reduced_oper = new ReducedSystemOperator(&M, &S, &H, ess_tdof_list);
604 if (!use_petsc)
605 {
606 HypreSmoother *J_hypreSmoother = new HypreSmoother;
607 J_hypreSmoother->SetType(HypreSmoother::l1Jacobi);
608 J_hypreSmoother->SetPositiveDiagonal(true);
609 J_prec = J_hypreSmoother;
610
611 MINRESSolver *J_minres = new MINRESSolver(f.GetComm());
612 J_minres->SetRelTol(rel_tol);
613 J_minres->SetAbsTol(0.0);
614 J_minres->SetMaxIter(300);
615 J_minres->SetPrintLevel(-1);
616 J_minres->SetPreconditioner(*J_prec);
617 J_solver = J_minres;
618
619 J_factory = NULL;
620
621 newton_solver.iterative_mode = false;
622 newton_solver.SetSolver(*J_solver);
623 newton_solver.SetOperator(*reduced_oper);
624 newton_solver.SetPrintLevel(1); // print Newton iterations
625 newton_solver.SetRelTol(rel_tol);
626 newton_solver.SetAbsTol(newton_abs_tol);
627 newton_solver.SetMaxIter(10);
628 }
629 else
630 {
631 // if using PETSc, we create the same solver (Newton + MINRES + Jacobi)
632 // by command line options (see rc_ex10p)
633 J_solver = NULL;
634 J_prec = NULL;
635 J_factory = NULL;
636 pnewton_solver = new PetscNonlinearSolver(f.GetComm(),
637 *reduced_oper);
638
639 // we can setup a factory to construct a "physics-based" preconditioner
640 if (use_petsc_factory)
641 {
642 J_factory = new PreconditionerFactory(*reduced_oper, "JFNK preconditioner");
643 pnewton_solver->SetPreconditionerFactory(J_factory);
644 }
645 pnewton_solver->SetPrintLevel(1); // print Newton iterations
646 pnewton_solver->SetRelTol(rel_tol);
647 pnewton_solver->SetAbsTol(newton_abs_tol);
648 pnewton_solver->SetMaxIter(10);
649 }
650}
651
652void HyperelasticOperator::Mult(const Vector &vx, Vector &dvx_dt) const
653{
654 // Create views to the sub-vectors v, x of vx, and dv_dt, dx_dt of dvx_dt
655 int sc = height/2;
656 Vector v(vx.GetData() + 0, sc);
657 Vector x(vx.GetData() + sc, sc);
658 Vector dv_dt(dvx_dt.GetData() + 0, sc);
659 Vector dx_dt(dvx_dt.GetData() + sc, sc);
660
661 H.Mult(x, z);
662 if (viscosity != 0.0)
663 {
665 z.SetSubVector(ess_tdof_list, 0.0);
666 }
667 z.Neg(); // z = -z
668 M_solver.Mult(z, dv_dt);
669
670 dx_dt = v;
671}
672
673void HyperelasticOperator::ImplicitSolve(const real_t dt,
674 const Vector &vx, Vector &dvx_dt)
675{
676 int sc = height/2;
677 Vector v(vx.GetData() + 0, sc);
678 Vector x(vx.GetData() + sc, sc);
679 Vector dv_dt(dvx_dt.GetData() + 0, sc);
680 Vector dx_dt(dvx_dt.GetData() + sc, sc);
681
682 // By eliminating kx from the coupled system:
683 // kv = -M^{-1}*[H(x + dt*kx) + S*(v + dt*kv)]
684 // kx = v + dt*kv
685 // we reduce it to a nonlinear equation for kv, represented by the
686 // reduced_oper. This equation is solved with the newton_solver
687 // object (using J_solver and J_prec internally).
688 reduced_oper->SetParameters(dt, &v, &x);
689 Vector zero; // empty vector is interpreted as zero r.h.s. by NewtonSolver
690 if (!pnewton_solver)
691 {
692 newton_solver.Mult(zero, dv_dt);
693 MFEM_VERIFY(newton_solver.GetConverged(),
694 "Newton solver did not converge.");
695 }
696 else
697 {
698 pnewton_solver->Mult(zero, dv_dt);
699 MFEM_VERIFY(pnewton_solver->GetConverged(),
700 "Newton solver did not converge.");
701 }
703}
704
705real_t HyperelasticOperator::ElasticEnergy(const ParGridFunction &x) const
706{
707 return H.GetEnergy(x);
708}
709
710real_t HyperelasticOperator::KineticEnergy(const ParGridFunction &v) const
711{
712 real_t energy = 0.5*M.ParInnerProduct(v, v);
713 return energy;
714}
715
716void HyperelasticOperator::GetElasticEnergyDensity(
717 const ParGridFunction &x, ParGridFunction &w) const
718{
719 ElasticEnergyCoefficient w_coeff(*model, x);
720 w.ProjectCoefficient(w_coeff);
721}
722
723HyperelasticOperator::~HyperelasticOperator()
724{
725 delete J_solver;
726 delete J_prec;
727 delete J_factory;
728 delete reduced_oper;
729 delete model;
730 delete Mmat;
731 delete pnewton_solver;
732}
733
734// This method gets called every time we need a preconditioner "oh"
735// contains the PetscParMatrix that wraps the operator constructed in
737// In this example, we just return a customizable PetscPreconditioner
738// using that matrix. However, the OperatorHandle argument can be
739// ignored, and any "physics-based" solver can be constructed since we
741Solver* PreconditionerFactory::NewPreconditioner(const mfem::OperatorHandle& oh)
742{
743 PetscParMatrix *pP;
744 oh.Get(pP);
745 return new PetscPreconditioner(*pP,"jfnk_");
746}
747
748real_t ElasticEnergyCoefficient::Eval(ElementTransformation &T,
749 const IntegrationPoint &ip)
750{
751 model.SetTransformation(T);
753 // return model.EvalW(J); // in reference configuration
754 return model.EvalW(J)/J.Det(); // in deformed configuration
755}
756
757
759{
760 // set the initial configuration to be the same as the reference,
761 // stress free, configuration
762 y = x;
763}
764
765void InitialVelocity(const Vector &x, Vector &v)
766{
767 const int dim = x.Size();
768 const real_t s = 0.1/64.;
769
770 v = 0.0;
771 v(dim-1) = s*x(0)*x(0)*(8.0-x(0));
772 v(0) = -s*x(0)*x(0);
773}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
Backward Euler ODE solver. L-stable.
Definition ode.hpp:412
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
A class to handle Vectors in a block fashion.
Definition solvers.hpp:513
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:718
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.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
real_t Det() const
Definition densemat.cpp:535
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
The classical forward Euler method.
Definition ode.hpp:118
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
Compute the vector gradient with respect to the physical element variable.
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition gridfunc.hpp:144
void MakeTRef(FiniteElementSpace *f, real_t *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
Definition gridfunc.cpp:221
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition gridfunc.hpp:150
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Abstract class for hyperelastic models.
virtual real_t EvalW(const DenseMatrix &Jpt) const =0
Evaluate the strain energy density function, W = W(Jpt).
void SetTransformation(ElementTransformation &Ttr_)
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition hypre.cpp:2356
Parallel smoothers in hypre.
Definition hypre.hpp:1020
void SetPositiveDiagonal(bool pos=true)
After computing l1-norms, replace them with their absolute values.
Definition hypre.hpp:1136
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Definition hypre.cpp:3550
@ l1Jacobi
l1-scaled Jacobi
Definition hypre.hpp:1080
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Implicit midpoint method. A-stable, not L-stable.
Definition ode.hpp:425
Class for integration point with weight.
Definition intrules.hpp:35
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
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
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
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
MINRES method.
Definition solvers.hpp:628
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.hpp:640
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 GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
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
virtual void Mult(const Vector &x, Vector &y) const
Evaluate the action of the NonlinearForm.
virtual real_t GetEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition ode.hpp:24
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition ode.cpp:18
virtual void Step(Vector &x, real_t &t, real_t &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in].
Pointer to an Operator of a specified type.
Definition handle.hpp:34
void Get(OpType *&A) const
Return the Operator pointer statically cast to a given OpType.
Definition handle.hpp:114
Abstract operator.
Definition operator.hpp:25
int height
Dimension of the output / number of rows in the matrix.
Definition operator.hpp:27
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.
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
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
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
Parallel non-linear operator on the true dofs.
Abstract class for PETSc's nonlinear solvers.
Definition petsc.hpp:907
virtual void Mult(const Vector &b, Vector &x) const
Application of the solver.
Definition petsc.cpp:4117
Wrapper for PETSc's matrix class.
Definition petsc.hpp:319
Abstract class for PETSc's preconditioners.
Definition petsc.hpp:805
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition ode.hpp:151
The classical explicit forth-order Runge-Kutta method, RK4.
Definition ode.hpp:164
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
Data type sparse matrix.
Definition sparsemat.hpp:51
void Add(const int i, const int j, const real_t val)
Base abstract class for first order time dependent operators.
Definition operator.hpp:317
virtual void SetTime(const real_t t_)
Set the current time.
Definition operator.hpp:360
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
void Neg()
(*this) = -(*this)
Definition vector.cpp:300
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:604
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:227
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
int dim
Definition ex24.cpp:53
real_t mu
Definition ex25.cpp:140
int main()
HYPRE_Int HYPRE_BigInt
const int visport
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:316
void MFEMInitializePetsc()
Convenience functions to initialize/finalize PETSc.
Definition petsc.cpp:201
void MFEMFinalizePetsc()
Definition petsc.cpp:241
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
const char vishost[]
RefCoord t[3]
RefCoord s[3]
void visualize(ostream &os, ParMesh *mesh, ParGridFunction *deformed_nodes, ParGridFunction *field, const char *field_name=NULL, bool init_vis=false)
Definition ex10p.cpp:475
void InitialDeformation(const Vector &x, Vector &y)
Definition ex10p.cpp:758
void InitialVelocity(const Vector &x, Vector &v)
Definition ex10p.cpp:765