Re: [Getfem-users] (no subject)

[Yves Renard]

Dear Bharat Bhushan,

I understand the problem, but do you really need that information. The best and 
simplest way to compute an integral in getFEM is to use the generic assembly 
language (see http://getfem.org/userdoc/gasm_high.html).
If you define the contour in term of element edges, this is not difficult to 
compute the contour integral with the generic assembly language and you don't 
need to know whose degree of freedom is an enriched one and who is not. 
Moreover, if you need to distinguish between the value on the two sides of the 
crack (but this is only to integrate on the crack itself) you can use the 
Xfem_plus and Xfem_minus facilities (see 
http://getfem.org/userdoc/gasm_high.html#xfem-discontinuity-evaluation-with-mesh-fem-level-set).

Of course, it would be possible to know the exact nature of each degree of 
freedom, but this is not so simple because of the operations such as direct sum 
of finite element spaces.
Moreover, you have to know that the implementation is such that for shape 
function completely cut by the crack, there is two degrees of freedom, one for 
each side of the crack. It differs a bit from the standard Xfem implementation 
where there is a degree of freedom for the continuous fem and a degree of 
freedom for the jump. The reason is that this strategy is far easier to 
generalize to the case of multi-cracks (which is supported by the 
implementation done in GetFEM).

Best regards,

Yves.




Le 02/05/2017 à 07:39, Bharat Bhushan a écrit :
> 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
>
> _______________________________________________
> Getfem-users mailing list
> address@hidden
> https://mail.gna.org/listinfo/getfem-users


-- 

  Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
  Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE
  http://math.univ-lyon1.fr/~renard

---------