[Getfem-users] (no subject)

[Bharat Bhushan]

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



Regards
Bharat Bhushan
IIT Bombay, Mumbai,India