33 #include "constitutive.h"
36 #include "meshes/rectangular_quadmesh.h"
39 #include "multi_physics/pseudo_elastic_preconditioner.h"
44 using namespace oomph;
53 template <
class ELEMENT>
55 public virtual ELEMENT
68 return 2*ELEMENT::dim();
85 std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list)
const
88 std::pair<unsigned,unsigned> dof_lookup;
91 const unsigned n_node = this->nnode();
94 const unsigned n_position_type = ELEMENT::nnodal_position_type();
95 const unsigned nodal_dim = ELEMENT::nodal_dimension();
101 for(
unsigned n=0;n<n_node;n++)
104 if (this->node_pt(n)->nvalue() != this->required_nvalue(n))
106 offset = ELEMENT::dim();
110 for(
unsigned k=0;k<n_position_type;k++)
113 for(
unsigned i=0;i<nodal_dim;i++)
116 local_unknown = ELEMENT::position_local_eqn(n,k,i);
117 if (local_unknown >= 0)
121 dof_lookup.first = this->eqn_number(local_unknown);
122 dof_lookup.second = offset+i;
125 dof_lookup_list.push_front(dof_lookup);
140 template<
class ELEMENT>
142 public virtual FaceGeometry<ELEMENT>
176 MatrixBasedDiagPreconditioner* prec_pt=
177 new MatrixBasedDiagPreconditioner;
200 bool operator()(FiniteElement*
const& el1_pt, FiniteElement*
const& el2_pt)
203 return el1_pt->node_pt(0)->x(0) < el2_pt->node_pt(0)->x(0);
227 BrokenCopy::broken_copy(
"WarpedLine");
233 BrokenCopy::broken_assign(
"WarpedLine");
241 void position(
const Vector<double>& zeta, Vector<double>& r)
const
244 r[0] = zeta[0]+5.0*Ampl*zeta[0]*(zeta[0]-1.0)*(zeta[0]-0.7);
245 r[1] = 1.0+Ampl*0.5*(1.0-cos(2.0*MathematicalConstants::Pi*zeta[0]));
251 void position(
const unsigned& t,
const Vector<double>& zeta,
252 Vector<double>& r)
const
304 template<
class ELEMENT>
321 {
return Solid_mesh_pt;}
343 ElasticRefineableRectangularQuadMesh<ELEMENT>* Solid_mesh_pt;
346 SolidMesh* Lagrange_multiplier_mesh_pt;
357 template<
class ELEMENT>
377 solid_mesh_pt() =
new ElasticRefineableRectangularQuadMesh<ELEMENT>(
381 solid_mesh_pt()->spatial_error_estimator_pt()=
new Z2ErrorEstimator;
385 unsigned n_element =solid_mesh_pt()->nelement();
386 for(
unsigned i=0;i<n_element;i++)
389 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(solid_mesh_pt()->element_pt(i));
392 el_pt->constitutive_law_pt()=&Global_Physical_Variables::Constitutive_law;
396 solid_mesh_pt()->refine_uniformly();
400 Lagrange_multiplier_mesh_pt=
new SolidMesh;
401 create_lagrange_multiplier_elements();
404 add_sub_mesh(solid_mesh_pt());
407 add_sub_mesh(Lagrange_multiplier_mesh_pt);
413 for (
unsigned b=0;b<4;b++)
417 unsigned n_side = solid_mesh_pt()->nboundary_node(b);
420 for(
unsigned i=0;i<n_side;i++)
422 solid_mesh_pt()->boundary_node_pt(b,i)->pin_position(0);
423 solid_mesh_pt()->boundary_node_pt(b,i)->pin_position(1);
429 PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
430 solid_mesh_pt()->element_pt());
433 cout <<
"Number of dofs: " << assign_eqn_numbers() << std::endl;
436 Doc_info.set_directory(
"RESLT");
439 IterativeLinearSolver* solver_pt=0;
442 #ifdef OOMPH_HAS_TRILINOS
445 solver_pt =
new TrilinosAztecOOSolver;
448 dynamic_cast<TrilinosAztecOOSolver*
>(solver_pt)->solver_type()
449 = TrilinosAztecOOSolver::GMRES;
454 solver_pt =
new GMRES<CRDoubleMatrix>;
459 linear_solver_pt() = solver_pt;
462 PseudoElasticPreconditioner * prec_pt =
new PseudoElasticPreconditioner;
465 prec_pt->set_elastic_mesh(Solid_mesh_pt);
466 prec_pt->set_lagrange_multiplier_mesh(Lagrange_multiplier_mesh_pt);
469 solver_pt->preconditioner_pt() = prec_pt;
472 if (CommandLineArgs::command_line_flag_has_been_set
473 (
"--block_upper_for_elastic_block"))
475 prec_pt->elastic_preconditioner_type()=
476 PseudoElasticPreconditioner::Block_upper_triangular_preconditioner;
479 #ifdef OOMPH_HAS_HYPRE
483 if (CommandLineArgs::command_line_flag_has_been_set
484 (
"--hypre_for_elastic_blocks"))
486 prec_pt->set_elastic_subsidiary_preconditioner
487 (Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper::
488 get_elastic_preconditioner_hypre);
494 #ifdef OOMPH_HAS_TRILINOS
498 if (CommandLineArgs::command_line_flag_has_been_set
499 (
"--trilinos_cg_for_lagrange_multiplier_blocks"))
501 prec_pt->set_lagrange_multiplier_subsidiary_preconditioner
502 (Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper::
503 get_lagrange_multiplier_preconditioner);
510 if (CommandLineArgs::command_line_flag_has_been_set
511 (
"--diagonal_scaling_for_elastic_blocks"))
513 prec_pt->set_elastic_subsidiary_preconditioner(
524 template<
class ELEMENT>
528 delete_lagrange_multiplier_elements();
531 rebuild_global_mesh();
541 template<
class ELEMENT>
547 create_lagrange_multiplier_elements();
550 rebuild_global_mesh();
553 PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
554 solid_mesh_pt()->element_pt());
563 template<
class ELEMENT>
571 unsigned n_element = solid_mesh_pt()->nboundary_element(b);
574 for(
unsigned e=0;e<n_element;e++)
577 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
578 solid_mesh_pt()->boundary_element_pt(b,e));
581 int face_index = solid_mesh_pt()->face_index_at_boundary(b,e);
584 Lagrange_multiplier_mesh_pt->add_element_pt(
585 new ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>(
586 bulk_elem_pt,face_index));
592 n_element=Lagrange_multiplier_mesh_pt->nelement();
593 for(
unsigned i=0;i<n_element;i++)
596 ImposeDisplacementByLagrangeMultiplierElement<ELEMENT> *el_pt =
597 dynamic_cast<ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>*
>
598 (Lagrange_multiplier_mesh_pt->element_pt(i));
603 el_pt->set_boundary_shape_geom_object_pt(
607 unsigned nnod=el_pt->nnode();
608 for (
unsigned j=0;j<nnod;j++)
610 Node* nod_pt = el_pt->node_pt(j);
613 if ((nod_pt->is_on_boundary(1))||(nod_pt->is_on_boundary(3)))
617 unsigned n_bulk_value=el_pt->nbulk_value(j);
620 unsigned nval=nod_pt->nvalue();
621 for (
unsigned j=n_bulk_value;j<nval;j++)
638 template<
class ELEMENT>
642 unsigned n_element = Lagrange_multiplier_mesh_pt->nelement();
645 for(
unsigned e=0;e<n_element;e++)
648 delete Lagrange_multiplier_mesh_pt->element_pt(e);
652 Lagrange_multiplier_mesh_pt->flush_element_and_node_storage();
661 template<
class ELEMENT>
674 sprintf(filename,
"%s/soln%i.dat",Doc_info.directory().c_str(),
676 some_file.open(filename);
677 solid_mesh_pt()->output(some_file,n_plot);
682 sprintf(filename,
"%s/lagr%i.dat",Doc_info.directory().c_str(),
684 some_file.open(filename);
688 std::vector<FiniteElement*> el_pt;
689 unsigned nelem=Lagrange_multiplier_mesh_pt->nelement();
690 for (
unsigned e=0;e<nelem;e++)
692 el_pt.push_back(Lagrange_multiplier_mesh_pt->finite_element_pt(e));
695 for (
unsigned e=0;e<nelem;e++)
697 el_pt[e]->output(some_file);
711 int main(
int argc,
char* argv[])
720 MPI_Helpers::init(argc,argv,
false);
725 CommandLineArgs::setup(argc,argv);
732 CommandLineArgs::specify_command_line_flag(
"--nel_1d",&nel_1d);
736 CommandLineArgs::specify_command_line_flag(
"--no_adapt");
739 CommandLineArgs::specify_command_line_flag(
"--block_upper_for_elastic_block");
743 CommandLineArgs::specify_command_line_flag(
"--hypre_for_elastic_blocks");
747 CommandLineArgs::specify_command_line_flag
748 (
"--trilinos_cg_for_lagrange_multiplier_blocks");
753 CommandLineArgs::specify_command_line_flag
754 (
"--diagonal_scaling_for_elastic_blocks");
757 CommandLineArgs::parse_and_assign();
760 CommandLineArgs::doc_specified_flags();
767 problem.doc_solution();
770 unsigned max_adapt=1;
771 if (CommandLineArgs::command_line_flag_has_been_set(
"--no_adapt"))
778 for(
unsigned i=0;i<nstep;i++)
785 problem.newton_solve(max_adapt);
788 problem.doc_solution();
795 problem.solid_mesh_pt()->set_lagrangian_nodal_coordinates();
Function-type-object to compare finite elements based on their x coordinate.
bool operator()(FiniteElement *const &el1_pt, FiniteElement *const &el2_pt) const
Comparison. Is x coordinate of el1_pt less than that of el2_pt?
Problem class for deformation of elastic block by prescribed boundary motion.
void delete_lagrange_multiplier_elements()
Delete elements that enforce prescribed boundary motion by Lagrange multiplilers.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of Lagrange multiplier elements.
void actions_after_newton_solve()
Update function (empty)
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of Lagrange multiplier elements.
void doc_solution()
Doc the solution.
PrescribedBoundaryDisplacementProblem()
Constructor:
ElasticRefineableRectangularQuadMesh< ELEMENT > *& solid_mesh_pt()
Access function for the solid mesh.
void actions_before_newton_solve()
Update function (empty)
void create_lagrange_multiplier_elements()
Create elements that enforce prescribed boundary motion by Lagrange multiplilers.
void position(const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
Parametrised position on object: r(zeta). Evaluated at previous timestep. t=0: current time; t>0: pre...
~WarpedLine()
Empty Destructor.
WarpedLine(const WarpedLine &dummy)
Broken copy constructor.
WarpedLine(const double &l)
Constructor: Specify amplitude of deflection from straight horizontal line.
unsigned ngeom_data() const
How many items of Data does the shape of the object depend on? None.
void position(const Vector< double > &zeta, Vector< double > &r) const
Position vector at Lagrangian coordinate zeta.
void operator=(const WarpedLine &)
Broken assignment operator.
double & ampl()
Access to amplitude.
FaceGeometry()
Constructor – required for more recent versions of gcc.
Pseudo-Elastic Solid element class to overload the block preconditioner methods ndof_types() and get_...
unsigned ndof_types() const
Returns the number of DOF types associated with this element: Twice the number of spatial dimensions ...
PseudoElasticBulkElement()
Constructor.
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Preconditioner * get_diagonal_preconditioner()
Create a matrix-based diagonal preconditioner for subsidiary linear systems.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
double Nu
Poisson's ratio.
WarpedLine Boundary_geom_object(0.0)
GeomObject specifying the shape of the boundary: Initially it's flat.
int main(int argc, char *argv[])
Driver code.