45 class BackwardEulerOperator;
73 BackwardEulerOperator *backward_euler_oper;
90 virtual void ImplicitSolve(
const double dt,
const Vector &x,
Vector &k);
96 virtual ~HyperelasticOperator();
103 class BackwardEulerOperator :
public Operator
116 void SetParameters(
double dt_,
const Vector *v_,
const Vector *x_);
119 virtual ~BackwardEulerOperator();
124 class ElasticEnergyCoefficient :
public Coefficient
133 : model(m), x(x_) { }
135 virtual ~ElasticEnergyCoefficient() { }
144 bool init_vis =
false);
147 int main(
int argc,
char *argv[])
151 MPI_Init(&argc, &argv);
152 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
153 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
156 const char *mesh_file =
"../data/beam-quad.mesh";
157 int ser_ref_levels = 2;
158 int par_ref_levels = 0;
160 int ode_solver_type = 3;
161 double t_final = 300.0;
164 bool visualization =
true;
168 args.
AddOption(&mesh_file,
"-m",
"--mesh",
169 "Mesh file to use.");
170 args.
AddOption(&ser_ref_levels,
"-rs",
"--refine-serial",
171 "Number of times to refine the mesh uniformly in serial.");
172 args.
AddOption(&par_ref_levels,
"-rp",
"--refine-parallel",
173 "Number of times to refine the mesh uniformly in parallel.");
175 "Order (degree) of the finite elements.");
176 args.
AddOption(&ode_solver_type,
"-s",
"--ode-solver",
177 "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
178 "\t 11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4.");
179 args.
AddOption(&t_final,
"-tf",
"--t-final",
180 "Final time; start time is 0.");
181 args.
AddOption(&dt,
"-dt",
"--time-step",
183 args.
AddOption(&visc,
"-v",
"--viscosity",
184 "Viscosity coefficient.");
185 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
186 "--no-visualization",
187 "Enable or disable GLVis visualization.");
188 args.
AddOption(&vis_steps,
"-vs",
"--visualization-steps",
189 "Visualize every n-th timestep.");
205 ifstream imesh(mesh_file);
209 cerr <<
"\nCan not open mesh file: " << mesh_file <<
'\n' << endl;
213 mesh =
new Mesh(imesh, 1, 1);
221 switch (ode_solver_type)
229 case 12: ode_solver =
new RK2Solver(0.5);
break;
231 case 14: ode_solver =
new RK4Solver;
break;
238 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
246 for (
int lev = 0; lev < ser_ref_levels; lev++)
254 for (
int lev = 0; lev < par_ref_levels; lev++)
268 cout <<
"Number of velocity/deformation unknowns: " << glob_size << endl;
272 true_offset[1] = true_size;
273 true_offset[2] = 2*true_size;
288 v_gf.ProjectCoefficient(velo);
301 HyperelasticOperator oper(fespace, ess_bdr, visc);
306 char vishost[] =
"localhost";
308 vis_v.
open(vishost, visport);
310 visualize(vis_v, pmesh, &x_gf, &v_gf,
"Velocity",
true);
314 vis_w.
open(vishost, visport);
317 oper.GetElasticEnergyDensity(x_gf, w_gf);
319 visualize(vis_w, pmesh, &x_gf, &w_gf,
"Elastic energy density",
true);
323 double ee0 = oper.ElasticEnergy(x_gf);
324 double ke0 = oper.KineticEnergy(v_gf);
327 cout <<
"initial elastic energy (EE) = " << ee0 << endl;
328 cout <<
"initial kinetic energy (KE) = " << ke0 << endl;
329 cout <<
"initial total energy (TE) = " << (ee0 + ke0) << endl;
334 ode_solver->
Init(oper);
337 bool last_step =
false;
338 for (
int ti = 1; !last_step; ti++)
340 if (t + dt >= t_final - dt/2)
343 ode_solver->
Step(vx, t, dt);
345 if (last_step || (ti % vis_steps) == 0)
350 double ee = oper.ElasticEnergy(x_gf);
351 double ke = oper.KineticEnergy(v_gf);
354 cout <<
"step " << ti <<
", t = " << t <<
", EE = " << ee
355 <<
", KE = " << ke <<
", ΔTE = " << (ee+ke)-(ee0+ke0) << endl;
362 oper.GetElasticEnergyDensity(x_gf, w_gf);
375 ostringstream mesh_name, velo_name, ee_name;
376 mesh_name <<
"deformed." << setfill(
'0') << setw(6) << myid;
377 velo_name <<
"velocity." << setfill(
'0') << setw(6) << myid;
378 ee_name <<
"elastic_energy." << setfill(
'0') << setw(6) << myid;
380 ofstream mesh_ofs(mesh_name.str().c_str());
381 mesh_ofs.precision(8);
382 pmesh->
Print(mesh_ofs);
384 ofstream velo_ofs(velo_name.str().c_str());
385 velo_ofs.precision(8);
387 ofstream ee_ofs(ee_name.str().c_str());
389 oper.GetElasticEnergyDensity(x_gf, w_gf);
414 out <<
"solution\n" << *mesh << *field;
420 out <<
"window_size 800 800\n";
421 out <<
"window_title '" << field_name <<
"'\n";
428 out <<
"autoscale value\n";
434 BackwardEulerOperator::BackwardEulerOperator(
436 :
Operator(M_->ParFESpace()->TrueVSize()), M(M_), S(S_), H(H_),
437 Jacobian(NULL), v(NULL), x(NULL), dt(0.0), w(height), z(height)
440 void BackwardEulerOperator::SetParameters(
double dt_,
const Vector *v_,
443 dt = dt_; v = v_; x = x_;
452 M->TrueAddMult(k, y);
453 S->TrueAddMult(w, y);
456 Operator &BackwardEulerOperator::GetGradient(
const Vector &k)
const
462 localJ->
Add(dt*dt, H->GetLocalGradient(z));
463 Jacobian = M->ParallelAssemble(localJ);
468 BackwardEulerOperator::~BackwardEulerOperator()
477 M(&fespace), S(&fespace), H(&fespace), M_solver(f.GetComm()),
478 newton_solver(f.GetComm()), z(height/2)
480 const double rel_tol = 1e-8;
481 const int skip_zero_entries = 0;
483 const double ref_density = 1.0;
486 M.Assemble(skip_zero_entries);
487 M.EliminateEssentialBC(ess_bdr);
488 M.Finalize(skip_zero_entries);
489 Mmat = M.ParallelAssemble();
491 M_solver.iterative_mode =
false;
492 M_solver.SetRelTol(rel_tol);
493 M_solver.SetAbsTol(0.0);
494 M_solver.SetMaxIter(30);
495 M_solver.SetPrintLevel(0);
496 M_prec.SetType(HypreSmoother::Jacobi);
497 M_solver.SetPreconditioner(M_prec);
498 M_solver.SetOperator(*Mmat);
504 H.SetEssentialBC(ess_bdr);
509 S.Assemble(skip_zero_entries);
510 S.EliminateEssentialBC(ess_bdr);
511 S.Finalize(skip_zero_entries);
513 backward_euler_oper =
new BackwardEulerOperator(&M, &S, &H);
516 J_hypreSmoother->
SetType(HypreSmoother::l1Jacobi);
517 J_prec = J_hypreSmoother;
528 newton_solver.SetSolver(*J_solver);
529 newton_solver.SetOperator(*backward_euler_oper);
530 newton_solver.SetPrintLevel(1);
531 newton_solver.SetRelTol(rel_tol);
532 newton_solver.SetAbsTol(0.0);
533 newton_solver.SetMaxIter(10);
546 if (viscosity != 0.0)
549 M_solver.Mult(z, dv_dt);
554 void HyperelasticOperator::ImplicitSolve(
const double dt,
569 backward_euler_oper->SetParameters(dt, &v, &x);
571 newton_solver.Mult(zero, dv_dt);
572 add(v, dt, dv_dt, dx_dt);
574 MFEM_VERIFY(newton_solver.GetConverged(),
"Newton Solver did not converge.");
579 return H.GetEnergy(x);
584 double loc_energy = 0.5*M.InnerProduct(v, v);
586 MPI_Allreduce(&loc_energy, &energy, 1, MPI_DOUBLE, MPI_SUM,
591 void HyperelasticOperator::GetElasticEnergyDensity(
594 ElasticEnergyCoefficient w_coeff(*model, x);
598 HyperelasticOperator::~HyperelasticOperator()
601 delete backward_euler_oper;
611 model.SetTransformation(T);
614 return model.EvalW(J)/J.Det();
627 const int dim = x.
Size();
628 const double s = 0.1/64.;
631 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)
double Eval(ElementTransformation &T, const IntegrationPoint &ip, double 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.
Subclass constant coefficient.
Base abstract class for time dependent operators: (x,t) -> f(x,t)
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
virtual void Step(Vector &x, double &t, double &dt)=0
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.
void ProjectCoefficient(Coefficient &coeff)
bool iterative_mode
If true, use the second argument of Mult as an initial guess.
virtual void Init(TimeDependentOperator &_f)
Hyperelastic integrator for any given HyperelasticModel.
Backward Euler ODE solver. L-stable.
void add(const Vector &v1, const Vector &v2, Vector &v)
void InitialVelocity(const Vector &x, Vector &v)
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
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 SetMaxIter(int max_it)
Parallel smoothers in hypre.
void PrintUsage(std::ostream &out) const
void GetTrueDofs(Vector &tv) const
Returns the true dofs in a HypreParVector.
int SpaceDimension() const
The classical explicit forth-order Runge-Kutta method, RK4.
void SetAbsTol(double atol)
Array< int > bdr_attributes
int main(int argc, char *argv[])
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 Distribute(const Vector *tv)
Implicit midpoint method. A-stable, not L-stable.
Class for integration point with weight.
void PrintOptions(std::ostream &out) const
Abstract class for hyperelastic models.
int open(const char hostname[], int port)
void GetNodes(Vector &node_coord) const
Arbitrary order H1-conforming (continuous) finite elements.
Class for parallel grid function.
The classical forward Euler method.
Wrapper for hypre's ParCSR matrix class.
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad)
Vector & GetBlock(int i)
Get the i-th vector in the block.
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.
void Neg()
(*this) = -(*this)
virtual void Print(std::ostream &out=std::cout) const