Dear users of getfem,I am trying to learn getfem for different applications, making individual steps to explore the possibility of the package with python. I am currently simulating linear elasticity, more precisely a simple 2D cantilever with Dirichlet conditions on the sides and an external force on the top (see the attached step2.pdf for illustration), however my geometry (see the attached .geo file) contains 2 materials. I described my problem specifying the region label when adding the linearized elasticity brick:
# adding the physical equations of the linear elasticity to the model
model.add_fem_variable('u', displacement)
# for the Steel
model.add_initialized_data('lambda_S', clambda_S)
model.add_initialized_data('mu_S', cmu_S)
model.add_filtered_fem_variable('lambdastar_S', clambdastar_S)
model.add_isotropic_linearized_elasticity_brick(mim,'u','lambda_S','mu_S',STEEL_CANT)
# for the Aluminium
model.add_initialized_data('lambda_Al', clambda_Al)
model.add_initialized_data('mu_Al', cmu_Al)
model.add_initialized_data('lambdastar_Al', clambdastar_Al)
model.add_isotropic_linearized_elasticity_brick(mim,'u','lambda_Al','mu_Al',ALUMINIUM_CANT)
I have a result for the displacement, however I am blocked to compute
the von Mises stress as I can't specify a region to the
compute_isotropic_linearized_Von_Mises_or_Tresca method... Is there any
solution to overcome this issue. (by the way, I am using gmsh to create
the mesh, the input and output files are attached as well as my python
script for getfem).
I thank you in advance, best regards, -- Dr. Florian Kolbl Senior Research Officer University of Essex School of Computer Science and Electronic Engineering Brain-Computer Interface Lab Room 3A.532D Telephone 01206 874384 Winvenhoe Park Colchester - Essex CO4 3SQ United Kingdom
Step2.pdf
Description: Adobe PDF document
lc=0.2;
//points of the rectangular cantilever
Point(1)={0,0,0,lc};
Point(2)={5,0,0,lc};
Point(3)={5,0.5,0,lc};
Point(4)={0,0.5,0,lc};
Point(5)={10,0.5,0,lc};
Point(6)={10,0,0,lc};
// lines forming the cantilever
Line(1)={1,2};
Line(2)={2,3};
Line(3)={3,4};
Line(4)={4,1};
Line(5)={3,5};
Line(6)={5,6};
Line(7)={6,2};
//generating the complete cantilever shape
Line Loop(1)={1,2,3,4};
Line Loop(2) = {5, 6, 7, 2};
//generating the cantilever cross-section
Plane Surface(1) = {1};
Plane Surface(2) = {2};
//labeling the boundaries and materials
Physical Line(10)={1,7}; //bottom line is region 10
Physical Line(11)={6}; //right line is region 11
Physical Line(12)={3,5}; //top line is region 12
Physical Line(13)={4}; //left line is region 13
Physical Surface(14) = {1}; //steel material is region 14
Physical Surface(15) = {2}; //steel material is region 15
//Line Loop(15) = {5, 6, 7, 2};
//Plane Surface(16) = {15};
problem_step2.msh
Description: Mesh model
# usefull packages
import numpy as np
import getfem as gf
import time
gf.util_trace_level(3)
#########################
# Constants Declaration #
#########################
# Constants for the Steel
E_S = 2e11 # Young's modulus
nu_S = 0.33 # Poisson's ratio
clambda_S = E_S*nu_S/((1+nu_S)*(1-2*nu_S)) # First Lame coefficient (N/cm^2)
cmu_S = E_S/(2*(1+nu_S)) # Second Lame coefficient (N/cm^2)
clambdastar_S = 2*clambda_S*cmu_S/(clambda_S+2*cmu_S); # Used to compute von Mises stress from displacement
# Constants for the Steel
E_Al = 5.7e10 # Young's modulus
nu_Al = 0.35 # Poisson's ratio
clambda_Al = E_Al*nu_Al/((1+nu_Al)*(1-2*nu_Al)) # First Lame coefficient (N/cm^2)
cmu_Al = E_Al/(2*(1+nu_Al)) # Second Lame coefficient (N/cm^2)
clambdastar_Al = 2*clambda_Al*cmu_Al/(clambda_Al+2*cmu_Al); # Used to compute von Mises stress from displacement
# Constants for the simulation
poly_order = 2 # Order of polynomial basis function on elements
integ_order = 5 # Order of integration
max_res = 1e-9 # Maximal resolution
max_iter = 30 # Maximal number of iterration to convergence
########################################
# Managing the geometry of the problem #
########################################
# importing the mesh
Mesh_file_name = 'problem_step2.msh'
mesh = gf.Mesh('import','gmsh',Mesh_file_name)
elements_degree = mesh.dim() # Degree of the finite element methods
# Recovering specific regions for boundary conditions
# see problem.geo file
BOTTOM_BOUND=10
RIGHT_BOUND=11
TOP_BOUND=12
LEFT_BOUND=13
STEEL_CANT=14
ALUMINIUM_CANT=15
##################################################
# Creating the FEM for the considered quantities #
# and managing the associated methods #
##################################################
# create a MeshFem objects for the quantities involvedin the problem
displacement = gf.MeshFem(mesh,elements_degree) # displacement
von_Mises = gf.MeshFem(mesh,1) # von Mises Stress
# assigning basis functions for the quantities on defined elements
method = 'FEM_PK(%d,%d)' %(elements_degree, poly_order)
displacement.set_fem(gf.Fem(method))
von_Mises.set_classical_discontinuous_fem(elements_degree)
# definition of the integration method
int_meth = 'IM_TRIANGLE(%d)'%integ_order
mim = gf.MeshIm(mesh, gf.Integ(int_meth))
######################
# Creating the model #
######################
model = gf.Model('real')
# adding the physical equations of the linear elasticity to the model
model.add_fem_variable('u', displacement)
# for the Steel
model.add_initialized_data('lambda_S', clambda_S)
model.add_initialized_data('mu_S', cmu_S)
model.add_filtered_fem_variable('lambdastar_S', clambdastar_S)
model.add_isotropic_linearized_elasticity_brick(mim,'u','lambda_S','mu_S',STEEL_CANT)
# for the Aluminium
model.add_initialized_data('lambda_Al', clambda_Al)
model.add_initialized_data('mu_Al', cmu_Al)
model.add_initialized_data('lambdastar_Al', clambdastar_Al)
model.add_isotropic_linearized_elasticity_brick(mim,'u','lambda_Al','mu_Al',ALUMINIUM_CANT)
#######################
# Boundary Conditions #
#######################
# adding a homogeneous Dirichlet condition
model.add_Dirichlet_condition_with_multipliers(mim, 'u', elements_degree-1, LEFT_BOUND)
model.add_Dirichlet_condition_with_multipliers(mim, 'u', elements_degree-1, RIGHT_BOUND)
# adding Neuman add_Dirichlet_condition_with_multipliers
vf = 1 # Vertical force
F = np.zeros(elements_degree)
F[elements_degree-1] = -vf
model.add_initialized_data('Fdata', F)
model.add_source_term_brick(mim, 'u', 'Fdata', TOP_BOUND)
###########
# Solving #
###########
print 'running solve...'
time_start = time.clock()
model.solve('max_res', max_res, 'max_iter', max_iter, 'noisy')
time_elapsed = (time.clock() - time_start)
print "solve done in {} second" .format(time_elapsed)
# visualisation of results
U = model.variable('u')
#VM_S = model.compute_isotropic_linearized_Von_Mises_or_Tresca('u', 'lambdastar', 'mu', von_Mises)
displacement.export_to_vtk('solution_step2.vtk', displacement, U, 'Displacement')_______________________________________________ Getfem-users mailing list [email protected] https://mail.gna.org/listinfo/getfem-users
