heated_linear_solid_contact_with_gravity.cc File Reference

#include <fenv.h>
#include "generic.h"
#include "solid.h"
#include "linear_elasticity.h"
#include "unsteady_heat.h"
#include "meshes/triangle_mesh.h"
#include "meshes/rectangular_quadmesh.h"
#include "contact_elements.h"

Classes
class	oomph::TemplateFreeHeatedPenetratorFluxElementBase
class	oomph::TLinearHeatAndElasticityElement< DIM, NNODE_1D >
	Bulk element which combines linear elasticity and unsteady heat. More...
class	oomph::FaceGeometry< TLinearHeatAndElasticityElement< DIM, NNODE_1D > >
	Face geometry of the customised element. More...
class	oomph::ProjectableLinearHeatAndElasticityElement< LINEAR_HEAT_AND_ELAST_ELEMENT >
	Linear heat and elasticity upgraded to become projectable. More...
class	oomph::FaceGeometry< ProjectableLinearHeatAndElasticityElement< ELEMENT > >
class	oomph::FaceGeometry< FaceGeometry< ProjectableLinearHeatAndElasticityElement< ELEMENT > > >
class	oomph::TemplateFreeHeatedLinearSurfaceContactElementBase
	Template-free base class for linear contact with a heated penetrator. More...
class	oomph::HeatedLinearSurfaceContactElement< ELEMENT >
class	oomph::HeatedPenetratorFluxElement< ELEMENT >
class	FiniteElementComp
class	WarpedLine
	Warped line in 2D space. More...
class	HeatedCircularPenetratorElement
class	ContactProblem< ELEMENT >
	Problem class. More...

Namespaces
	oomph
	DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
	ProblemParameters
	Namespace for exact solution and problem parameters.

Functions
void	ProblemParameters::unit_flux (const double &time, const Vector< double > &x, double &flux)
	hierher temp flux. More...
int	main (int argc, char *argv[])
	Driver code. More...

Variables
double	ProblemParameters::T_contact =0.0
	hierher More...
double	ProblemParameters::X_contact_end_left =0.3
	Left end of contact region (for unstructured mesh only) More...
double	ProblemParameters::X_contact_end_right =0.7
	Right end of contact region (for unstructured mesh only) More...
WarpedLine	ProblemParameters::Boundary_geom_object_left (1.0e-10, 0.0, X_contact_end_left)
WarpedLine	ProblemParameters::Boundary_geom_object_contact (1.0e-10, X_contact_end_left, X_contact_end_right)
WarpedLine	ProblemParameters::Boundary_geom_object_right (1.0e-10, X_contact_end_right, 1.0)
bool	ProblemParameters::Impose_position_of_centre =true
IsotropicElasticityTensor	ProblemParameters::E (Nu)
	The elasticity tensor. More...
double	ProblemParameters::Weight =0.0
	NOTE: WE IMPOSE EITHER THESE ... More...
double	ProblemParameters::Horizontal_force =0.0
	Horizontal force of penetrator. More...
double	ProblemParameters::Y_c =0.0
	... OR THESE... More...
double	ProblemParameters::Rotation_angle =0.0
	Target rotation angle about control node. More...
double	ProblemParameters::El_area =0.02
	Initial/max element area. More...
double	ProblemParameters::Element_length_factor =0.01
	Factor for element length on contact boundary. More...

Function Documentation

main()

int main	(	int argc	,
		char *	argv[]
	)

Driver code.

///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////

{
 //feenableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW);
 
 FiniteElement::Accept_negative_jacobian=true;
 
 // Store command line arguments
 CommandLineArgs::setup(argc,argv);
 
 // Define possible command line arguments and parse the ones that
 // were actually specified
 
 // Suppress adaptation
 CommandLineArgs::specify_command_line_flag("--no_adapt");
    
 // Initial element size
 CommandLineArgs::specify_command_line_flag("--el_area",
                                            &ProblemParameters::El_area);
    
 // Suppress adaptation
 CommandLineArgs::specify_command_line_flag("--validate");
    
 // Parse command line
 CommandLineArgs::parse_and_assign(); 
 
 // Doc what has actually been specified on the command line
 CommandLineArgs::doc_specified_flags();
 
 // Create generalised Hookean constitutive equations
 ProblemParameters::Constitutive_law_pt = 
  new GeneralisedHookean(&ProblemParameters::Nu);
      
 
 // Define centre of penetrator
 ProblemParameters::Centre.resize(2);
 ProblemParameters::Centre[0]=0.5;
 ProblemParameters::Centre[1]=1.0+ProblemParameters::Radius;
 
 // Build penetrator
 ProblemParameters::Penetrator_pt =
  new HeatedCircularPenetrator(&ProblemParameters::Centre,
                               ProblemParameters::Radius);
 
 
#ifdef STRUCTURED_MESH
 
 // Build problem
 ContactProblem<RefineablePseudoSolidNodeUpdateElement<
  RefineableQLinearElasticityElement<2,3>,
  RefineableQPVDElement<2,3> > >
  problem;
 
#else
 
 // Build problem
 ContactProblem
  <ProjectableLinearHeatAndElasticityElement<
   PseudoSolidNodeUpdateElement<TLinearHeatAndElasticityElement<2,3>,
                                TPVDElement<2,3> > > > problem;
 
#endif
 
 
 // Initialise timestep -- also sets the weights for all timesteppers
 // in the problem.
 double dt=0.01;
 problem.initialise_dt(dt);
 problem.assign_initial_values_impulsive();
 
 double max_residuals = 100;
 problem.max_residuals() = max_residuals;
 unsigned max_iterations = 10;
 
  //Output initial condition
 problem.doc_solution();
 
 unsigned max_adapt=1; // hierher
 if (CommandLineArgs::command_line_flag_has_been_set("--no_adapt"))
  {
   max_adapt=0;
  }
 
 // Pure imposed displacement of upper surface -- no contact
 //---------------------------------------------------------
 {
  //Parameter incrementation
  unsigned nstep=2; 
  for(unsigned i=0;i<nstep;i++)
   {
    
    double d_ampl=0.05;
    double ampl=0.0;
 
#ifdef STRUCTURED_MESH
    
    // Increment imposed boundary displacement
    ProblemParameters::Boundary_geom_object.ampl()+=d_ampl;
    ampl=ProblemParameters::Boundary_geom_object.ampl();
 
#else
    
    // Increment imposed boundary displacement
    ProblemParameters::Boundary_geom_object_left.ampl()+=d_ampl;
    ProblemParameters::Boundary_geom_object_contact.ampl()+=d_ampl;
    ProblemParameters::Boundary_geom_object_right.ampl()+=d_ampl;
    ampl=ProblemParameters::Boundary_geom_object_contact.ampl();
        
#endif
    
    oomph_info << "Re-solving for prescr displ amplitude = "
               << ampl << std::endl;
  
    // Solve the problem with Newton's method, allowing
    // up to max_adapt mesh adaptations after every solve.
    bool first=false;
    problem.unsteady_newton_solve(dt,max_adapt,first); 
    
    // Doc solution
    problem.doc_solution();
   }
 }
 
#ifdef STRUCTURED_MESH
 ProblemParameters::Centre[1]=
  ProblemParameters::Boundary_geom_object.y_c()- 
  ProblemParameters::Boundary_geom_object.radius()+
  ProblemParameters::Radius;
#else
 ProblemParameters::Centre[1]=
  ProblemParameters::Boundary_geom_object_contact.y_c()- 
  ProblemParameters::Boundary_geom_object_contact.radius()+
  ProblemParameters::Radius;
#endif
 
 
 // Switch on temperature
 ProblemParameters::T_contact=1.0;
 
 // Move position of centre directly
 //----------------------------------
 unsigned nstep=30;
 for (unsigned i=0;i<nstep;i++)
  {
 
   double dyc=0.00024; 
   ProblemParameters::Centre[1]-=dyc;
   oomph_info << "Re-solving imposed circle pos for yc=" 
              << ProblemParameters::Centre[1]
              << std::endl;
 
   // Solve
   //bool first=false;
   // problem.unsteady_newton_solve(dt,max_adapt,first);
   try
    {      
     max_iterations = 10;
     problem.max_newton_iterations() = max_iterations;
     problem.unsteady_newton_solve(dt);
    }
   catch(OomphLibError& error) //This will catch any error, 
//but actually we only want to catch failure to converge
    {
     max_iterations = 100;
     problem.max_newton_iterations() = max_iterations;
     problem.enable_globally_convergent_newton_method();
     problem.unsteady_newton_solve(dt,false);
     problem.disable_globally_convergent_newton_method();
    }
 
   //Output solution
   problem.doc_solution();
  }
 
 
// exit(0);
 
 // // Now increase resolution
 // //------------------------ 
 // nstep=3;
 // for (unsigned i=0;i<nstep;i++)
 // {
 
 //  oomph_info << "Re-solving imposed circle pos for element length factor: "
 //             << ProblemParameters::Element_length_factor
 //             << std::endl;
 
 //  // Re-solve
 //  problem.newton_solve(max_adapt);
  
 //  //Output solution
 //  problem.doc_solution();
 
 //  // Refine 
 //  ProblemParameters::Element_length_factor/=2.0;
 // }
 
 // Switch to node control
 ProblemParameters::Impose_position_of_centre=false;
 
 // Use displacement control initially
 //-----------------------------------
 problem.adapt();
 problem.switch_to_displ_control();
 double dyc=0.0003; 
 nstep=1;
 for (unsigned i=0;i<nstep;i++)
  {
   ProblemParameters::Y_c-=dyc;
   oomph_info << "Re-solving for yc=" 
              << ProblemParameters::Y_c
              << std::endl;
 
   // Solve
   
   // Try using normal newton solver, if it fails to converge, try a more stabl
   // but slower globally convergent method
   try
    {
     //This problem can have quite high initial residuals
     max_residuals = 100;
 
     //Set max number of iterations to low number
     max_iterations = 10;
     problem.max_newton_iterations() = max_iterations;
 
     problem.newton_solve(max_adapt);
    }
   catch(OomphLibError& error) //This will catch any error, 
       //but actually we only want to catch failure to converge
    {   
     //Increase number of iterations to give it chance to converge
     max_iterations = 100;
     problem.max_newton_iterations() = max_iterations;
 
     //Activate globally convergent method, then deactivate after
     problem.enable_globally_convergent_newton_method();
     problem.newton_solve(max_adapt);
     problem.disable_globally_convergent_newton_method();
    }
 
 
   //Output solution
   problem.doc_solution();
  }
 
 
 // Switch to force control
 problem.switch_to_force_control();
 
 // Balance it in the horizontal direction
 //---------------------------------------
 ProblemParameters::Horizontal_force=0.0;
 
 oomph_info << "RE-solving for weight=" 
            << dynamic_cast<HeatedCircularPenetratorElement*>(
             ProblemParameters::Penetrator_pt)->target_weight() 
            << " and horizontal force: "
            << dynamic_cast<HeatedCircularPenetratorElement*>(
             ProblemParameters::Penetrator_pt)->target_horizontal_force() 
            << std::endl;
 
 // Re-solve
 // problem.unsteady_newton_solve(dt); // max_adapt);
 
 // Try using normal newton solver, if it fails to converge, try a more stable
 // but slower globally convergent method
 try
  {
   //Set max number of iterations to low number
   max_iterations = 10;
   problem.max_newton_iterations() = max_iterations;
   
   problem.unsteady_newton_solve(dt);
  }
 catch(OomphLibError& error) //This will catch any error, 
//but actually we only want to catch failure to converge
  {
   //Increase number of iterations to give it chance to converge
   max_iterations = 100;
   problem.max_newton_iterations() = max_iterations;
   
   //Activate globally convergent method, then deactivate after
   problem.enable_globally_convergent_newton_method();
   problem.unsteady_newton_solve(dt,false);
   problem.disable_globally_convergent_newton_method(); 
  }
 
 
 
 
 //Output solution
 problem.doc_solution();
 
 
 // // Now increase weight
 // //-------------------- 
 // double dweight=0.0001; 
 // nstep=1;
 // for (unsigned i=0;i<nstep;i++)
 //  {
 //   oomph_info << "Re-solving for weight=" 
 //              << dynamic_cast<HeatedCircularPenetratorElement*>(
 //               ProblemParameters::Penetrator_pt)->target_weight() 
 //              << " and horizontal force: "
 //              << dynamic_cast<HeatedCircularPenetratorElement*>(
 //               ProblemParameters::Penetrator_pt)->target_horizontal_force() 
 //              << std::endl;
   
 //   // Re-solve
 //   problem.newton_solve(max_adapt);
   
 //   //Output solution
 //   problem.doc_solution();
   
 //   // Increase weight
 //   ProblemParameters::Weight+=dweight;
 //  }
 
 
// Now detach elastic body
//------------------------
 nstep=3;
 double d_lift_off_ampl=0.0001;
 for (unsigned i=0;i<nstep;i++)
  {
 
   double lift_off=0.0;
 
#ifdef STRUCTURED_MESH
    
    // Increment imposed boundary displacement
    ProblemParameters::Boundary_geom_object.lift_off_amplitude()+=
     d_lift_off_ampl;
 
    lift_off=ProblemParameters::Boundary_geom_object.lift_off_amplitude();
    
#else
    
    // Increment imposed boundary displacement
    ProblemParameters::Boundary_geom_object_left.lift_off_amplitude()+= 
     d_lift_off_ampl;
    ProblemParameters::Boundary_geom_object_contact.lift_off_amplitude()+=
     d_lift_off_ampl;
    ProblemParameters::Boundary_geom_object_right.lift_off_amplitude()+=
     d_lift_off_ampl;
 
    lift_off=ProblemParameters::Boundary_geom_object_contact.
     lift_off_amplitude();
    
#endif
 
 
  oomph_info << "Re-solving for weight=" 
              << dynamic_cast<HeatedCircularPenetratorElement*>(
               ProblemParameters::Penetrator_pt)->target_weight() 
              << " and horizontal force: "
              << dynamic_cast<HeatedCircularPenetratorElement*>(
               ProblemParameters::Penetrator_pt)->target_horizontal_force() 
              << " and lift off: "
              << lift_off
              << std::endl;
  
  // Re-solve
   try
    {
     max_residuals = 100;
     max_iterations = 10;
     problem.max_newton_iterations() = max_iterations;
     problem.max_residuals() = max_residuals;
     problem.newton_solve(max_adapt);
    }
   catch(OomphLibError& error) //This will catch any error, 
//but actually we only want to catch failure to converge
    {
     max_iterations = 100;
     problem.max_newton_iterations() = max_iterations;
     problem.enable_globally_convergent_newton_method();
     problem.newton_solve(max_adapt);
     problem.disable_globally_convergent_newton_method();
    }
 
 
  
  //Output solution
  problem.doc_solution();
  
 }
 
 
 
 
}

References WarpedLine::ampl(), Global_Physical_Variables::Boundary_geom_object, ProblemParameters::Boundary_geom_object_contact, ProblemParameters::Boundary_geom_object_left, ProblemParameters::Boundary_geom_object_right, ProblemParameters::Centre, oomph::CommandLineArgs::command_line_flag_has_been_set(), ProblemParameters::Constitutive_law_pt, oomph::CommandLineArgs::doc_specified_flags(), ProblemParameters::El_area, calibrate::error, ProblemParameters::Horizontal_force, i, ProblemParameters::Impose_position_of_centre, WarpedLine::lift_off_amplitude(), ProblemParameters::Nu, oomph::oomph_info, oomph::CommandLineArgs::parse_and_assign(), ProblemParameters::Penetrator_pt, problem, ProblemParameters::Radius, WarpedLine::radius(), Flag_definition::setup(), oomph::CommandLineArgs::specify_command_line_flag(), ProblemParameters::T_contact, WarpedLine::y_c(), and ProblemParameters::Y_c.