Detailed documentation to be written. Here's the already fairly well documented driver code...
#include "generic.h"
#include "solid.h"
#include "meshes/simple_cubic_mesh.template.h"
using namespace std;
using namespace oomph;
template<class ELEMENT>
public virtual SolidMesh
{
public:
const double &a, const double &b, const double &c,
TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper) :
SimpleCubicMesh<ELEMENT>(nx,ny,nz,-a,a,-b,b,-c,c,time_stepper_pt)
{
set_lagrangian_nodal_coordinates();
}
};
{
const double& t,
Vector<double>& b)
{
b[0]=0.0;
}
}
template<class ELEMENT>
{
double Shear;
void set_incompressible(ELEMENT *el_pt,const bool &incompressible);
public:
void run(const std::string &dirname);
void doc_solution(DocInfo& doc_info);
void actions_after_newton_solve() {}
void actions_before_newton_solve()
{
apply_boundary_conditions();
bool update_all_solid_nodes=true;
mesh_pt()->node_update(update_all_solid_nodes);
}
void apply_boundary_conditions()
{
unsigned ibound = 5;
unsigned num_nod=mesh_pt()->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
{
SolidNode *solid_nod_pt = static_cast<SolidNode*>(
mesh_pt()->boundary_node_pt(ibound,inod));
solid_nod_pt->x(0) = solid_nod_pt->xi(0) + Shear*
solid_nod_pt->xi(2);
}
}
};
template<class ELEMENT>
: Shear(0.0)
{
double a = 1.0, b = 1.0, c = 1.0;
unsigned nx = 5, ny = 5, nz = 5;
for(unsigned b=0;b<6;b++)
{
unsigned n_node = mesh_pt()->nboundary_node(b);
for(unsigned n=0;n<n_node;n++)
{
for(unsigned i=1;i<3;i++)
{
mesh_pt()->boundary_node_pt(b,n)->pin_position(i);
}
if((b==0) || (b==5))
{
mesh_pt()->boundary_node_pt(b,n)->pin_position(0);
}
}
}
unsigned n_element =mesh_pt()->nelement();
for(unsigned i=0;i<n_element;i++)
{
ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
el_pt->constitutive_law_pt() =
set_incompressible(el_pt,incompressible);
}
cout << assign_eqn_numbers() << std::endl;
}
template<class ELEMENT>
{
ofstream some_file;
char filename[100];
unsigned npts = 5;
sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
doc_info.number());
some_file.open(filename);
mesh_pt()->output(some_file,npts);
some_file.close();
sprintf(filename,"%s/stress%i.dat", doc_info.directory().c_str(),
doc_info.number());
some_file.open(filename);
Vector<double> s(3,0.0);
Vector<double> x(3);
DenseMatrix<double> sigma(3,3);
unsigned n_element = mesh_pt()->nelement();
for(unsigned e=0;e<n_element;e++)
{
ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
el_pt->interpolated_x(s,x);
el_pt->get_stress(s,sigma);
for(unsigned i=0;i<3;i++)
{
some_file << x[i] << " ";
}
for(unsigned i=0;i<3;i++)
{
for(unsigned j=0;j<3;j++)
{
some_file << sigma(i,j) << " ";
}
}
some_file << std::endl;
}
some_file.close();
}
template<class ELEMENT>
{
DocInfo doc_info;
doc_info.set_directory(dirname);
doc_info.number()=0;
unsigned nstep=2;
for(unsigned i=0;i<nstep;i++)
{
newton_solve();
doc_solution(doc_info);
doc_info.number()++;
Shear += 0.5;
}
}
template<>
QPVDElement<3,3> *el_pt, const bool &incompressible)
{
}
template<>
QPVDElementWithPressure<3> *el_pt, const bool &incompressible)
{
else {el_pt->set_compressible();}
}
template<>
set_incompressible(
QPVDElementWithContinuousPressure<3> *el_pt, const bool &incompressible)
{
else {el_pt->set_compressible();}
}
{
for (unsigned i=0;i<2;i++)
{
new IsotropicStrainEnergyFunctionConstitutiveLaw(
{
problem.run("RESLT");
}
{
problem.run("RESLT_pres");
}
{
problem.run("RESLT_cont_pres");
}
}
}
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void doc_solution(DocInfo &doc_info)
Doc the solution.
void set_incompressible(ELEMENT *el_pt, const bool &incompressible)
void run(const std::string &dirname)
Run simulation.
SimpleShearProblem(const bool &incompressible)
Constructor:
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void body_force(const Vector< double > &xi, const double &t, Vector< double > &b)
Body force vector: Vertically downwards with magnitude Gravity.
double Nu
Poisson's ratio.
ConstitutiveLaw * Constitutive_law_pt
Pointer to constitutive law.
double C1
"Mooney Rivlin" coefficient for generalised Mooney Rivlin law
double Gravity
Body force.
StrainEnergyFunction * Strain_energy_function_pt
Pointer to strain energy function.
int main()
Driver for simple elastic problem.