Skip to content

Commit e446b75

Browse files
ShiWeiHuHShiWeiHuH
committed
Realization of simple ablation and ablation boundary updates
1 parent 466328a commit e446b75

File tree

8 files changed

+716
-141
lines changed

8 files changed

+716
-141
lines changed

‎src/HyperFSI.jl‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@ using WriteVTK
1515
using Dates
1616

1717
abstract type AbstractFlowTimeSolver end
18+
abstract type AbstractPDGeometry end
1819

1920
# Pre processing
2021
export Post2D, Post3D
@@ -33,9 +34,11 @@ export hsource_bc!, hsource_databc!, temperature_ic!, temperature_bc!, temperatu
3334
export Thermstep, Thermomechstep, Dualstep, Thermstep_ablation, Flowstep, FSI_job, FSI_submit, IBM2D, Bcstruct
3435

3536

37+
include("IBM/boundary_counter.jl")
3638
include("IBM/post_2d.jl")
3739
include("IBM/post_3d.jl")
3840
include("IBM/ibm.jl")
41+
include("IBM/find_new_bc_edges.jl")
3942

4043
include("structure/physics/modify.jl")
4144
include("structure/physics/new_boundaries.jl")

‎src/IBM/boundary_counter.jl‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
2+
mutable struct BoundaryCounter
3+
counts::Dict{Tuple{Vararg{Int}}, Int} # Key without direction → count
4+
records::Dict{Tuple{Vararg{Int}}, Vector{Tuple{Vector{Int},Int}}} # Key without direction → [(face_origional, pd_id)]
5+
end
6+
7+
function BoundaryCounter()
8+
BoundaryCounter(Dict{Tuple{Vararg{Int}}, Int}(),
9+
Dict{Tuple{Vararg{Int}}, Vector{Tuple{Vector{Int},Int}}}())
10+
end
11+
12+
function add_element!(bc::BoundaryCounter, element::Vector{Int}, celltype::String ,pd_id::Int)
13+
for face in extract_faces(element, celltype) # Assuming 2D triangles; modify as needed
14+
key = sort(face) |> Tuple
15+
bc.counts[key] = get(bc.counts, key, 0) + 1
16+
recs = get!(bc.records, key, Vector{Tuple{Vector{Int},Int}}())
17+
push!(recs, (face, pd_id)) # Keep original direction
18+
end
19+
end
20+
21+
function remove_element!(bc::BoundaryCounter, element::Vector{Int}, celltype::String, pd_id::Int)
22+
for face in extract_faces(element, celltype)
23+
key = sort(face) |> Tuple
24+
if !haskey(bc.counts, key)
25+
error("remove_element!: face $face not found")
26+
end
27+
28+
bc.counts[key] -= 1
29+
recs = bc.records[key]
30+
31+
for i in eachindex(recs)
32+
if recs[i][1] == face && recs[i][2] == pd_id
33+
recs[i] = recs[end]
34+
pop!(recs)
35+
break
36+
end
37+
end
38+
39+
if bc.counts[key] == 0
40+
delete!(bc.counts, key)
41+
delete!(bc.records, key)
42+
end
43+
end
44+
end
45+
46+
function get_boundary(bc::BoundaryCounter)
47+
result = Vector{Tuple{Vector{Int},Int}}()
48+
for (key,cnt) in bc.counts
49+
if cnt == 1
50+
append!(result, bc.records[key])
51+
end
52+
end
53+
return result
54+
end
55+
56+
function extract_faces(element::Vector{Int}, celltype::String)
57+
n = length(element)
58+
celltype_upper = uppercase(celltype)
59+
60+
if celltype_upper == "CPS3" && n == 3
61+
return [[element[1],element[2]],
62+
[element[2],element[3]],
63+
[element[3],element[1]]]
64+
65+
elseif celltype_upper == "CPS4" && n == 4
66+
return [[element[1],element[2]],
67+
[element[2],element[3]],
68+
[element[3],element[4]],
69+
[element[4],element[1]]]
70+
71+
elseif celltype_upper == "C3D4" && n == 4
72+
return [[element[1],element[2],element[3]],
73+
[element[1],element[4],element[2]],
74+
[element[2],element[4],element[3]],
75+
[element[3],element[4],element[1]]]
76+
77+
elseif celltype_upper == "C3D8" && n == 8
78+
return [[element[1],element[2],element[3],element[4]],
79+
[element[5],element[6],element[7],element[8]],
80+
[element[1],element[5],element[6],element[2]],
81+
[element[2],element[6],element[7],element[3]],
82+
[element[3],element[7],element[8],element[4]],
83+
[element[4],element[8],element[5],element[1]]]
84+
85+
# 3D（Prism）
86+
elseif (celltype_upper == "C3D6" || celltype_upper == "PRISM") && n == 6
87+
return [
88+
[element[1], element[2], element[3]], # Bottom
89+
[element[4], element[5], element[6]], # Top
90+
[element[1], element[2], element[5], element[4]], # Side 1
91+
[element[2], element[3], element[6], element[5]], # Side 2
92+
[element[3], element[1], element[4], element[6]] # Side 3
93+
]
94+
95+
else
96+
error("Unsupported cell type or element length mismatch: $celltype with $n nodes")
97+
end
98+
end

‎src/IBM/find_new_bc_edges.jl‎

Lines changed: 247 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,247 @@
1+
2+
function update_new_edges!(geo::Post2D, new_ablation_point_idx::Vector{Int})
3+
4+
new_ablation_elements = Vector{Vector{Int}}()
5+
for i in new_ablation_point_idx
6+
mesh_id = geo.pd_id_to_surface[i]
7+
push!(new_ablation_elements, geo.surface_elements[mesh_id])
8+
end
9+
10+
for (i, elem) in enumerate(new_ablation_elements)
11+
pd_id = new_ablation_point_idx[i]
12+
cell_type = geo.element_type[geo.pd_id_to_surface[pd_id]]
13+
remove_element!(geo.bc_counter, elem, cell_type, pd_id)
14+
end
15+
16+
boundary_keys = get_boundary(geo.bc_counter)
17+
18+
new_bc_edges = Dict{Int, Vector{Vector{Vector{Float64}}}}()
19+
20+
for (face, pd_id) in boundary_keys
21+
coords = [geo.mesh_nodes[n] for n in face]
22+
if haskey(new_bc_edges, pd_id)
23+
push!(new_bc_edges[pd_id], coords)
24+
else
25+
new_bc_edges[pd_id] = [coords]
26+
end
27+
end
28+
29+
geo.bc_edge = new_bc_edges
30+
end
31+
32+
function update_new_edges!(geo::Post3D, new_ablation_point_idx::Vector{Int})
33+
34+
new_ablation_elements = Vector{Vector{Int}}()
35+
for i in new_ablation_point_idx
36+
mesh_id = geo.pd_id_to_volume[i]
37+
push!(new_ablation_elements, geo.volume_elements[mesh_id])
38+
end
39+
40+
for (i, elem) in enumerate(new_ablation_elements)
41+
pd_id = new_ablation_point_idx[i]
42+
cell_type = geo.element_type[geo.pd_id_to_volume[pd_id]]
43+
remove_element!(geo.bc_counter, elem, cell_type, pd_id)
44+
end
45+
46+
boundary_keys = get_boundary(geo.bc_counter)
47+
48+
new_bc_surfaces = Dict{Int, Vector{Vector{Vector{Float64}}}}()
49+
50+
for (face, pd_id) in boundary_keys
51+
coords = [geo.mesh_nodes[n] for n in face]
52+
if haskey(new_bc_surfaces, pd_id)
53+
push!(new_bc_surfaces[pd_id], coords)
54+
else
55+
new_bc_surfaces[pd_id] = [coords]
56+
end
57+
end
58+
59+
geo.bc_surface = new_bc_surfaces
60+
end
61+
62+
function save_bc_edges_vtk(geo::Post2D, root::String, step::Int, Δt::Float64; pvd_name="bc_edges.pvd")
63+
64+
mkpath(root)
65+
new_bc_edges = geo.bc_edge
66+
vtkfile = joinpath(root, "bc_edges_" * lpad(string(step), 4, '0') * ".vtk")
67+
68+
all_points = Vector{Vector{Float64}}()
69+
polylines = Vector{Vector{Int}}()
70+
polyline_bc_ids = Vector{Int}()
71+
72+
point_count = 0
73+
for (bc_id, edges) in new_bc_edges
74+
for edge_points in edges
75+
if length(edge_points) >= 2
76+
edge_indices = Vector{Int}()
77+
for point in edge_points
78+
push!(all_points, point)
79+
push!(edge_indices, point_count)
80+
point_count += 1
81+
end
82+
push!(polylines, edge_indices)
83+
push!(polyline_bc_ids, bc_id)
84+
end
85+
end
86+
end
87+
88+
open(vtkfile, "w") do io
89+
println(io, "# vtk DataFile Version 3.0")
90+
println(io, "BC Edges step=$step, time=$(step*Δt)")
91+
println(io, "ASCII")
92+
println(io, "DATASET POLYDATA")
93+
94+
# Points
95+
println(io, "POINTS $(length(all_points)) float")
96+
for point in all_points
97+
if length(point) == 2
98+
println(io, "$(point[1]) $(point[2]) 0.0")
99+
else
100+
println(io, "$(point[1]) $(point[2]) $(point[3])")
101+
end
102+
end
103+
104+
# Lines
105+
total_size = sum(length(poly) + 1 for poly in polylines)
106+
println(io, "LINES $(length(polylines)) $total_size")
107+
for poly in polylines
108+
print(io, "$(length(poly))")
109+
for idx in poly
110+
print(io, " $idx")
111+
end
112+
println(io)
113+
end
114+
115+
# Cell data
116+
println(io, "CELL_DATA $(length(polylines))")
117+
println(io, "SCALARS bc_id int 1")
118+
println(io, "LOOKUP_TABLE default")
119+
for bc_id in polyline_bc_ids
120+
println(io, bc_id)
121+
end
122+
end
123+
124+
pvdfile = joinpath(root, pvd_name)
125+
current_time = step * Δt
126+
127+
if step == 0 || !isfile(pvdfile)
128+
open(pvdfile, "w") do io
129+
println(io, """<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">""")
130+
println(io, " <Collection>")
131+
println(io, " <DataSet timestep=\"" *
132+
string(round(current_time, digits=6)) *
133+
"\" group=\"\" part=\"0\" file=\"" *
134+
basename(vtkfile) * "\"/>")
135+
println(io, " </Collection>")
136+
println(io, "</VTKFile>")
137+
end
138+
else
139+
lines = readlines(pvdfile)
140+
insert!(lines, length(lines)-1,
141+
" <DataSet timestep=\"" *
142+
string(round(current_time, digits=6)) *
143+
"\" group=\"\" part=\"0\" file=\"" *
144+
basename(vtkfile) * "\"/>")
145+
open(pvdfile, "w") do io
146+
for l in lines
147+
println(io, l)
148+
end
149+
end
150+
end
151+
152+
return vtkfile
153+
end
154+
155+
function save_bc_edges_vtk(geo::Post3D, root::String, step::Int, Δt::Float64; pvd_name="bc_surfaces.pvd")
156+
157+
mkpath(root)
158+
new_bc_edges = geo.bc_surface
159+
vtkfile = joinpath(root, "bc_surfaces_" * lpad(string(step), 4, '0') * ".vtk")
160+
161+
all_points = Vector{Vector{Float64}}()
162+
polylines = Vector{Vector{Int}}()
163+
polyline_bc_ids = Vector{Int}()
164+
165+
point_count = 0
166+
for (bc_id, edges) in new_bc_edges
167+
for edge_points in edges
168+
if length(edge_points) >= 2
169+
edge_indices = Vector{Int}()
170+
for point in edge_points
171+
push!(all_points, point)
172+
push!(edge_indices, point_count)
173+
point_count += 1
174+
end
175+
push!(polylines, edge_indices)
176+
push!(polyline_bc_ids, bc_id)
177+
end
178+
end
179+
end
180+
181+
open(vtkfile, "w") do io
182+
println(io, "# vtk DataFile Version 3.0")
183+
println(io, "BC Edges step=$step, time=$(step*Δt)")
184+
println(io, "ASCII")
185+
println(io, "DATASET POLYDATA")
186+
187+
# Points
188+
println(io, "POINTS $(length(all_points)) float")
189+
for point in all_points
190+
if length(point) == 2
191+
println(io, "$(point[1]) $(point[2]) 0.0")
192+
else
193+
println(io, "$(point[1]) $(point[2]) $(point[3])")
194+
end
195+
end
196+
197+
# Lines
198+
total_size = sum(length(poly) + 1 for poly in polylines)
199+
println(io, "LINES $(length(polylines)) $total_size")
200+
for poly in polylines
201+
print(io, "$(length(poly))")
202+
for idx in poly
203+
print(io, " $idx")
204+
end
205+
println(io)
206+
end
207+
208+
# Cell data
209+
println(io, "CELL_DATA $(length(polylines))")
210+
println(io, "SCALARS bc_id int 1")
211+
println(io, "LOOKUP_TABLE default")
212+
for bc_id in polyline_bc_ids
213+
println(io, bc_id)
214+
end
215+
end
216+
217+
pvdfile = joinpath(root, pvd_name)
218+
current_time = step * Δt
219+
220+
if step == 0 || !isfile(pvdfile)
221+
open(pvdfile, "w") do io
222+
println(io, """<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">""")
223+
println(io, " <Collection>")
224+
println(io, " <DataSet timestep=\"" *
225+
string(round(current_time, digits=6)) *
226+
"\" group=\"\" part=\"0\" file=\"" *
227+
basename(vtkfile) * "\"/>")
228+
println(io, " </Collection>")
229+
println(io, "</VTKFile>")
230+
end
231+
else
232+
lines = readlines(pvdfile)
233+
insert!(lines, length(lines)-1,
234+
" <DataSet timestep=\"" *
235+
string(round(current_time, digits=6)) *
236+
"\" group=\"\" part=\"0\" file=\"" *
237+
basename(vtkfile) * "\"/>")
238+
open(pvdfile, "w") do io
239+
for l in lines
240+
println(io, l)
241+
end
242+
end
243+
end
244+
245+
return vtkfile
246+
end
247+

0 commit comments

Comments
 (0)