Dear all,
I was wondering whether there an existing package that allows to generate shapes such as rhombic dodecahedron.
I am trying to make a 3D model (Cartesian regular mesh) in which I would like to identify cell centroids as either inside or outside the volume of a rhombic dodecahedron.
So far, I could only find this package (Oscar.jl) but I’m not sure if I can use it for my purpose.
Thanks in advance for the hints .
That’s what I ended up doing in the end. Just in case it can be useful to anyone:
using Plots, WriteVTK, LinearAlgebra
```Set of vertices and vertex 2 face connectivity ```
function RhombicDodecahedronGeometry( )
nfac = 12
p = zeros(14,3) # vertices
f2v = zeros(Int64,nfac,3) # face 2 vertices
# Vertices set 1
p[1,:] = [1,1,1]
p[2,:] = [1,1,−1]
p[3,:] = [1,−1,1]
p[4,:] = [1,−1,−1]
p[5,:] = [−1,1,1]
p[6,:] = [−1,1,−1]
p[7,:] = [−1,−1,1]
p[8,:] = [−1,−1,−1]
# Vertices set 2
p[9,:] = [0,0,2]
p[10,:] = [0,0,−2]
p[11,:] = [0,2,0]
p[12,:] = [0,−2,0]
p[13,:] = [2,0,0]
p[14,:] = [−2,0,0]
# Face connectivity
f2v[ 1,:] = [14, 11, 5]
f2v[ 2,:] = [9, 14, 5]
f2v[ 3,:] = [7, 9, 12]
f2v[ 4,:] = [13, 9, 1 ]
f2v[ 5,:] = [11, 13, 1 ]
f2v[ 6,:] = [11, 1, 5 ]
f2v[ 7,:] = [11, 10, 2 ]
f2v[ 8,:] = [10, 13, 2 ]
f2v[ 9,:] = [10, 6, 14 ]
f2v[10,:] = [10, 8, 12 ]
f2v[11,:] = [4, 12, 13 ]
f2v[12,:] = [14, 12, 8 ]
return p, f2v
end
```Compute fit```
function FitRhombicDocdecahedron( p, f2v )
nfac = size(f2v,1)
f = zeros(nfac,4) # fit for each face
for ifac=1:nfac
P1, P2, P3 = p[f2v[ifac,1],:], p[f2v[ifac,2],:], p[f2v[ifac,3],:]
v1 = P2 .- P1
v2 = P2 .- P3
X = cross(v1,v2)
k = -P1'*X
f[ifac,1], f[ifac,2], f[ifac,3], f[ifac,4] = X[1], X[2], X[3], k
end
return f
end
```Identify if a point is inside the rhombic dodecahedron```
function IsInRhombicDodecahedron( X, Y, Z, f )
nfac = size(f,1)
F = true
for ifac=1:nfac
F = F && (f[ifac,1]*X + f[ifac,2]*Y + f[ifac,3]*Z + f[ifac,4]) < 0.0
end
return F
end
function main( n )
# Domain
xmin, xmax = -2.5e-2, 2.5e-2
ymin, ymax = -2.5e-2, 2.5e-2
zmin, zmax = -2.5e-2, 2.5e-2
# Centroid location and embedding radius
x0, y0, z0 = 1.e-3, -2e-3, 3.0e-3
r = 2.0e-2
# Spatial discretisation
nx, ny, nz = n*8+1, n*8+1, n*8+1
dx = (xmax-xmin)/nx
dy = (ymax-ymin)/ny
dz = (zmax-zmin)/nz
xc = LinRange(xmin+dx/2, xmax-dx/2, nx)
yc = LinRange(ymin+dy/2, ymax-dy/2, ny)
zc = LinRange(zmin+dz/2, zmax-dz/2, nz)
# Arrays
phase = zeros(Float64,nx,ny,nz) # phase type
# Geometric data for rhombic dodecahedron
p, f2v = RhombicDodecahedronGeometry( )
# Scale to desired radius
p .*= r/2.0
# Call fitting function
f = FitRhombicDocdecahedron( p, f2v )
# Loop through cells and check if in or out
@Threads.threads for k=1:nz
@inbounds for j=1:ny
for i=1:nx
if IsInRhombicDodecahedron( xc[i]-x0, yc[j]-y0, zc[k]-z0, f )
phase[i,j,k] = 1.0
end
end
end
end
# Quick check with Plots
xmid, ymid, zmid = nx÷2+1, ny÷2+1, nz÷2+1
p1 = heatmap(xc,zc,phase[:,ymid,:]') # xz
p2 = heatmap(zc,yc,phase[xmid,:,:]') # yz
p3 = heatmap(xc,yc,phase[:,:,zmid]') # xy
display(plot(p1,p2,p3))
# Output for Paraview
vtk_grid("RhombicDodecahedron", xc, yc, zc) do vtk
vtk["phase"] = phase
end
end
main( 12 )