[Getfem-users] Non-homogeneous linear elastic problem with non-constant load

[Tom Haeck]

Hi all,

I am dealing with a very basic non-homogeneous linear elastic problem.  I’ve constructed a tetrahedral mesh of the material.   I need to find the displacements of the nodes of the mesh.  I’ve added two mesh_regions to the mesh.  In one region, the forces on the nodes are provided (and the forces are different for different nodes).  In the other region, the displacements are assumed to be zero. 

In order to tackle this problem in C++, I started off from an example from Daniel Burkhart, posted on the getfem mailing list in 2009 (see below).  This example from Daniel is very similar to mine and I only need to adapt two things:

• In Daniel’s example, the force is constant on the TOP region.  What data structure do I use to parse the non-constant force on the TOP region to the mdbrick_source_term<>?
• In Daniel’s example, the lame coefficients are constant throughout the material. Same question here: what data structure do I use to parse the non-constant Lame coefficients to the mdbrick_isotropic_linearized_elasticity<>?

Thanks!

Tom 

> 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);