• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

a guest Jun 12th, 2019 57 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)
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.
72.     nnode0 = 0
73.     nelt0 = 0
74.     nodes = Vector{Vector{Float64}}()
75.     elements = Vector{Vector{Int}}()
76.     labels = Vector{Int}()
79.         w = split(s)
80.         spacedim = parse(Int, w[1])
81.         nnode_per_elt = parse(Int, w[2])
83.         w = split(s)
84.         nnode0 = parse(Int, w[1])
85.         nelt0 = parse(Int, w[2])
86.         for nn = 1 : nnode0
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
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.
108.     spacedim = 2
109.     xyr, polyv1, shapefunc, shapefuncdxi1, shapefuncdxi2 = t6_shapefunc()
111.     NIxidxi = Vector{Array{Float64,2}}()
112.     for s = 1 : ngau2
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
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.
176.     points::Array{Float64,2}
177.     weights::Vector{Float64}
178. end
179.
180.
181. # 3-point quadrature for T6 elements
182.
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.

Top