Skip to content

Commit 25929b7

Browse files
boriskausboriskaus
committed
add project_FEData_CartData
1 parent 8f5e20a commit 25929b7

File tree

4 files changed

+108
-51
lines changed

4 files changed

+108
-51
lines changed

‎ext/Gmsh_utils.jl‎

Lines changed: 1 addition & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ module Gmsh_utils
44

55
using GridapGmsh, StaticArrays
66

7-
import GeophysicalModelGenerator: import_Gmsh, FEData
7+
import GeophysicalModelGenerator: import_Gmsh, FEData, CartData
88

99

1010
"""
@@ -66,49 +66,4 @@ end
6666

6767

6868

69-
70-
"""
71-
inside = point_in_tetrahedron(p::_T, a::_T, b::_T, c::_T, d::_T, tol=1e-10)
72-
Determines if a point `p` is inside a tetrahedron specified by `a`,`b`,`c`,`d` or not
73-
"""
74-
function point_in_tetrahedron(p::_T, a::_T, b::_T, c::_T, d::_T, tol=1e-10) where _T<:Vector{Float64}
75-
76-
# check bounding box
77-
xmin = min(a[1],b[1],c[1],d[1])
78-
xmax = max(a[1],b[1],c[1],d[1])
79-
ymin = min(a[2],b[2],c[2],d[2])
80-
ymax = max(a[2],b[2],c[2],d[2])
81-
zmin = min(a[3],b[3],c[3],d[3])
82-
zmax = max(a[3],b[3],c[3],d[3])
83-
84-
inside = true
85-
if p[1] < xmin || p[1] > xmax
86-
inside = false
87-
end
88-
if (p[2] < ymin || p[2] > ymax) && inside
89-
inside = false
90-
end
91-
if (p[3] < zmin || p[3] > zmax) && inside
92-
inside = false
93-
end
94-
95-
if inside
96-
v0 = @SVector [d[i] - a[i] for i in 1:3]
97-
v1 = @SVector [b[i] - a[i] for i in 1:3]
98-
v2 = @SVector [c[i] - a[i] for i in 1:3]
99-
v3 = @SVector [p[i] - a[i] for i in 1:3]
100-
101-
denom = dot(v0, cross(v1, v2))
102-
103-
u = dot(v3, cross(v1, v2)) / denom
104-
v = dot(v0, cross(v3, v2)) / denom
105-
w = dot(v0, cross(v1, v3)) / denom
106-
107-
inside = (u >= -tol) && (v >= -tol) && (w >= -tol) && (u + v + w <= 1 + tol)
108-
end
109-
110-
return inside
111-
end
112-
113-
11469
end

‎src/data_types.jl‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1454,9 +1454,9 @@ function Base.show(io::IO, d::FEData{dim, points_per_cell}) where {dim, points_p
14541454
println(io,"FEData{$dim,$points_per_cell} ")
14551455
println(io," elements : $(size(d.connectivity,2))")
14561456
println(io," vertices : $(size(d.vertices,2))")
1457-
println(io," x ϵ [ $(minimum(d.vertices)[1]) : $(maximum(d.vertices)[1])]")
1458-
println(io," y ϵ [ $(minimum(d.vertices)[2]) : $(maximum(d.vertices)[2])]")
1459-
println(io," z ϵ [ $(minimum(d.vertices)[3]) : $(maximum(d.vertices)[3])]")
1457+
println(io," x ϵ [ $(minimum(d.vertices,dims=2)[1]) : $(maximum(d.vertices,dims=2)[1])]")
1458+
println(io," y ϵ [ $(minimum(d.vertices,dims=2)[2]) : $(maximum(d.vertices,dims=2)[2])]")
1459+
println(io," z ϵ [ $(minimum(d.vertices,dims=2)[3]) : $(maximum(d.vertices,dims=2)[3])]")
14601460
println(io," fields : $(keys(d.fields))")
14611461
println(io," cellfields : $(keys(d.cellfields))")
14621462
end

‎src/transformation.jl‎

Lines changed: 95 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,8 @@
11
# This provides various transformations (GeoData <=> Cartesian; UTMData <=> Cartesian)
22
#
33

4-
export project_CartData
4+
export project_CartData, project_FEData_CartData
5+
56

67

78
"""
@@ -99,4 +100,96 @@ function project_CartData(d_cart::CartData, d::UTMData, p::ProjectionPoint)
99100
end
100101

101102
return d_cart
102-
end
103+
end
104+
105+
106+
107+
"""
108+
inside = point_in_tetrahedron(p::_T, a::_T, b::_T, c::_T, d::_T, tol=1e-10)
109+
Determines if a point `p` is inside a tetrahedron specified by `a`,`b`,`c`,`d` or not
110+
"""
111+
function point_in_tetrahedron(p::_T, a::_T, b::_T, c::_T, d::_T, tol=1e-10) where _T<:Vector{Float64}
112+
113+
# check bounding box
114+
xmin = min(a[1],b[1],c[1],d[1])
115+
xmax = max(a[1],b[1],c[1],d[1])
116+
ymin = min(a[2],b[2],c[2],d[2])
117+
ymax = max(a[2],b[2],c[2],d[2])
118+
zmin = min(a[3],b[3],c[3],d[3])
119+
zmax = max(a[3],b[3],c[3],d[3])
120+
121+
inside = true
122+
if p[1] < xmin || p[1] > xmax
123+
inside = false
124+
end
125+
if (p[2] < ymin || p[2] > ymax) && inside
126+
inside = false
127+
end
128+
if (p[3] < zmin || p[3] > zmax) && inside
129+
inside = false
130+
end
131+
132+
if inside
133+
v0 = @SVector [d[i] - a[i] for i in 1:3]
134+
v1 = @SVector [b[i] - a[i] for i in 1:3]
135+
v2 = @SVector [c[i] - a[i] for i in 1:3]
136+
v3 = @SVector [p[i] - a[i] for i in 1:3]
137+
138+
denom = dot(v0, cross(v1, v2))
139+
140+
u = dot(v3, cross(v1, v2)) / denom
141+
v = dot(v0, cross(v3, v2)) / denom
142+
w = dot(v0, cross(v1, v3)) / denom
143+
144+
inside = (u >= -tol) && (v >= -tol) && (w >= -tol) && (u + v + w <= 1 + tol)
145+
end
146+
147+
return inside
148+
end
149+
150+
"""
151+
data_cart = project_FEData_CartData(data_cart::CartData, data_fe::FEData)
152+
153+
Projects a FEData object with tetrahedrons (e.g., from Gmsh) to a Cartesian grid
154+
"""
155+
function project_FEData_CartData(data_cart::CartData, data_fe::FEData)
156+
157+
cellfields_regions = data_fe.cellfields.regions
158+
regions = zeros(Int64, size(data_cart.x.val))
159+
160+
for i = 1:size(data_fe.connectivity,2) # loop over tetrahedrons
161+
tetra = data_fe.connectivity[:,i]
162+
163+
a = data_fe.vertices[:,tetra[1]]
164+
b = data_fe.vertices[:,tetra[2]]
165+
c = data_fe.vertices[:,tetra[3]]
166+
d = data_fe.vertices[:,tetra[4]]
167+
168+
xmin = min(a[1],b[1],c[1],d[1])
169+
xmax = max(a[1],b[1],c[1],d[1])
170+
ymin = min(a[2],b[2],c[2],d[2])
171+
ymax = max(a[2],b[2],c[2],d[2])
172+
zmin = min(a[3],b[3],c[3],d[3])
173+
zmax = max(a[3],b[3],c[3],d[3])
174+
175+
ind = findall( data_cart.x.val .>= xmin .&& data_cart.x.val .<= xmax .&&
176+
data_cart.y.val .>= ymin .&& data_cart.y.val .<= ymax .&&
177+
data_cart.z.val .>= zmin .&& data_cart.z.val .<= zmax);
178+
179+
for I in ind
180+
x = data_cart.x.val[I]
181+
y = data_cart.y.val[I]
182+
z = data_cart.z.val[I]
183+
p = [x,y,z]
184+
if point_in_tetrahedron(p,a,b,c,d)
185+
regions[I] = cellfields_regions[i]
186+
end
187+
end
188+
189+
end
190+
191+
return addfield(data_cart, (regions=regions,))
192+
end
193+
194+
195+

‎test/test_Gmsh.jl‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,4 +7,13 @@ fname="test_files/subduction_ptatin.msh"
77
fe_data, tag_names = import_Gmsh(fname)
88
@test sum(fe_data.cellfields.regions) == 830
99

10+
# Swap y and z dimensions (as pTatin uses a differenty definition)
11+
data_fe = swap_yz_dims(fe_data)
1012

13+
# Define a CartData set with the same dimensions as the Gmsh file
14+
bbox = extrema(fe_data_gmg);
15+
nx,ny,nz = 100,50,80
16+
data_cart = CartData( xyz_grid(range(bbox[1]...,length=nx),range(bbox[2]...,length=ny),range(bbox[3]...,length=nz) ))
17+
18+
data_cart1 = project_FEData_CartData(data_cart, data_fe)
19+
@test extrema(data_cart1.fields.regions) == (2, 11)

0 commit comments

Comments
 (0)