[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Getfem-users] Displacement fields in uniaxial tensile test

From: Konstantinos Poulios
Subject: Re: [Getfem-users] Displacement fields in uniaxial tensile test
Date: Fri, 24 May 2019 09:29:35 +0200

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


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

On Fri, May 24, 2019 at 3:29 AM Xin Liu <address@hidden> 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::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

getfem::mesh_im mim(m); // Integration method

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);
(md, mim, "u", "clambdastar", "cmu");
(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;
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
School of Aeronautics & Astronautics
Purdue University
West Lafayette, Indiana 47907-2045

reply via email to

[Prev in Thread] Current Thread [Next in Thread]