solid_contact_with_gravity.cc File Reference

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

Classes
class	CircularPenetratorElement
class	ContactProblem< ELEMENT >
	Problem class. More...

Namespaces
	ProblemParameters
	Namespace for exact solution and problem parameters.

Functions
void	ProblemParameters::body_force (const double &time, const Vector< double > &x, Vector< double > &result)
	The body force function. More...
int	main (int argc, char *argv[])
	Driver code. More...

Variables
double	ProblemParameters::Radius_of_elastic_body =2.0
	Radius of elastic body. More...

Function Documentation

main()

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

Driver code.

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

{
 
 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;
 
#ifdef STRUCTURED_MESH
 
 // Build problem
 ContactProblem<RefineableQPVDElement<2,3> > problem;
 
#else
 
 // Build problem
 ContactProblem<ProjectablePVDElement<TPVDElement<2,3> > > problem;
 
#endif
 
 
 // Update position of centre: Y_c has been determined as the 
 // vertical position of a (control) node on the contact boundary.
 // We place the penetrator vertically above this.
 ProblemParameters::Centre[1]=ProblemParameters::Y_c
  +ProblemParameters::Radius;
 
  //Output initial condition
 problem.doc_solution();
 
 unsigned max_adapt=1;
 if (CommandLineArgs::command_line_flag_has_been_set("--no_adapt"))
  {
   max_adapt=0;
  }
 
 
 // Move position of centre directly
 //----------------------------------
 double dyc=0.00024;
 unsigned nstep=1;
 for (unsigned i=0;i<nstep;i++)
  {
   ProblemParameters::Centre[1]-=dyc;
   oomph_info << "Re-solving imposed circle pos for yc=" 
              << ProblemParameters::Centre[1]
              << std::endl;
 
   // Solve
   problem.newton_solve(max_adapt);
   
   //Output solution
   problem.doc_solution();
  }
 
 // Now increase resolution
 //------------------------ 
 nstep=1;
 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();
 dyc=0.0003; // 0.00006; 
 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
   problem.newton_solve(max_adapt);
   
   //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<CircularPenetratorElement*>(
             ProblemParameters::Penetrator_pt)->target_weight() 
            << " and horizontal force: "
            << dynamic_cast<CircularPenetratorElement*>(
             ProblemParameters::Penetrator_pt)->target_horizontal_force() 
            << std::endl;
 
 // Re-solve
 problem.newton_solve(max_adapt);
 
 //Output solution
 problem.doc_solution();
 
 
 // Now increase weight
 //-------------------- 
 double dweight=0.0003; // 0.0001; 
 nstep=1;
 for (unsigned i=0;i<nstep;i++)
  {
   oomph_info << "Re-solving for weight=" 
              << dynamic_cast<CircularPenetratorElement*>(
               ProblemParameters::Penetrator_pt)->target_weight() 
              << " and horizontal force: "
              << dynamic_cast<CircularPenetratorElement*>(
               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 increase body force
//------------------------
 nstep=4;
 for (unsigned i=0;i<nstep;i++)
  {
  // Bump up
  ProblemParameters::Body_force_amplitude+=1.0;
 
  oomph_info << "Re-solving for weight=" 
              << dynamic_cast<CircularPenetratorElement*>(
               ProblemParameters::Penetrator_pt)->target_weight() 
              << " and horizontal force: "
              << dynamic_cast<CircularPenetratorElement*>(
               ProblemParameters::Penetrator_pt)->target_horizontal_force() 
              << " and body force ampl: "
              << ProblemParameters::Body_force_amplitude
              << std::endl;
  
  // Re-solve
  problem.newton_solve(max_adapt);
  
  //Output solution
  problem.doc_solution();
  
 }
 
 exit(0);
 
 // Now increase resolution
 //-------------------- 
 nstep=3;
 for (unsigned i=0;i<nstep;i++)
 {
 
  oomph_info << "Re-solving for weight=" 
             << dynamic_cast<CircularPenetratorElement*>(
              ProblemParameters::Penetrator_pt)->target_weight() 
             << " and horizontal force: "
             << dynamic_cast<CircularPenetratorElement*>(
              ProblemParameters::Penetrator_pt)->target_horizontal_force() 
             << " and body force ampl: "
             << ProblemParameters::Body_force_amplitude
             << " and 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;
 }
 
 
}

References ProblemParameters::Body_force_amplitude, ProblemParameters::Centre, oomph::CommandLineArgs::command_line_flag_has_been_set(), ProblemParameters::Constitutive_law_pt, oomph::CommandLineArgs::doc_specified_flags(), ProblemParameters::El_area, ProblemParameters::Element_length_factor, ProblemParameters::Horizontal_force, i, ProblemParameters::Impose_position_of_centre, ProblemParameters::Nu, oomph::oomph_info, oomph::CommandLineArgs::parse_and_assign(), ProblemParameters::Penetrator_pt, problem, ProblemParameters::Radius, Flag_definition::setup(), oomph::CommandLineArgs::specify_command_line_flag(), ProblemParameters::Weight, and ProblemParameters::Y_c.