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 > [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 ---------
_______________________________________________ Getfem-users mailing list [email protected] https://mail.gna.org/listinfo/getfem-users
