45 class BackwardEulerOperator;
72 BackwardEulerOperator *backward_euler_oper;
89 virtual void ImplicitSolve(
const double dt,
const Vector &x,
Vector &k);
91 double ElasticEnergy(
Vector &x)
const;
92 double KineticEnergy(
Vector &v)
const;
95 virtual ~HyperelasticOperator();
102 class BackwardEulerOperator :
public Operator
114 void SetParameters(
double dt_,
const Vector *v_,
const Vector *x_);
117 virtual ~BackwardEulerOperator();
122 class ElasticEnergyCoefficient :
public Coefficient
131 : model(m), x(x_) { }
133 virtual ~ElasticEnergyCoefficient() { }
142 bool init_vis =
false);
145 int main(
int argc,
char *argv[])
148 const char *mesh_file =
"../data/beam-quad.mesh";
151 int ode_solver_type = 3;
152 double t_final = 300.0;
155 bool visualization =
true;
159 args.
AddOption(&mesh_file,
"-m",
"--mesh",
160 "Mesh file to use.");
161 args.
AddOption(&ref_levels,
"-r",
"--refine",
162 "Number of times to refine the mesh uniformly.");
164 "Order (degree) of the finite elements.");
165 args.
AddOption(&ode_solver_type,
"-s",
"--ode-solver",
166 "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
167 "\t 11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4.");
168 args.
AddOption(&t_final,
"-tf",
"--t-final",
169 "Final time; start time is 0.");
170 args.
AddOption(&dt,
"-dt",
"--time-step",
172 args.
AddOption(&visc,
"-v",
"--viscosity",
173 "Viscosity coefficient.");
174 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
175 "--no-visualization",
176 "Enable or disable GLVis visualization.");
177 args.
AddOption(&vis_steps,
"-vs",
"--visualization-steps",
178 "Visualize every n-th timestep.");
190 ifstream imesh(mesh_file);
193 cerr <<
"Can not open mesh: " << mesh_file << endl;
196 mesh =
new Mesh(imesh, 1, 1);
204 switch (ode_solver_type)
212 case 12: ode_solver =
new RK2Solver(0.5);
break;
214 case 14: ode_solver =
new RK4Solver;
break;
220 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
227 for (
int lev = 0; lev < ref_levels; lev++)
240 cout <<
"Number of velocity/deformation unknowns: " << fe_size << endl;
243 fe_offset[1] = fe_size;
244 fe_offset[2] = 2*fe_size;
271 HyperelasticOperator oper(fespace, ess_bdr, visc);
276 char vishost[] =
"localhost";
278 vis_v.
open(vishost, visport);
280 visualize(vis_v, mesh, &x, &v,
"Velocity",
true);
281 vis_w.
open(vishost, visport);
284 oper.GetElasticEnergyDensity(x, w);
286 visualize(vis_w, mesh, &x, &w,
"Elastic energy density",
true);
290 double ee0 = oper.ElasticEnergy(x);
291 double ke0 = oper.KineticEnergy(v);
292 cout <<
"initial elastic energy (EE) = " << ee0 << endl;
293 cout <<
"initial kinetic energy (KE) = " << ke0 << endl;
294 cout <<
"initial total energy (TE) = " << (ee0 + ke0) << endl;
298 ode_solver->
Init(oper);
301 bool last_step =
false;
302 for (
int ti = 1; !last_step; ti++)
304 if (t + dt >= t_final - dt/2)
307 ode_solver->
Step(vx, t, dt);
309 if (last_step || (ti % vis_steps) == 0)
311 double ee = oper.ElasticEnergy(x);
312 double ke = oper.KineticEnergy(v);
314 cout <<
"step " << ti <<
", t = " << t <<
", EE = " << ee <<
", KE = "
315 << ke <<
", ΔTE = " << (ee+ke)-(ee0+ke0) << endl;
322 oper.GetElasticEnergyDensity(x, w);
334 ofstream mesh_ofs(
"deformed.mesh");
335 mesh_ofs.precision(8);
336 mesh->
Print(mesh_ofs);
338 ofstream velo_ofs(
"velocity.sol");
339 velo_ofs.precision(8);
341 ofstream ee_ofs(
"elastic_energy.sol");
343 oper.GetElasticEnergyDensity(x, w);
355 GridFunction *field,
const char *field_name,
bool init_vis)
365 out <<
"solution\n" << *mesh << *field;
371 out <<
"window_size 800 800\n";
372 out <<
"window_title '" << field_name <<
"'\n";
379 out <<
"autoscale value\n";
385 BackwardEulerOperator::BackwardEulerOperator(
387 :
Operator(M_->Height()), M(M_), S(S_), H(H_), Jacobian(NULL),
388 v(NULL), x(NULL), dt(0.0), w(height), z(height)
391 void BackwardEulerOperator::SetParameters(
double dt_,
const Vector *v_,
394 dt = dt_; v = v_; x = x_;
407 Operator &BackwardEulerOperator::GetGradient(
const Vector &k)
const
410 Jacobian =
Add(1.0, M->SpMat(), dt, S->SpMat());
414 Jacobian->
Add(dt*dt, *grad_H);
418 BackwardEulerOperator::~BackwardEulerOperator()
427 M(&fespace), S(&fespace), H(&fespace), z(height/2)
429 const double rel_tol = 1e-8;
430 const int skip_zero_entries = 0;
432 const double ref_density = 1.0;
435 M.Assemble(skip_zero_entries);
436 M.EliminateEssentialBC(ess_bdr);
437 M.Finalize(skip_zero_entries);
439 M_solver.iterative_mode =
false;
440 M_solver.SetRelTol(rel_tol);
441 M_solver.SetAbsTol(0.0);
442 M_solver.SetMaxIter(30);
443 M_solver.SetPrintLevel(0);
444 M_solver.SetPreconditioner(M_prec);
445 M_solver.SetOperator(M.SpMat());
451 H.SetEssentialBC(ess_bdr);
456 S.Assemble(skip_zero_entries);
457 S.EliminateEssentialBC(ess_bdr);
458 S.Finalize(skip_zero_entries);
460 backward_euler_oper =
new BackwardEulerOperator(&M, &S, &H);
462 #ifndef MFEM_USE_SUITESPARSE
477 newton_solver.SetSolver(*J_solver);
478 newton_solver.SetOperator(*backward_euler_oper);
479 newton_solver.SetPrintLevel(1);
480 newton_solver.SetRelTol(rel_tol);
481 newton_solver.SetAbsTol(0.0);
482 newton_solver.SetMaxIter(10);
495 if (viscosity != 0.0)
498 M_solver.Mult(z, dv_dt);
503 void HyperelasticOperator::ImplicitSolve(
const double dt,
518 backward_euler_oper->SetParameters(dt, &v, &x);
520 newton_solver.Mult(zero, dv_dt);
521 add(v, dt, dv_dt, dx_dt);
523 MFEM_VERIFY(newton_solver.GetConverged(),
"Newton Solver did not converge.");
526 double HyperelasticOperator::ElasticEnergy(
Vector &x)
const
528 return H.GetEnergy(x);
531 double HyperelasticOperator::KineticEnergy(
Vector &v)
const
533 return 0.5*M.InnerProduct(v, v);
536 void HyperelasticOperator::GetElasticEnergyDensity(
539 ElasticEnergyCoefficient w_coeff(*model, x);
543 HyperelasticOperator::~HyperelasticOperator()
546 delete backward_euler_oper;
555 model.SetTransformation(T);
558 return model.EvalW(J)/J.Det();
571 const int dim = x.
Size();
572 const double s = 0.1/64.;
575 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.
Data type for scaled Jacobi-type smoother of sparse matrix.
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)
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)
Direct sparse solver using UMFPACK.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
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)
void PrintUsage(std::ostream &out) const
int SpaceDimension() const
virtual void Print(std::ostream &out=std::cout) const
Print the mesh to the given stream using the default MFEM mesh format.
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)
Abstract finite element space.
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.
Implicit midpoint method. A-stable, not L-stable.
Class for integration point with weight.
void PrintOptions(std::ostream &out) const
void ProjectCoefficient(Coefficient &coeff)
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.
The classical forward Euler method.
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad)
Vector & GetBlock(int i)
Get the i-th vector in the block.
Arbitrary order "L2-conforming" discontinuous finite elements.
void Neg()
(*this) = -(*this)