Dear Tom,
You can iterate on the mesh faces as you do and obtain the dof
indices on the face with the function
mf.ind_basic_dof_of_face_of_element(cv(), f()).
Of course, this works only for a Lagrange fem.
You ca also use a normal_source_term brick if the force is normal
one. A third mean is to use the generic assembly (the unit normal
vector is available there, see
http://download.gna.org/getfem/html/homepage/userdoc/gasm_high.html).
An _expression_ such as "sin(x[0])*Normal" is allowed to define a
source term in the generic assembly.
Yves.
Le 17/10/2014 17:40, Tom Haeck a écrit :
Dear Yves,
Thank you for your fast response. tests/ elastostatic.cc
is indeed a nice example of how the linear elastic brick works.
However, I encounter troubles when I want to define the
source forces perpendicular to the a certain region (face) and
write it to a plain_vector in order to parse it to a source term
brick together with a getfem::mesh_fem. The problem is that the
numbering of the nodes of the mesh is different from the
numbering of the dofs of the getfem::mesh_fem. With the
numbering of the nodes of the mesh, it is easy to write the
normal vectors in the right place in the plain_vector (see
below). How could I get the same plain_vector, but with the
numbering of the dofs of the getfem::mesh_fem?
plain_vector normal_forces;
for (getfem::mr_visitor i(mesh.region(1)); !i.finished(); ++i) {
bgeot::base_node dir =
mesh.normal_of_face_of_convex(i.cv(), i.f());
for (int j =
0; j < mesh.dim(); j++){
for (int k =
0; k < mesh.dim(); k++)
normal_forces[mesh.ind_points_of_face_of_convex(i.cv(),
i.f())[j]*mesh.dim() + k] += dir[k];
}
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
---------
--
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
---------
|