GETFEM中xfem程序节点自由度的提取++

2024-03-28 13:02:02 发布

您现在位置:Python中文网/ 问答频道 /正文

但是我需要用有限元方法来解决这个问题,这样我就可以用有限元方法来解决。我使用的是python接口,getfem++的python接口中没有直接计算应力强度因子的函数,所以我必须编写一个计算交互积分的代码。为了编写代码,我必须提取节点自由度的值。但是在提取的时候,我只能得到一个列向量,这个列向量包含了一个特定节点的所有自由度。所以我不知道列矩阵的哪一个元素是哪一种丰富的自由度(它是非丰富的,重边的还是反符号丰富的)。请建议如何提取节点自由度。非常感谢。你知道吗

下面是我的getfem++python代码:

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
CT_Node=[] # crack tip nodes inside the specified radius rad
for i in range(nnodes):
    if ((nodes[i,0]-cx)**2+(nodes[i,1]-cy)**2)<rad**2:
        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
mls.add(ls)
mls.adapt()

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')
md.add_fem_variable('u', mfu)

# Add the lame's coefficients to the model
Lambda = 1.25E10 # Lame coefficient #
Mu = 1.875E10    # Lame coefficient #
md.add_initialized_data('lambda', [Lambda])
md.add_initialized_data('mu', [Mu])

# Add an isotropic elasticity matrix to the model
md.add_isotropic_linearized_elasticity_brick(mim,'u','lambda','mu')

# Apply the Neumann condition.
md.add_initialized_fem_data('NeumannData', mfrhs, F2)
md.add_normal_source_term_brick(mim, 'u', 'NeumannData',top)

# 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)
md.add_initialized_data('cpoints',cpoints)
md.add_initialized_data('cunitv',cunitv)
md.add_pointwise_constraints_with_multipliers('u','cpoints','cunitv')


# 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)

Tags: andtheinfromnodeifnpnodes