Hi everyone,
I have tried to convert the tripod example (from
http://download.gna.org/getfem/doc/getfem_matlab/gfm_12.html) from the
matlab interface to the C++ interface.
The attached Screenshots show input mesh (green region is TOP, boundary),
and the (very noisy) output mesh)
Any suggestion about errors in my code would be great. Maybe the last lines,
where I update the positions of the nodes are incorrect, but I don't know
the correct choice!
...
And here my simulation Code:
Precondition: The tripod mesh is loaded to m_mesh, and the boundaries
(TOP,BOTTOM) are initialized. This is done with a GUI.
std::cout << "MESH:" << std::endl;
std::cout << " Nodes: " << m_mesh.nb_points() << std::endl;
std::cout << " Convex: " << m_mesh.nb_convex() << std::endl;
bool LINEAR = true;
bool INCOMPRESSIBLE = false;
int TOP = 0;
int BOTTOM = 1;
getfem::mesh_fem mfu(m_mesh, 3); //mesh-fem supporting a 3D-vector
field
getfem::mesh_fem mfd(m_mesh, 1); //scalar mesh_fem, for data fields.
//the mesh_im stores the integration methods for each tetrahedron
getfem::mesh_im mim(m_mesh);
mim.set_integration_method(getfem::int_method_descriptor("IM_TETRAHEDRON(5)"));
//we choose a P2 fem for the main unknown
mfu.set_finite_element(getfem::fem_descriptor("FEM_PK(3,2)"));
//the material is homogeneous, hence we use a P0 fem for the data
mfd.set_finite_element(getfem::fem_descriptor("FEM_PK(3,0)"));
//boundary stuff
//ensure top is boundary 0, bottom is boundary 1
//Lamé coefficient
double E = 1000;
double Nu = 0.3;
//set the Lame coefficients
double lambda = E*Nu/((1+Nu)*(1-2*Nu));
double mu = E/(2*(1+Nu));
//create a meshfem for the pressure field (used if incompressible = 0)
getfem::mesh_fem mfp(m_mesh);
mfp.set_finite_element(getfem::fem_descriptor("FEM_PK_DISCONTINUOUS(3,0)"));
//mesh stats:
cout << "Number of dofs for u:" << mfu.nb_dof() << endl;
cout << "Number of dofs for p:" << mfp.nb_dof() << endl;
//bricks stuff
// the linearized elasticity , for small displacements
getfem::mdbrick_abstract<> *b0=0;
getfem::mdbrick_abstract<> *b1=0;
getfem::mdbrick_abstract<> *b2=0;
getfem::mdbrick_abstract<> *b3=0;
//getfem::mdbrick_isotropic_linearized_elasticity<> *b0 = 0;
//getfem::mdbrick_linear_incomp<> *b1 = 0;
if(LINEAR)
{
b0 = new getfem::mdbrick_isotropic_linearized_elasticity<>(mim, mfu,
lambda, mu);
if(INCOMPRESSIBLE)
b1 = new getfem::mdbrick_linear_incomp<>(*b0,mfp);
else
b1=b0;
}
else
{
bgeot::base_vector p(2); p[0] = 0.02; p[1] = 0.5;
if(INCOMPRESSIBLE)
{
getfem::Mooney_Rivlin_hyperelastic_law HypLaw;
b0 = new getfem::mdbrick_nonlinear_elasticity<>(HypLaw, mim,
mfu, p);
b1 = new getfem::mdbrick_nonlinear_incomp<>(*b0,mfp);
}
else
{
getfem::SaintVenant_Kirchhoff_hyperelastic_law HypLaw;
b0 = new getfem::mdbrick_nonlinear_elasticity<>(HypLaw, mim,
mfu, p);
b1 = b0;
}
}
std::cout << "Init boundary conditions..." << std::endl;
//boundary conditions:
//set a vertical force on the top of the tripod
bgeot::base_vector q(3); q[0] = 0.0; q[1] = -10.0; q[2]=0.0;
b2 = new getfem::mdbrick_source_term<>(*b1,mfd, q, TOP);
//atach to the ground
b3 = new getfem::mdbrick_Dirichlet<> (*b2, BOTTOM);
((getfem::mdbrick_Dirichlet<>*)b3)->set_constraints_type(getfem::PENALIZED_CONSTRAINTS);
//solve
std::cout << "Calculation start\n";
getfem::standard_model_state MS(*b3);
getfem::modeling_standard_plain_vector F(mfd.nb_dof()*3);
((getfem::mdbrick_Dirichlet<>*)b3)->rhs().set(mfd,F);
gmm::iteration iter(1e-10, 2, 5000);
getfem::standard_solve(MS, *b3, iter);
// Solution extraction
getfem::modeling_standard_plain_vector U(mfu.nb_dof());
gmm::copy(((getfem::mdbrick_isotropic_linearized_elasticity<>*)b0)->get_solution(MS),
U);
std::cout << "Finished: " << iter.converged() << std::endl;
int j=0;
for (dal::bv_visitor i(m_mesh.points_index()); !i.finished(); ++i)
{
m_mesh.points()[i][0] -= U[j++];
m_mesh.points()[i][1] -= U[j++];
m_mesh.points()[i][2] -= U[j++];
}
_viewer->SetMesh(&m_mesh);
Thanks for your help!
Cheers,
Daniel