Dear Xin Liu The boundary condition that you have applied on the left side of your block is not the one you are showing in your figure.
You need to use the function add_normal_Dirichlet_condition_with_multipliers and somehow also restrict your structure in the vertical direction. You can for example apply the same boundary condition to the bottom side. Best regards Kostas On Fri, May 24, 2019 at 3:29 AM Xin Liu <[email protected]> wrote: > Dear GetFEM-Users, > > I wrote a simple 2D uniaxial tensile problem (see the mesh and coordinates > in the attachment) mainly based on the example > "thermo_elasticity_electrical_coupling". > The code is: > " > #include "getfem/getfem_mesh.h" > #include "getfem/getfem_mesh_fem.h" > #include "getfem/getfem_import.h" > #include "getfem/getfem_generic_assembly.h" > #include "gmm/gmm.h" > #include "getfem/getfem_model_solvers.h" > #include "getfem/getfem_mesher.h" > #include "getfem/getfem_fem.h" > > using std::endl; using std::cout; > using bgeot::base_small_vector; > typedef getfem::model_real_plain_vector plain_vector; > > > int main(){ > getfem::mesh m; > getfem::import_mesh("gmsh:./Simple.msh",m); > getfem::mesh_region border_faces; > getfem::outer_faces_of_mesh(m, border_faces); > getfem::mesh_region fb1 = getfem::select_faces_of_normal(m, border_faces, > base_small_vector( 1., 0.), 0.01); > getfem::mesh_region fb2 = getfem::select_faces_of_normal(m, border_faces, > base_small_vector( -1., 0.), 0.01); > m.region(1) = fb1; m.region(2) = fb2; > > // Set finite element methods > getfem::pfem pf = getfem::fem_descriptor("FEM_QK(2,1)"); > getfem::mesh_fem mfu(m); // Finite element for the elastic displacement > mfu.set_qdim(2); > mfu.set_finite_element(pf); > > getfem::mesh_im mim(m); // Integration method > mim.set_integration_method(bgeot::dim_type(2)); > > getfem::model md; > md.add_fem_variable("u", mfu); // Displacement of the structure > > // Membrane elastic deformation > double E = 21E6; > double nu = 0.3; > double F = 100E2; > double clambda = E*nu/((1+nu)*(1-2*nu)); // First Lame coefficient > (N/cm^2) > double cmu = E/(2*(1+nu)); // Second Lame coefficient (N/cm^2) > double clambdastar = 2*clambda*cmu/(clambda+2*cmu); // Lame coefficient > md.add_initialized_scalar_data("cmu", cmu); > md.add_initialized_scalar_data("clambdastar", clambdastar); > getfem::add_isotropic_linearized_elasticity_brick > (md, mim, "u", "clambdastar", "cmu"); > getfem::add_Dirichlet_condition_with_multipliers > (md, mim, "u", bgeot::dim_type(2-1), 2); > md.add_initialized_fixed_size_data("Fdata", base_small_vector(F,0.)); > getfem::add_source_term_brick(md, mim, "u", "Fdata", 1); > gmm::iteration iter(1E-9, 1, 100); > getfem::standard_solve(md, iter); > > // > // Solution export > plain_vector U(mfu.nb_dof()); > gmm::copy(md.real_variable("u"), U); > > bgeot::base_vector coeff; > getfem::slice_vector_on_basic_dof_of_element(mfu,U,1,coeff); > std::vector<std::string> v = { "Point a ", "Point b ", "Point c ", "Point > d "}; > cout << "Displacements for element 1: \n"; > for (int i=0; i < coeff.size()/2; i++) > cout << v[i] << "u1: " << coeff[2*i] << "; u2： " << coeff[2*i+1] << endl; > > return 0; > } > " > > I output the displacement at points a, b, c, and d (see attachment). > Point a u1: 0.000145356; u2： 1.79962e-05 > Point b u1: 0.000161904; u2： 6.42766e-05 > Point c u1: 0.000306614; u2： 2.5169e-05 > Point d u1: 0.000313698; u2： 7.23111e-05 > > I expect the u1 value of a and b is the same and c and d is the same in > this problem since this is just an uniaxial tensile test. So does the > difference come from numerical errors or some mistakes in the code? > Any comments/suggestions will be greatly appreicated. > > > Best regards, > > -- > ====================================== > Xin Liu > PhD Candidate, Graduate Research Assistant > Multiscale Structural Mechanics Group > https://sites.google.com/view/xinliu8809/ > <https://cdmhub.org/members/1776> > School of Aeronautics & Astronautics > Purdue University > West Lafayette, Indiana 47907-2045 > ====================================== >
