adaptive_unstructured_two_d_solid.cc File Reference

#include "generic.h"
#include "solid.h"
#include "constitutive.h"
#include "meshes/triangle_mesh.h"
#include "meshes/quad_from_triangle_mesh.h"

Classes
class	MyMesh< ELEMENT >
	Triangle-based mesh upgraded to become a solid mesh. More...
class	UnstructuredSolidProblem< ELEMENT, MESH >
	Unstructured solid problem. More...

Namespaces
	Global_Physical_Variables
	Global variables.

Functions
void	Global_Physical_Variables::constant_pressure (const Vector< double > &xi, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
	Constant pressure load. More...
int	main (int argc, char **argv)
	Demonstrate how to solve an unstructured solid problem. More...

Function Documentation

main()

int main	(	int argc	,
		char **	argv
	)

Demonstrate how to solve an unstructured solid problem.

{
 // Store command line arguments
 CommandLineArgs::setup(argc,argv);
 
 // Define case to be run
 unsigned i_case=0;
 CommandLineArgs::specify_command_line_flag("--case",&i_case);
 
 // Parse command line
 CommandLineArgs::parse_and_assign(); 
 
 // Doc what has actually been specified on the command line
 CommandLineArgs::doc_specified_flags();
 
 // Number of adaptations
 unsigned max_adapt=2;
 
 // Number of steps
 unsigned nstep=3;
 
 // Create generalised Hookean constitutive equations
 Global_Physical_Variables::Constitutive_law_pt = 
  new GeneralisedHookean(&Global_Physical_Variables::Nu);
 
 // Open an output file for the strain data
 std::ofstream strain("RESLT/s_energy.dat");
 
 if (i_case==0)
 {
  std::cout << "Running with pure displacement formulation\n";
 
  typedef RefineableQPVDElement<2,3> ELEMENT;
  typedef MyMesh<ELEMENT> MESH;
  
  // Set up the problem
  UnstructuredSolidProblem<ELEMENT,MESH> problem;
  
  // Output initial configuration
  problem.doc_solution();
  
  // Parameter study
  Global_Physical_Variables::P=0.0;
  double pressure_increment=0.1e-2;
  
  for (unsigned istep=0;istep<nstep;istep++)
  {
   // Solve the problem with one round of adaptivity
   problem.newton_solve(max_adapt);
 
   double strain_energy = problem.get_strain_energy();
   std::cout << "Strain energy is " << strain_energy << "\n";
    
   // Output strain energy to file
   strain << Global_Physical_Variables::P << " " << strain_energy << std::endl;
 
   // Output solution
   problem.doc_solution();
    
   // On the 3rd run
   if (istep==2)
   {
    pressure_increment*=-1.0;
   }
 
   // Increase (or decrease) load
   Global_Physical_Variables::P+=pressure_increment;
  }
 
  strain.close();
 } //end_displacement_formulation
 else if (i_case==1)
  // Repeat for displacement/pressure formulation 
 {
  std::cout << "Running with pressure/displacement formulation\n";
  
  typedef RefineableQPVDElementWithContinuousPressure<2> ELEMENT;
  typedef MyMesh<ELEMENT> MESH;
  
  // Set up the problem
  UnstructuredSolidProblem<ELEMENT,MESH> problem;
  
  // Output initial configuration
  problem.doc_solution();
  
  // Parameter study
  Global_Physical_Variables::P=0.0;
  double pressure_increment=0.1e-2;
 
  // Loop over the steps
  for (unsigned istep=0;istep<nstep;istep++)
  {
   // Solve the problem
   problem.newton_solve(max_adapt);
 
   double strain_energy = problem.get_strain_energy();
   std::cout << "Strain energy is "<< strain_energy << "\n";
    
   // Output strain energy to file
   strain << Global_Physical_Variables::P << " " << strain_energy << std::endl;
    
   // Output solution
   problem.doc_solution();
 
   // On the 3rd run
   if (istep==2)
   {
    pressure_increment *= -1.0;
   }
    
   // Increase (or decrease) pressure load
   Global_Physical_Variables::P+=pressure_increment;
  }
 
  strain.close();
 }
 else if (i_case==2)
  // Repeat for displacement/pressure formulation enforcing incompressibility
 {
  std::cout << "Running with pressure/displacement "
        << "formulation (incompressible) " << std::endl;
  
  typedef RefineableQPVDElementWithContinuousPressure<2> ELEMENT;
  typedef MyMesh<ELEMENT> MESH;
  
  // Set up the problem
  UnstructuredSolidProblem<ELEMENT,MESH> problem;
  
  // Loop over all elements and set incompressibility flag
  {
   const unsigned n_element = problem.mesh_pt()->nelement();
   for(unsigned e=0;e<n_element;e++)
   {
    // Cast the element to the equation base of our 2D elasticity elements
    PVDEquationsWithPressure<2> *cast_el_pt =
     dynamic_cast<PVDEquationsWithPressure<2>*>(
      problem.mesh_pt()->element_pt(e));
     
    // If the cast was successful, it's a bulk element, 
    // so set the incompressibilty flag
    if (cast_el_pt)
    {
     cast_el_pt->set_incompressible();
    }
   }
  }
 
  // Turn on the incompressibity flag so that elements stay incompressible
  // after refinement
  problem.set_incompressible();
 
  // Output initial configuration
  problem.doc_solution();
  
  // Parameter study
  Global_Physical_Variables::P=0.0;
  double pressure_increment=0.1e-2;
  
  for (unsigned istep=0;istep<nstep;istep++)
  {
   // Solve the problem
   problem.newton_solve(max_adapt);
    
   double strain_energy=problem.get_strain_energy();
   std::cout << "Strain energy is " << strain_energy << "\n";
    
   // Output strain energy to file
   strain << Global_Physical_Variables::P << " " << strain_energy << std::endl;
    
   // Output solution
   problem.doc_solution();
 
   // On the 3rd run
   if (istep==2)
   {
    pressure_increment *= -1.0;
   }
    
   // Increase (or decrease) pressure load
   Global_Physical_Variables::P+=pressure_increment;
  }
 
  strain.close();
 }
 
} // end main

References Global_Physical_Variables::Constitutive_law_pt, oomph::CommandLineArgs::doc_specified_flags(), e(), Global_Physical_Variables::Nu, Global_Physical_Variables::P, oomph::CommandLineArgs::parse_and_assign(), problem, oomph::PVDEquationsWithPressure< DIM >::set_incompressible(), Flag_definition::setup(), and oomph::CommandLineArgs::specify_command_line_flag().