unstructured_adaptive_solid.cc File Reference

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

Classes
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.

{
 
 //Doc info object
 DocInfo doc_info;
 
 // Output directory
 doc_info.set_directory("RESLT");
 
 // Create generalised Hookean constitutive equations
 Global_Physical_Variables::Constitutive_law_pt = 
  new GeneralisedHookean(&Global_Physical_Variables::Nu);
 
 {
  std::ofstream strain("RESLT/s_energy.dat");
  std::cout << "Running with pure displacement formulation\n";
 
  //Set up the problem
  UnstructuredSolidProblem<ProjectablePVDElement<TPVDElement<2,3> > > problem;
  
  //Output initial configuration
  problem.doc_solution(doc_info);
  doc_info.number()++;
  
  // Parameter study
  Global_Physical_Variables::P=0.0;
  double pressure_increment=0.1e-2;
  
  unsigned nstep=5;
 
  for (unsigned istep=0;istep<nstep;istep++)
   {
    // Solve the problem with one round of adaptivity
    problem.newton_solve(1);
 
    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(doc_info);
    doc_info.number()++;
    
    //Reverse direction of increment 
    if(istep==2) {pressure_increment *= -1.0;}
 
    // Increase (or decrease) load
    Global_Physical_Variables::P+=pressure_increment;
   }
 
  strain.close();
 } //end_displacement_formulation
 
 
 //Repeat for displacement/pressure formulation 
 {
  std::ofstream strain("RESLT_pres_disp/s_energy.dat");
  std::cout << "Running with pressure/displacement formulation\n";
 
  // Change output directory
  doc_info.set_directory("RESLT_pres_disp");
  //Reset doc_info number
  doc_info.number() = 0;
  
  //Set up the problem
  UnstructuredSolidProblem<
  ProjectablePVDElementWithContinuousPressure<
  TPVDElementWithContinuousPressure<2> > > problem;
  
  //Output initial configuration
  problem.doc_solution(doc_info);
  doc_info.number()++;
  
  // Parameter study
  Global_Physical_Variables::P=0.0;
  double pressure_increment=0.1e-2;
  
  unsigned nstep=5;
  for (unsigned istep=0;istep<nstep;istep++)
   {
    // Solve the problem
    problem.newton_solve(1);
 
    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(doc_info);
    doc_info.number()++;
 
    if(istep==2) {pressure_increment *= -1.0;}    
    // Increase (or decrease) pressure load
    Global_Physical_Variables::P+=pressure_increment;
   }
 
  strain.close();
 }
 
 
 //Repeat for displacement/pressure formulation 
 //enforcing incompressibility
 {
  std::ofstream strain("RESLT_pres_disp_incomp/s_energy.dat");
  std::cout << 
   "Running with pressure/displacement formulation (incompressible) \n";
 
  // Change output directory
  doc_info.set_directory("RESLT_pres_disp_incomp");
  //Reset doc_info number
  doc_info.number() = 0;
  
  //Set up the problem
  UnstructuredSolidProblem<
  ProjectablePVDElementWithContinuousPressure<
  TPVDElementWithContinuousPressure<2> > > 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 elastiticy 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(doc_info);
  doc_info.number()++;
  
  // Parameter study
  Global_Physical_Variables::P=0.0;
  double pressure_increment=0.1e-2;
  
  unsigned nstep=5;
  for (unsigned istep=0;istep<nstep;istep++)
   {
    // Solve the problem
    problem.newton_solve(1);
    
    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(doc_info);
    doc_info.number()++;
 
    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, e(), Global_Physical_Variables::Nu, oomph::DocInfo::number(), Global_Physical_Variables::P, problem, oomph::DocInfo::set_directory(), and oomph::PVDEquationsWithPressure< DIM >::set_incompressible().