56 class JacobianPreconditioner :
public Solver
83 virtual void SetOperator(
const Operator &op);
85 virtual ~JacobianPreconditioner();
92 class RubberOperator :
public Operator
121 Array<int> &block_trueOffsets,
double rel_tol,
double abs_tol,
129 void Solve(
Vector &xp)
const;
131 virtual ~RubberOperator();
137 bool init_vis =
false);
144 int main(
int argc,
char *argv[])
148 MPI_Init(&argc, &argv);
149 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
150 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
153 const char *mesh_file =
"../data/beam-tet.mesh";
154 int ser_ref_levels = 0;
155 int par_ref_levels = 0;
157 bool visualization =
true;
158 double newton_rel_tol = 1e-4;
159 double newton_abs_tol = 1e-6;
160 int newton_iter = 500;
164 args.
AddOption(&mesh_file,
"-m",
"--mesh",
165 "Mesh file to use.");
166 args.
AddOption(&ser_ref_levels,
"-rs",
"--refine-serial",
167 "Number of times to refine the mesh uniformly in serial.");
168 args.
AddOption(&par_ref_levels,
"-rp",
"--refine-parallel",
169 "Number of times to refine the mesh uniformly in parallel.");
171 "Order (degree) of the finite elements.");
172 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
173 "--no-visualization",
174 "Enable or disable GLVis visualization.");
175 args.
AddOption(&newton_rel_tol,
"-rel",
"--relative-tolerance",
176 "Relative tolerance for the Newton solve.");
177 args.
AddOption(&newton_abs_tol,
"-abs",
"--absolute-tolerance",
178 "Absolute tolerance for the Newton solve.");
179 args.
AddOption(&newton_iter,
"-it",
"--newton-iterations",
180 "Maximum iterations for the Newton solve.");
181 args.
AddOption(&mu,
"-mu",
"--shear-modulus",
182 "Shear modulus for the neo-Hookean material.");
201 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
207 for (
int lev = 0; lev < ser_ref_levels; lev++)
217 for (
int lev = 0; lev < par_ref_levels; lev++)
235 spaces[0] = &R_space;
236 spaces[1] = &W_space;
252 ess_bdr[0] = &ess_bdr_u;
253 ess_bdr[1] = &ess_bdr_p;
258 std::cout <<
"***********************************************************\n";
259 std::cout <<
"dim(u) = " << glob_R_size <<
"\n";
260 std::cout <<
"dim(p) = " << glob_W_size <<
"\n";
261 std::cout <<
"dim(u+p) = " << glob_R_size + glob_W_size <<
"\n";
262 std::cout <<
"***********************************************************\n";
267 block_trueOffsets[0] = 0;
268 block_trueOffsets[1] = R_space.
TrueVSize();
269 block_trueOffsets[2] = W_space.
TrueVSize();
293 RubberOperator oper(spaces, ess_bdr, block_trueOffsets,
294 newton_rel_tol, newton_abs_tol, newton_iter, c_mu);
310 char vishost[] =
"localhost";
312 vis_u.
open(vishost, visport);
314 visualize(vis_u, pmesh, &x_gf, &x_def,
"Deformation",
true);
318 vis_p.
open(vishost, visport);
320 visualize(vis_p, pmesh, &x_gf, &p_gf,
"Pressure",
true);
329 ostringstream mesh_name, pressure_name, deformation_name;
330 mesh_name <<
"mesh." << setfill(
'0') << setw(6) << myid;
331 pressure_name <<
"pressure." << setfill(
'0') << setw(6) << myid;
332 deformation_name <<
"deformation." << setfill(
'0') << setw(6) << myid;
334 ofstream mesh_ofs(mesh_name.str().c_str());
335 mesh_ofs.precision(8);
336 pmesh->
Print(mesh_ofs);
338 ofstream pressure_ofs(pressure_name.str().c_str());
339 pressure_ofs.precision(8);
340 p_gf.
Save(pressure_ofs);
342 ofstream deformation_ofs(deformation_name.str().c_str());
343 deformation_ofs.precision(8);
344 x_def.
Save(deformation_ofs);
360 :
Solver(offsets[2]), block_trueOffsets(offsets), pressure_mass(&mass)
371 mass_prec = mass_prec_amg;
382 mass_pcg = mass_pcg_iter;
394 block_trueOffsets[1]-block_trueOffsets[0]);
396 block_trueOffsets[2]-block_trueOffsets[1]);
399 block_trueOffsets[1]-block_trueOffsets[0]);
401 block_trueOffsets[2]-block_trueOffsets[1]);
403 Vector temp(block_trueOffsets[1]-block_trueOffsets[0]);
404 Vector temp2(block_trueOffsets[1]-block_trueOffsets[0]);
407 mass_pcg->Mult(pres_in, pres_out);
410 jacobian->GetBlock(0,1).Mult(pres_out, temp);
413 stiff_pcg->Mult(temp2, disp_out);
416 void JacobianPreconditioner::SetOperator(
const Operator &op)
421 if (stiff_prec == NULL)
427 stiff_prec = stiff_prec_amg;
437 stiff_pcg = stiff_pcg_iter;
445 JacobianPreconditioner::~JacobianPreconditioner()
461 :
Operator(fes[0]->TrueVSize() + fes[1]->TrueVSize()),
462 newton_solver(fes[0]->GetComm()), mu(c_mu), block_trueOffsets(trueOffsets)
476 Hform->SetEssentialBC(ess_bdr, rhs);
488 mass.SetOperatorOwner(
false);
489 pressure_mass = mass.Ptr();
492 JacobianPreconditioner *jac_prec =
493 new JacobianPreconditioner(fes, *pressure_mass, block_trueOffsets);
508 newton_solver.SetSolver(*j_solver);
509 newton_solver.SetOperator(*
this);
510 newton_solver.SetPrintLevel(1);
511 newton_solver.SetRelTol(rel_tol);
512 newton_solver.SetAbsTol(abs_tol);
513 newton_solver.SetMaxIter(iter);
517 void RubberOperator::Solve(
Vector &xp)
const
520 newton_solver.Mult(zero, xp);
521 MFEM_VERIFY(newton_solver.GetConverged(),
522 "Newton Solver did not converge.");
537 RubberOperator::~RubberOperator()
540 delete pressure_mass;
561 out <<
"solution\n" << *mesh << *field;
567 out <<
"window_size 800 800\n";
568 out <<
"window_title '" << field_name <<
"'\n";
575 out <<
"autoscale value\n";
591 y[1] = x[1] + 0.25*x[0];
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
void InitialDeformation(const Vector &x, Vector &y)
Conjugate gradient method.
void ReferenceConfiguration(const Vector &x, Vector &y)
Class for grid function - Vector with associated FE space.
A class to handle Vectors in a block fashion.
Subclass constant coefficient.
virtual Operator & GetGradient(const Vector &x) const
Evaluate the gradient operator at the point x. The default behavior in class Operator is to generate ...
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.
void Copy(Array ©) const
Create a copy of the current array.
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[])
double * GetData() const
Return a pointer to the beginning of the Vector data.
The BoomerAMG solver in hypre.
void SetPrintLevel(int print_lvl)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Mesh * GetMesh() const
Returns the mesh.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
void SetPrintLevel(int print_level)
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
void SetElasticityOptions(ParFiniteElementSpace *fespace)
void PrintUsage(std::ostream &out) const
int SpaceDimension() const
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 subtract(const Vector &x, const Vector &y, Vector &z)
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)
void PartialSum()
Partial Sum.
void Distribute(const Vector *tv)
void PrintOptions(std::ostream &out) const
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
int open(const char hostname[], int port)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
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...
A class to handle Block systems in a matrix-free implementation.
Class for parallel meshes.
Vector & GetBlock(int i)
Get the i-th vector in the block.