# Rhombic Dodecahedron

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, X, X, 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( )
p    .*= r/2.0
# Call fitting function
f  = FitRhombicDocdecahedron( p, f2v )
# Loop through cells and check if in or out
@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 )

``````