getfem-users
[Top][All Lists]

## Re: [Getfem-users] Basic example problem

 From: Renard Yves Subject: Re: [Getfem-users] Basic example problem Date: Tue, 24 Mar 2009 20:55:30 +0100 User-agent: Dynamic Internet Messaging Program (DIMP) H3 (1.1)

```
Dear Daniel,

```
I think this is because of your last loop. You move the point of your mesh with U, but U is the value of dofs which have their own numbering (moreover you use a P2 finite element). So you have first to establish the correspondance between the node indices and the dof indices. You can use for that the function mf.ind_dof_of_element(...) or identify the point using mf.point_of_dof(...) or use an interpolation on a P1 finite element if you want to make a generic function ...
```
Yves.

Daniel Burkhart <address@hidden> a écrit :

```
```Hi everyone,
I have tried to convert the tripod example (from
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);

Cheers,
Daniel

```
```

```