/
preprocess.jl
274 lines (233 loc) · 8.47 KB
/
preprocess.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
#=
- read meshes from different formats
- reorder connectivity, create element sets, node sets, ...
- create partitions for parallel runs
- renumber elements / nodes
- maybe precheck for bad elements
- check surface normal direction in boundary elements
- orientation of 2d elements
- etc only topology related stuff
=#
import Base: copy
using JuliaFEM
type Mesh
nodes :: Dict{Int, Vector{Float64}}
node_sets :: Dict{Symbol, Set{Int}}
elements :: Dict{Int, Vector{Int}}
element_types :: Dict{Int, Symbol}
element_codes :: Dict{Int, Symbol}
element_sets :: Dict{Symbol, Set{Int}}
surface_sets :: Dict{Symbol, Vector{Tuple{Int, Symbol}}}
surface_types :: Dict{Symbol, Symbol}
end
function Mesh()
return Mesh(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict())
end
"""
Mesh(m::Dict)
Create new `Mesh` using data `m`. It is assumed that `m` is in format what
`abaqus_read_mesh` in `AbaqusReader.jl` is returning.
"""
function Mesh(m::Dict)
mesh = Mesh()
mesh.nodes = m["nodes"]
mesh.elements = m["elements"]
mesh.element_types = m["element_types"]
mesh.surface_sets = m["surface_sets"]
mesh.surface_types = m["surface_types"]
for (nset_name, node_ids) in m["node_sets"]
mesh.node_sets[Symbol(nset_name)] = Set(node_ids)
end
for (elset_name, element_ids) in m["element_sets"]
mesh.element_sets[Symbol(elset_name)] = Set(element_ids)
end
for (surfset_name, surfaces) in m["surface_sets"]
mesh.surface_sets[Symbol(surfset_name)] = surfaces
end
return mesh
end
function add_node!(mesh::Mesh, nid::Int, ncoords::Vector{Float64})
mesh.nodes[nid] = ncoords
end
function add_nodes!(mesh::Mesh, nodes::Dict{Int, Vector{Float64}})
for (nid, ncoords) in nodes
add_node!(mesh, nid, ncoords)
end
end
function add_node_to_node_set!(mesh::Mesh, set_name, nids...)
if !haskey(mesh.node_sets, set_name)
mesh.node_sets[set_name] = Set{Int}()
end
push!(mesh.node_sets[set_name], nids...)
return
end
""" Create a new node set from nodes in element set. """
function create_node_set_from_element_set!(mesh::Mesh, set_names::String...)
for set_name in set_names
set_name = Symbol(set_name)
info("Creating node set $set_name from element set")
node_ids = Set{Int}()
for elid in mesh.element_sets[set_name]
push!(node_ids, mesh.elements[elid]...)
end
mesh.node_sets[set_name] = node_ids
end
return
end
function create_node_set_from_element_set!(mesh::Mesh, set_name::Symbol)
create_node_set_from_element_set!(mesh, string(set_name))
end
function add_element!(mesh::Mesh, elid::Int, eltype::Symbol, connectivity::Vector{Int})
mesh.elements[elid] = connectivity
mesh.element_types[elid] = eltype
end
function add_elements!(mesh::Mesh, elements::Dict{Int, Tuple{Symbol, Vector{Int}}})
for (elid, (eltype, elcon)) in elements
add_element!(mesh, elid, eltype, elcon)
end
end
function add_element_to_element_set!(mesh::Mesh, set_name, elids...)
if !haskey(mesh.element_sets, set_name)
mesh.element_sets[set_name] = Set{Int}()
end
push!(mesh.element_sets[set_name], elids...)
end
function copy(mesh::Mesh)
mesh2 = Mesh()
mesh2.nodes = copy(mesh.nodes)
mesh2.node_sets = copy(mesh.node_sets)
mesh2.elements = copy(mesh.elements)
mesh2.element_types = copy(mesh.element_types)
mesh2.element_sets = copy(mesh.element_sets)
return mesh2
end
function filter_by_element_id(mesh::Mesh, element_ids::Vector{Int})
mesh2 = copy(mesh)
mesh2.elements = Dict()
for elid in element_ids
if haskey(mesh.elements, elid)
mesh2.elements[elid] = mesh.elements[elid]
end
end
return mesh2
end
function filter_by_element_set(mesh::Mesh, set_name)
filter_by_element_id(mesh::Mesh, collect(mesh.element_sets[set_name]))
end
function create_element(mesh::Mesh, id::Int)
connectivity = mesh.elements[id]
element_type = getfield(JuliaFEM, mesh.element_types[id])
element = Element(element_type, connectivity)
element.id = id
update!(element, "geometry", mesh.nodes)
return element
end
function create_elements(mesh::Mesh; element_type=nothing)
element_ids = collect(keys(mesh.elements))
if element_type != nothing
filter!(id -> mesh.element_types[id] == element_type, element_ids)
end
elements = [create_element(mesh, id) for id in element_ids]
return elements
end
function create_elements(mesh::Mesh, element_sets::Symbol...; element_type=nothing)
if isempty(element_sets)
element_ids = collect(keys(mesh.elements))
else
element_ids = Set{Int}()
for set_name in element_sets
element_ids = union(element_ids, mesh.element_sets[set_name])
end
end
if element_type != nothing
filter!(id -> mesh.element_types[id] == element_type, element_ids)
end
elements = [create_element(mesh, id) for id in element_ids]
return elements
end
function create_elements(mesh::Mesh, element_sets::AbstractString...; element_type=nothing)
element_sets = map(parse, element_sets)
return create_elements(mesh, element_sets...; element_type=element_type)
end
""" find npts nearest nodes from mesh and return id numbers as list. """
function find_nearest_nodes(mesh::Mesh, coords::Vector{Float64}, npts::Int=1; node_set=nothing)
dist = Dict{Int, Float64}()
for (nid, c) in mesh.nodes
if node_set != nothing && !(nid in mesh.node_sets[Symbol(node_set)])
continue
end
dist[nid] = norm(coords-c)
end
s = sort(collect(dist), by=x->x[2])
nd = s[1:npts] # [(id1, dist1), (id2, dist2), ..., (id_npts, dist_npts)]
node_ids = [n[1] for n in nd]
return node_ids
end
function find_nearest_node(mesh::Mesh, coords::Vector{Float64}; node_set=nothing)
return first(find_nearest_nodes(mesh, coords, 1; node_set=node_set))
end
"""
Apply new node ordering to elements. In JuliaFEM same node ordering is used
than in ABAQUS and if mesh is parsed from FEM format with other node ordering
this can be used to reorder nodes.
Parameters
----------
mapping :: Dict{Symbol, Vector{Int}}
e.g. :Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
"""
function reorder_element_connectivity!(mesh::Mesh, mapping::Dict{Symbol, Vector{Int}})
for (elid, eltype) in mesh.element_types
haskey(mapping, eltype) || continue
new_order = mapping[eltype]
element_connectivity = mesh.elements[elid]
new_element_connectivity = element_connectivity[new_order]
mesh.elements[elid] = new_element_connectivity
end
end
function JuliaFEM.Problem{P<:FieldProblem}(mesh::Mesh, ::Type{P}, name::AbstractString, dimension::Int)
problem = Problem(P, name, dimension)
problem.elements = create_elements(mesh, name)
return problem
end
function JuliaFEM.Problem{P<:BoundaryProblem}(mesh::Mesh, ::Type{P}, name, dimension, parent_field_name)
problem = Problem(P, name, dimension, parent_field_name)
problem.elements = create_elements(mesh, name)
return problem
end
using AbaqusReader
function abaqus_read_mesh(fn::String)
m = AbaqusReader.abaqus_read_mesh(fn)
return Mesh(m)
end
"""
create_surface_elements(mesh::Mesh, surface_name::Symbol)
Create a set of surface elements from solid elements.
Notation follow what is defined in ABAQUS. For example, if solid elements
are Tet10, surface elements will be Tri6 and they can be used to define
boundary conditions.
"""
function create_surface_elements(mesh::Mesh, surface_name::Symbol)
elements = Element[]
for (elid, elsi) in mesh.surface_sets[surface_name]
elty = mesh.element_types[elid]
elco = mesh.elements[elid]
chel, chcon = AbaqusReader.create_surface_element(elty, elsi, elco)
ch = Element(getfield(JuliaFEM, chel), chcon)
push!(elements, ch)
end
update!(elements, "geometry", mesh.nodes)
return elements
end
"""
create_surface_elements(mesh::Mesh, surface_name::String)
Create a set of surface elements from solid elements.
Notation follow what is defined in ABAQUS. For example, if solid elements
are Tet10, surface elements will be Tri6 and they can be used to define
boundary conditions.
"""
function create_surface_elements(mesh::Mesh, surface_name::String)
return create_surface_elements(mesh, Symbol(surface_name))
end
export abaqus_read_mesh, create_surface_elements