getfem-users
[Top][All Lists]

## [Getfem-users] (no subject)

 From: Bharat Bhushan Subject: [Getfem-users] (no subject) Date: Tue, 2 May 2017 11:09:08 +0530

```Dear sir,

I have solved an edge crack problem using GETFEM++ with XFEM
enrichments(code is attached) but i need to know the stress intensity
factor so that i can compare it with the analytical solution. I am
using the python interface and there is no function to directly
calculate the stress intensity factor in python interface of the
getfem++ so i have to write a code to calculate interaction integral.
To write the code i have to extract the values of nodal degrees of
freedom. But when extracting, i can only get a column vector which
contains all the degrees of freedom of a particular node in random
order. So i don't know which element of the column matrix is which
kind of enriched degree of freedom ( whether it is nonenriched,
heavyside or assymptotic enriched). please suggest how can i extract
the nodal degrees of freedom. Thank you.

Here is my getfem++ python code:

from __future__ import division
import getfem as gf
import numpy as np
import matplotlib.pyplot as plt

w=178.0 # Width of cracked plate
l=407.0  # Length of cracked plate

######################################
## Build the Mesh ####################
######################################

# bulding cartesian mesh with 24*48 elements ( 4 - noded quadrilaterals)
m=gf.Mesh('cartesian',np.linspace(0.0,w,25),np.linspace(0.0,l,49))

# extracting the nodal and element connectivity matrices from the mesh
nodes=m.pts().transpose() # Nodal coordinates matrix
elem=m.pid_from_cvid()[0].reshape(m.nbcvs(),4) # Element connectivity matrix
nnodes=len(nodes)
nelem=len(elem)
elem1=elem+1 # adding one to connectivity matrix to later import
in MATLAB for SIF calculation

# saving the nodes and elem matrix in csv format
np.savetxt('nodes.csv',nodes,delimiter=',')
np.savetxt('elem.csv',elem1,delimiter=',')

##############################################################
# detect the nodes to be enriched by crack-tip enrichments ###
##############################################################

# define the Crack tip locations
cx=85.29 # x-coordinate of crack-tip
cy=207.7395 # y-coordinate of crack-tip
rad=8  # Radius of enrichment, taken close to element length so as
to enrich only one layer of nodes
for i in range(nnodes):
CT_Node.append(i)

##############################################
### XFEM-Enrichments #########################
##############################################

# 4-crack tip enrichment functions defined in GetFem++
ck0 = gf.GlobalFunction('crack',0)
ck1 = gf.GlobalFunction('crack',1)
ck2 = gf.GlobalFunction('crack',2)
ck3 = gf.GlobalFunction('crack',3)

# Levelset definition:
ls = gf.LevelSet(m,1,'y-207.7395','x-85.29') # m is the Mesh, no.
of crack=1, equations of crack are y-cy=0 and x-cx=0

mls = gf.MeshLevelSet(m)
# mls is the meshlevelset which will be used for defining the crack

mf_pre_u = gf.MeshFem(m)
mf_pre_u.set_fem(gf.Fem('FEM_QK(2,1)'))

#mf_pre_u is the MeshFem object without crack, with Q4 elements

mfls_u = gf.MeshFem('levelset',mls,mf_pre_u)

# mfls_u is the MeshFem object with the Heavyside enrichments

################################################################
#### define MeshFem with crack-tip enrichments #################
###############################################################

mf_sing_u = gf.MeshFem('global function',m,ls,[ck0,ck1,ck2,ck3],1)
mf_part_unity = gf.MeshFem(m)
mf_part_unity.set_classical_fem(1)
DOFpts = mf_part_unity.basic_dof_nodes()
Idofs_center=np.array(CT_Node)
mf_xfem_sing = gf.MeshFem('product', mf_part_unity, mf_sing_u)
# mf_xfem_sing is the MeshFem object with crack-tip enrichments
mf_xfem_sing.set_enriched_dofs(Idofs_center)
mfu = gf.MeshFem('sum', mfls_u,mf_xfem_sing) # summing up the
MeshFems with heavyside and crack-tip
mfu.set_qdim(2) # set dimenstion eqaual to 2

####################################################################
#### Integration Method ###########################################
####################################################################

mim = gf.MeshIm('levelset', mls, 'all',
gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),9)'),

gf.Integ('IM_STRUCTURED_COMPOSITE(IM_GAUSS_PARALLELEPIPED(2,6),9)'),gf.Integ('IM_STRUCTURED_COMPOSITE(IM_GAUSS_PARALLELEPIPED(2,6),9)'))

# Boundary selection
flst  = m.outer_faces()
fnor  = m.normal_of_faces(flst)
tbot = abs(fnor[1,:]+1) < 1e-14
ttop  = abs(fnor[1,:]-1) < 1e-14
fbot = np.compress(tbot, flst, axis=1)
ftop  = np.compress(ttop, flst, axis=1)
fneum = np.compress(True - ttop - tbot, flst, axis=1)

# Mark it as boundary
bottom = 1
top = 2
N = 3
m.set_region(bottom, fbot)
m.set_region(top, ftop)
m.set_region(N, fneum)

# define the load to be applied on the top boundary
mfrhs = gf.MeshFem(m, 2)
mfrhs.set_fem(gf.Fem('FEM_QK(2,1)'))
F2 = mfrhs.eval('[0,0,0,1000000]')

# Define the Model
md = gf.Model('real')

# Add the lame's coefficients to the model
Lambda = 1.25E10 # Lame coefficient #
Mu = 1.875E10    # Lame coefficient #

# Add an isotropic elasticity matrix to the model

# Apply the Neumann condition.

# Apply the Dirichlet condition so that its simply supported at the bottom

cpoints=[]
cunitv=[]
for i in range(nnodes):
if nodes[i,1]<1e-14:
cpoints.append(nodes[i,0])
cpoints.append(nodes[i,1])
cunitv.append(0)
cunitv.append(1)
cpoints.append(0)
cpoints.append(0)
cunitv.append(1)
cunitv.append(0)

# Assembly of the linear system and solve.
gf.memstats()
md.solve()
U = md.variable('u')

#####################################################
## Export data to Matlab to find SIF's ##############
####################################################

ele_cut=np.array([],dtype=int) #Elements cut by the crack
ele_cut0=np.array([],dtype=int) #Elements cut by the crack and
enriched by only heavyside
ele_cut1=np.array([],dtype=int) #Elements enriched by all
heavyside and one cracktip
ele_cut2= np.array([],dtype=int) #Elements enriched by all
heavyside and two cracktip
ele_cut3= np.array([],dtype=int) #Elements enriched by no
heavyside and all cracktip and all heavyside+3 crack
ele_cut4= np.array([],dtype=int) #Elements enriched by one
heavyside and all cracktip
ele_cut5= np.array([],dtype=int) #Elements enriched by two
heavyside and all cracktip
helem= np.array([],dtype=int)
sumdof=np.zeros([nelem,1],dtype=int)

h_node=[]
nn=np.zeros((nnodes),dtype=int)

for i in range(nnodes):
if (cy-nodes[i,1])>0:
nn[i]=-1
if (cy-nodes[i,1])<0:
nn[i]=1

for i in range(nelem):
if nn[elem[i,0]]+nn[elem[i,1]]+nn[elem[i,2]]+nn[elem[i,3]]==0
and nodes[elem[i,1],0]<cx-8:
h_node=np.append(h_node,i)
if nn[elem[i,0]]+nn[elem[i,1]]+nn[elem[i,2]]+nn[elem[i,3]]==0
and nodes[elem[i,0],0]<cx:
ele_cut=np.append(ele_cut,i)
if len(mfu.basic_dof_from_cv(i))==16:
ele_cut0=np.append(ele_cut0,i)
if len(mfu.basic_dof_from_cv(i))==24:
ele_cut1=np.append(ele_cut1,i)
if len(mfu.basic_dof_from_cv(i))==32:
ele_cut2=np.append(ele_cut2,i)
if len(mfu.basic_dof_from_cv(i))==40:
ele_cut3=np.append(ele_cut3,i)
if len(mfu.basic_dof_from_cv(i))==42:
ele_cut4=np.append(ele_cut4,i)
if len(mfu.basic_dof_from_cv(i))==44:
ele_cut5=np.append(ele_cut5,i)

sumdof[i]=len(mfu.basic_dof_from_cv(i))

#####################################################################3

Node_DOF=np.zeros((nnodes,2),dtype=int)

for i in range(nnodes):
Node_DOF[i,0]=2
if i in CT_Node:
Node_DOF[i,0]=Node_DOF[i,0]+8
if i in np.array(m.pid_from_cvid(h_node)).tolist()[0].tolist():
Node_DOF[i,0]=Node_DOF[i,0]+2
if i!=0:
Node_DOF[i,1]=Node_DOF[i,0]+Node_DOF[i-1,1]
if i==0:
Node_DOF[i,1]=Node_DOF[i,0]

#####################################################################
Ele_Node_DOF=np.zeros((nelem,5),dtype=int)
aa=[0,1,2,3]
for i in range(nelem):

Ele_Node_DOF[i,aa]=2
for j in range(len(elem[i])):
if elem[i][j] in CT_Node:
Ele_Node_DOF[i,j]=Ele_Node_DOF[i,j]+8
if elem[i][j] in
np.array(m.pid_from_cvid(h_node)).tolist()[0].tolist():
Ele_Node_DOF[i,j]=Ele_Node_DOF[i,j]+2
if i!=0:
Ele_Node_DOF[i,4]=np.sum(Ele_Node_DOF[i,aa])+Ele_Node_DOF[i-1,4]
if i==0:
Ele_Node_DOF[i,4]=np.sum(Ele_Node_DOF[i,aa])

H_Node=np.array([],dtype=int)
H_node=np.array(m.pid_from_cvid(h_node)).tolist()[0].tolist()
for i in H_node:
if i not in H_Node:
H_Node=np.append(H_Node,i)

Regards
Bharat Bhushan
IIT Bombay, Mumbai,India

```