Dear Daniel, I think this is because of your last loop. You move the point of your mesh with U, but U is the value of dofs which have their own numbering (moreover you use a P2 finite element). So you have first to establish the correspondance between the node indices and the dof indices. You can use for that the function mf.ind_dof_of_element(...) or identify the point using mf.point_of_dof(...) or use an interpolation on a P1 finite element if you want to make a generic function ...
Yves. Daniel Burkhart <[email protected]> a écrit : > 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); > > > > Thanks for your help! > Cheers, > Daniel > _______________________________________________ Getfem-users mailing list [email protected] https://mail.gna.org/listinfo/getfem-users
