44 #ifndef MFEM_USE_PETSC
45 #error This example requires that MFEM is built with MFEM_USE_PETSC=YES
51 class ReducedSystemOperator;
80 ReducedSystemOperator *reduced_oper;
99 double visc,
double mu,
double K,
100 bool use_petsc,
bool petsc_use_jfnk);
106 virtual void ImplicitSolve(
const double dt,
const Vector &x,
Vector &k);
113 virtual ~HyperelasticOperator();
120 class ReducedSystemOperator :
public Operator
136 void SetParameters(
double dt_,
const Vector *v_,
const Vector *x_);
144 virtual ~ReducedSystemOperator();
155 PreconditionerFactory(
const ReducedSystemOperator& op_,
const string& name_)
158 virtual ~PreconditionerFactory() {}
163 class ElasticEnergyCoefficient :
public Coefficient
172 : model(m), x(x_) { }
174 virtual ~ElasticEnergyCoefficient() { }
183 bool init_vis =
false);
186 int main(
int argc,
char *argv[])
190 MPI_Init(&argc, &argv);
191 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
192 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
195 const char *mesh_file =
"../../data/beam-quad.mesh";
196 int ser_ref_levels = 2;
197 int par_ref_levels = 0;
199 int ode_solver_type = 3;
200 double t_final = 300.0;
205 bool visualization =
true;
207 bool use_petsc =
true;
208 const char *petscrc_file =
"";
209 bool petsc_use_jfnk =
false;
212 args.
AddOption(&mesh_file,
"-m",
"--mesh",
213 "Mesh file to use.");
214 args.
AddOption(&ser_ref_levels,
"-rs",
"--refine-serial",
215 "Number of times to refine the mesh uniformly in serial.");
216 args.
AddOption(&par_ref_levels,
"-rp",
"--refine-parallel",
217 "Number of times to refine the mesh uniformly in parallel.");
219 "Order (degree) of the finite elements.");
220 args.
AddOption(&ode_solver_type,
"-s",
"--ode-solver",
221 "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
222 " 11 - Forward Euler, 12 - RK2,\n\t"
223 " 13 - RK3 SSP, 14 - RK4.");
224 args.
AddOption(&t_final,
"-tf",
"--t-final",
225 "Final time; start time is 0.");
226 args.
AddOption(&dt,
"-dt",
"--time-step",
228 args.
AddOption(&visc,
"-v",
"--viscosity",
229 "Viscosity coefficient.");
230 args.
AddOption(&mu,
"-mu",
"--shear-modulus",
231 "Shear modulus in the Neo-Hookean hyperelastic model.");
232 args.
AddOption(&K,
"-K",
"--bulk-modulus",
233 "Bulk modulus in the Neo-Hookean hyperelastic model.");
234 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
235 "--no-visualization",
236 "Enable or disable GLVis visualization.");
237 args.
AddOption(&vis_steps,
"-vs",
"--visualization-steps",
238 "Visualize every n-th timestep.");
239 args.
AddOption(&use_petsc,
"-usepetsc",
"--usepetsc",
"-no-petsc",
241 "Use or not PETSc to solve the nonlinear system.");
242 args.
AddOption(&petscrc_file,
"-petscopts",
"--petscopts",
243 "PetscOptions file to use.");
244 args.
AddOption(&petsc_use_jfnk,
"-jfnk",
"--jfnk",
"-no-jfnk",
246 "Use JFNK with user-defined preconditioner factory.");
271 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
278 switch (ode_solver_type)
286 case 12: ode_solver =
new RK2Solver(0.5);
break;
288 case 14: ode_solver =
new RK4Solver;
break;
296 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
305 for (
int lev = 0; lev < ser_ref_levels; lev++)
315 for (
int lev = 0; lev < par_ref_levels; lev++)
332 cout <<
"Number of velocity/deformation unknowns: " << glob_size << endl;
337 true_offset[1] = true_size;
338 true_offset[2] = 2*true_size;
342 v_gf.
MakeTRef(&fespace, vx, true_offset[0]);
343 x_gf.
MakeTRef(&fespace, vx, true_offset[1]);
369 HyperelasticOperator *oper =
new HyperelasticOperator(fespace, ess_bdr, visc,
376 char vishost[] =
"localhost";
378 vis_v.
open(vishost, visport);
380 visualize(vis_v, pmesh, &x_gf, &v_gf,
"Velocity",
true);
384 vis_w.
open(vishost, visport);
387 oper->GetElasticEnergyDensity(x_gf, w_gf);
389 visualize(vis_w, pmesh, &x_gf, &w_gf,
"Elastic energy density",
true);
393 double ee0 = oper->ElasticEnergy(x_gf);
394 double ke0 = oper->KineticEnergy(v_gf);
397 cout <<
"initial elastic energy (EE) = " << ee0 << endl;
398 cout <<
"initial kinetic energy (KE) = " << ke0 << endl;
399 cout <<
"initial total energy (TE) = " << (ee0 + ke0) << endl;
404 ode_solver->
Init(*oper);
408 bool last_step =
false;
409 for (
int ti = 1; !last_step; ti++)
411 double dt_real = min(dt, t_final - t);
413 ode_solver->
Step(vx, t, dt_real);
415 last_step = (t >= t_final - 1e-8*dt);
417 if (last_step || (ti % vis_steps) == 0)
421 double ee = oper->ElasticEnergy(x_gf);
422 double ke = oper->KineticEnergy(v_gf);
426 cout <<
"step " << ti <<
", t = " << t <<
", EE = " << ee
427 <<
", KE = " << ke <<
", ΔTE = " << (ee+ke)-(ee0+ke0) << endl;
435 oper->GetElasticEnergyDensity(x_gf, w_gf);
449 ostringstream mesh_name, velo_name, ee_name;
450 mesh_name <<
"deformed." << setfill(
'0') << setw(6) << myid;
451 velo_name <<
"velocity." << setfill(
'0') << setw(6) << myid;
452 ee_name <<
"elastic_energy." << setfill(
'0') << setw(6) << myid;
454 ofstream mesh_ofs(mesh_name.str().c_str());
455 mesh_ofs.precision(8);
456 pmesh->
Print(mesh_ofs);
458 ofstream velo_ofs(velo_name.str().c_str());
459 velo_ofs.precision(8);
461 ofstream ee_ofs(ee_name.str().c_str());
463 oper->GetElasticEnergyDensity(x_gf, w_gf);
494 out <<
"solution\n" << *mesh << *field;
500 out <<
"window_size 800 800\n";
501 out <<
"window_title '" << field_name <<
"'\n";
508 out <<
"autoscale value\n";
515 ReducedSystemOperator::ReducedSystemOperator(
518 :
Operator(M_->ParFESpace()->TrueVSize()), M(M_), S(S_), H(H_),
519 Jacobian(NULL), dt(0.0), v(NULL), x(NULL), w(height), z(height),
520 ess_tdof_list(ess_tdof_list_)
523 void ReducedSystemOperator::SetParameters(
double dt_,
const Vector *v_,
526 dt = dt_; v = v_; x = x_;
535 M->TrueAddMult(k, y);
536 S->TrueAddMult(w, y);
540 Operator &ReducedSystemOperator::GetGradient(
const Vector &k)
const
546 localJ->
Add(dt*dt, H->GetLocalGradient(z));
549 Jacobian = M->ParallelAssemble(localJ);
556 ReducedSystemOperator::~ReducedSystemOperator()
564 double mu,
double K,
bool use_petsc,
565 bool use_petsc_factory)
567 M(&fespace), S(&fespace), H(&fespace),
568 viscosity(visc), M_solver(f.GetComm()),
569 newton_solver(f.GetComm()), pnewton_solver(NULL), z(height/2)
571 const double rel_tol = 1e-8;
572 const int skip_zero_entries = 0;
574 const double ref_density = 1.0;
577 M.Assemble(skip_zero_entries);
578 M.Finalize(skip_zero_entries);
579 Mmat = M.ParallelAssemble();
580 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
584 M_solver.iterative_mode =
false;
585 M_solver.SetRelTol(rel_tol);
586 M_solver.SetAbsTol(0.0);
587 M_solver.SetMaxIter(30);
588 M_solver.SetPrintLevel(0);
589 M_prec.SetType(HypreSmoother::Jacobi);
590 M_solver.SetPreconditioner(M_prec);
591 M_solver.SetOperator(*Mmat);
595 H.SetEssentialTrueDofs(ess_tdof_list);
599 S.Assemble(skip_zero_entries);
600 S.Finalize(skip_zero_entries);
602 reduced_oper =
new ReducedSystemOperator(&M, &S, &H, ess_tdof_list);
606 J_hypreSmoother->
SetType(HypreSmoother::l1Jacobi);
608 J_prec = J_hypreSmoother;
621 newton_solver.SetSolver(*J_solver);
622 newton_solver.SetOperator(*reduced_oper);
623 newton_solver.SetPrintLevel(1);
624 newton_solver.SetRelTol(rel_tol);
625 newton_solver.SetAbsTol(0.0);
626 newton_solver.SetMaxIter(10);
639 if (use_petsc_factory)
641 J_factory =
new PreconditionerFactory(*reduced_oper,
"JFNK preconditioner");
642 pnewton_solver->SetPreconditionerFactory(J_factory);
644 pnewton_solver->SetPrintLevel(1);
645 pnewton_solver->SetRelTol(rel_tol);
646 pnewton_solver->SetAbsTol(0.0);
647 pnewton_solver->SetMaxIter(10);
661 if (viscosity != 0.0)
667 M_solver.Mult(z, dv_dt);
672 void HyperelasticOperator::ImplicitSolve(
const double dt,
687 reduced_oper->SetParameters(dt, &v, &x);
691 newton_solver.Mult(zero, dv_dt);
692 MFEM_VERIFY(newton_solver.GetConverged(),
693 "Newton solver did not converge.");
697 pnewton_solver->Mult(zero, dv_dt);
698 MFEM_VERIFY(pnewton_solver->GetConverged(),
699 "Newton solver did not converge.");
701 add(v, dt, dv_dt, dx_dt);
704 double HyperelasticOperator::ElasticEnergy(
const ParGridFunction &x)
const
706 return H.GetEnergy(x);
709 double HyperelasticOperator::KineticEnergy(
const ParGridFunction &v)
const
711 double loc_energy = 0.5*M.InnerProduct(v, v);
713 MPI_Allreduce(&loc_energy, &energy, 1, MPI_DOUBLE, MPI_SUM,
718 void HyperelasticOperator::GetElasticEnergyDensity(
721 ElasticEnergyCoefficient w_coeff(*model, x);
725 HyperelasticOperator::~HyperelasticOperator()
733 delete pnewton_solver;
753 model.SetTransformation(T);
756 return model.EvalW(J)/J.Det();
770 const double s = 0.1/64.;
773 v(dim-1) = s*x(0)*x(0)*(8.0-x(0));
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
double Eval(ElementTransformation &T, const IntegrationPoint &ip, double t)
Evaluate the coefficient in the element described by T at the point ip at time t. ...
void Add(const int i, const int j, const double a)
void InitialDeformation(const Vector &x, Vector &y)
Conjugate gradient method.
Class for grid function - Vector with associated FE space.
Abstract class for PETSc's preconditioners.
A class to handle Vectors in a block fashion.
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Subclass constant coefficient.
Wrapper for PETSc's matrix class.
Base abstract class for time dependent operators.
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Pointer to an Operator of a specified type.
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Data type dense matrix using column-major storage.
int Size() const
Returns the size of the vector.
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
virtual void Save(std::ostream &out) const
Abstract parallel finite element space.
virtual void ProjectCoefficient(Coefficient &coeff)
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
int main(int argc, char *argv[])
Backward Euler ODE solver. L-stable.
double * GetData() const
Return a pointer to the beginning of the Vector data.
void add(const Vector &v1, const Vector &v2, Vector &v)
void InitialVelocity(const Vector &x, Vector &v)
void SetPositiveDiagonal(bool pos=true)
After computing l1-norms, replace them with their absolute values.
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
void SetPrintLevel(int print_lvl)
Mesh * GetMesh() const
Returns the mesh.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
void MakeTRef(FiniteElementSpace *f, double *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
void SetMaxIter(int max_it)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
HYPRE_Int GlobalTrueVSize() const
Newton's method for solving F(x)=b for a given operator F.
virtual void Print(std::ostream &out=mfem::out) const
Parallel smoothers in hypre.
void PrintUsage(std::ostream &out) const
Abstract class for PETSc's nonlinear solvers.
int SpaceDimension() const
The classical explicit forth-order Runge-Kutta method, RK4.
void SetAbsTol(double atol)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void SetRelTol(double rtol)
Base class Coefficient that may optionally depend on time.
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)
Third-order, strong stability preserving (SSP) Runge-Kutta method.
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad) const
Implicit midpoint method. A-stable, not L-stable.
Class for integration point with weight.
void MFEMInitializePetsc()
Convenience functions to initialize/finalize PETSc.
void PrintOptions(std::ostream &out) const
Abstract class for hyperelastic models.
int open(const char hostname[], int port)
void Get(OpType *&A) const
Return the Operator pointer statically cast to a given OpType.
void GetNodes(Vector &node_coord) const
Arbitrary order H1-conforming (continuous) finite elements.
Class for parallel grid function.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
The classical forward Euler method.
Wrapper for hypre's ParCSR matrix class.
Class for parallel meshes.
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Arbitrary order "L2-conforming" discontinuous finite elements.