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 > <http://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]; > } > > > > Thank you, > > Tom > > On 16 Oct 2014, at 16:53, Yves Renard <[email protected] > <mailto:[email protected]>> wrote: > >> >> >> 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 <http://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 >> <http://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 >>> [email protected] >>> https://mail.gna.org/listinfo/getfem-users >> >> >> -- >> >> Yves Renard ([email protected]) 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 ([email protected]) 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 ---------
_______________________________________________ Getfem-users mailing list [email protected] https://mail.gna.org/listinfo/getfem-users
