SHARE
TWEET

Untitled

a guest Jun 12th, 2019 55 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. module test1
  2. import LinearAlgebra.det
  3. import LinearAlgebra.I
  4. eye(n) = Matrix{Float64}(I,n,n)
  5.  
  6.  
  7. struct BulkMesh
  8.     nodepos::Array{Float64, 2}
  9.     eltnodes::Array{Int, 2}
  10. end
  11.  
  12. struct BulkMeshAux
  13.     B_matrix::Vector{Vector{Array{Float64,2}}}
  14. end    
  15.  
  16. function problemsetup0(displacement_delta)
  17.     spacedim = 2
  18.     meshfilename = "testmesh.dat"
  19.     nnode0, nelt0, nodes0, elements0, labels = readmesh(meshfilename)
  20.     nodes0_sample = Vector{Vector{Float64}}()
  21.     nodes0_pin1 = Vector{Vector{Float64}}()
  22.     nodes0_pin3 = Vector{Vector{Float64}}()
  23.     node0lab = zeros(Int, nnode0)
  24.     node0renum = zeros(Int, nnode0)
  25.    
  26.     for j = 1 : nnode0
  27.         push!(nodes0_sample, nodes0[j])
  28.         node0lab[j] = 10
  29.         node0renum[j] = length(nodes0_sample)
  30.  
  31.     end
  32.  
  33.     elements0_sample = Vector{Vector{Int}}()
  34.     elements0_pin1 = Vector{Vector{Int}}()
  35.     elements0_pin3 = Vector{Vector{Int}}()
  36.  
  37.     for j = 1 : nelt0
  38.         erenum = node0renum[elements0[j]]
  39.         push!(elements0_sample, erenum)
  40.     end
  41.    
  42.     bulkmesh0_sample = make_bulkmesh0(nodes0_sample,
  43.                                       elements0_sample)
  44.     nnode0_sample = size(bulkmesh0_sample.nodepos, 1)
  45.     bulkmesh0_pin1 = make_bulkmesh0(nodes0_pin1,
  46.                                     elements0_pin1)
  47.     nnode_pin1 = size(bulkmesh0_pin1.nodepos, 1)
  48.     bulkmesh0_pin3 = make_bulkmesh0(nodes0_pin3,
  49.                                     elements0_pin3)
  50.     bulkmesh_sample = deepcopy(bulkmesh0_sample)
  51.     nnode_sample = size(bulkmesh_sample.nodepos,1)
  52.     allnodes = vcat(bulkmesh_sample.nodepos,
  53.                     bulkmesh0_pin1.nodepos,
  54.                     bulkmesh0_pin3.nodepos)
  55.     allelts = vcat(bulkmesh_sample.eltnodes,
  56.                    bulkmesh0_pin1.eltnodes .+ nnode_sample,
  57.                    bulkmesh0_pin3.eltnodes .+ (nnode_sample + nnode_pin1))
  58.    
  59.     completemesh = BulkMesh(allnodes,
  60.                             allelts)
  61.     bulkmesh_aux = make_bulkmesh_aux(completemesh, triquad3pt())
  62.     numstep = Int(ceil(2.5 / displacement_delta))
  63.     if nnode0_sample == 3
  64.         println("nnode0_sample has changed")
  65.         error("W")
  66.     end
  67.  
  68.     completemesh
  69. end
  70.  
  71. function readmesh(filename)
  72.     nnode0 = 0
  73.     nelt0 = 0
  74.     nodes = Vector{Vector{Float64}}()
  75.     elements = Vector{Vector{Int}}()
  76.     labels = Vector{Int}()
  77.     open(filename, "r") do readhandle
  78.         s = readline(readhandle)
  79.         w = split(s)
  80.         spacedim = parse(Int, w[1])
  81.         nnode_per_elt = parse(Int, w[2])
  82.         s = readline(readhandle)
  83.         w = split(s)
  84.         nnode0 = parse(Int, w[1])
  85.         nelt0 = parse(Int, w[2])
  86.         for nn = 1 : nnode0
  87.             s = readline(readhandle)
  88.             w = split(s)
  89.             push!(nodes, zeros(0))
  90.             for d = 1 : spacedim
  91.                 push!(nodes[nn], parse(Float64, w[d]))
  92.             end
  93.             push!(labels, parse(Int, w[spacedim+1]))
  94.         end
  95.         for ne = 1 : nelt0
  96.             s = readline(readhandle)
  97.             w = split(s)
  98.             push!(elements, zeros(Int,0))
  99.             for i = 1 : nnode_per_elt
  100.                 push!(elements[ne], parse(Int, w[i]))
  101.             end
  102.         end
  103.     end
  104.     return nnode0, nelt0, nodes, elements, labels
  105. end
  106.  
  107. function make_bulkmesh_aux(bulkmesh0, quadrule)
  108.     spacedim = 2
  109.     xyr, polyv1, shapefunc, shapefuncdxi1, shapefuncdxi2 = t6_shapefunc()
  110.     ngau2 = size(quadrule.points,1)
  111.     NIxidxi = Vector{Array{Float64,2}}()
  112.     for s = 1 : ngau2
  113.         qpxi = quadrule.points[s,1]
  114.         qpeta = quadrule.points[s,2]
  115.         monomialvec = [1.0 qpxi  qpeta]
  116.         push!(NIxidxi,[monomialvec * shapefuncdxi1
  117.                        monomialvec * shapefuncdxi2])
  118.     end
  119.     nelement = size(bulkmesh0.eltnodes,1)
  120.     B_matrix = Vector{Vector{Array{Float64,2}}}()
  121.     eye2 = eye(spacedim)
  122.     for e = 1 : nelement
  123.         push!(B_matrix, Vector{Array{Float64,2}}())
  124.         mynodes = bulkmesh0.nodepos[vec(bulkmesh0.eltnodes[e,:]),:]
  125.         for s = 1 : ngau2
  126.             qpxi = quadrule.points[s,1]
  127.             qpeta = quadrule.points[s,2]
  128.             monomialvec = [1 qpxi qpeta]
  129.             dXdxi = [monomialvec * shapefuncdxi1 * mynodes
  130.                      monomialvec * shapefuncdxi2 * mynodes]'
  131.             invdXdxi = inv(dXdxi)
  132.             detdXdxi = det(dXdxi)
  133.             B = invdXdxi' * NIxidxi[s]
  134.             push!(B_matrix[e], B)
  135.         end
  136.     end
  137.     bulkmesh_aux = BulkMeshAux(B_matrix)
  138.     return bulkmesh_aux
  139. end
  140.  
  141.  
  142. function t6_shapefunc()
  143.     xyr = [0.0 0.0
  144.            0.5 0.0
  145.            1.0 0.0
  146.            0.5 0.5
  147.            0.0 1.0
  148.            0.0 0.5]
  149.     polyv1 = [ones(6,1) xyr[:,1] xyr[:,2] xyr[:,1].^2  xyr[:,1].*xyr[:,2] xyr[:,2].^2]
  150.     shapefunc = inv(polyv1)
  151.     shapefuncdxi1 = vcat(vec(shapefunc[2,:])', 2.0*vec(shapefunc[4,:])', vec(shapefunc[5,:])')
  152.     shapefuncdxi2 = vcat(vec(shapefunc[3,:])', vec(shapefunc[5,:])', 2.0*vec(shapefunc[6,:])')
  153.     xyr, polyv1, shapefunc, shapefuncdxi1, shapefuncdxi2
  154. end
  155.  
  156.  
  157. function make_bulkmesh0(nodes0,
  158.                         elements0)
  159.     nnode0 = length(nodes0)
  160.     nelt0 = length(elements0)
  161.     spacedim  = 2
  162.     nnode_per_element = 6
  163.     bulkmesh0 = BulkMesh(zeros(nnode0, spacedim),
  164.                          zeros(Int, nelt0, nnode_per_element))
  165.  
  166.     for nn = 1 : nnode0
  167.         bulkmesh0.nodepos[nn,:] .= nodes0[nn][:]
  168.     end
  169.     for ne = 1 : nelt0
  170.         bulkmesh0.eltnodes[ne,:] .= elements0[ne][:]
  171.     end
  172.     bulkmesh0
  173. end
  174.  
  175. struct QuadRule
  176.     points::Array{Float64,2}
  177.     weights::Vector{Float64}
  178. end
  179.  
  180.  
  181. # 3-point quadrature for T6 elements
  182.  
  183. function triquad3pt()
  184.     QuadRule([1.0/6.0 1.0/6.0; 1.0/6.0 2.0/3.0; 2.0/3.0 1.0/6.0],
  185.              [1.0/6.0,  1.0/6.0,  1.0/6.0])
  186. end
  187.  
  188.  
  189. end
  190.  
  191. for k = 1 : 500
  192.     bulkmesh = test1.problemsetup0(.1)
  193.     println(bulkmesh.nodepos[1,1])
  194. end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top