|
Dear Tom,
First, the program from Daniel Burkhart uses the old Getfem
bricks. You would better use the new ones and take example on the
program /tests/elastostatic.cc for instance (however, the simpler
is probably to use the Python interface, see examples in
interface/tests/python).
In order to define a non-constant source term of coefficient, you
can either interpolate it on a Lagrange fem with
getfem::interpolation_function function (see
/tests/elastostatic.cc) or you can just add the term on a
particular region (the new bricks accept a region number where to
add the term) or use mf.basic_dof_on_region method
method of the mesh_fem object. Once a non constant coefficient is
defined (i.e. the vector of its dof is computed on a particular
finite element method) it can be given to the linear elsaticity
brick.
Yves.
Le 15/10/2014 17:19, Tom Haeck a écrit :
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);
_______________________________________________
Getfem-users mailing list
address@hidden
https://mail.gna.org/listinfo/getfem-users
--
Yves Renard (address@hidden) tel : (33) 04.72.43.87.08
Pole de Mathematiques, INSA-Lyon fax : (33) 04.72.43.85.29
20, rue Albert Einstein
69621 Villeurbanne Cedex, FRANCE
http://math.univ-lyon1.fr/~renard
---------
|