# Regular Hex Mesh with Materials from surface regions
# This mesh is truncated by the top surface 
# This will result in stair-stepped geometries
# But will connect into a Delaunay tet mesh


#-------------------------------------------------
# Create the hex mesh
# This can be done with a single call to createpts
# Or with multiple calls with various resolutions

define/XMIN/ 0.
define/YMIN/ 0.
define/ZMIN/ 0.
define/XMAX/ 20.
define/YMAX/ 10.
define NX 21
define NY 11

# Vertical spacing
define/ZMAX/ 18. 
define/Z2  / 10. 
define/Z1  / 5.
define/ZMIN/ 0.

cmo/create/cmohex///hex
cmo/select/cmohex

# Vertical spacing from middle to top 
createpts/brick/xyz/NX NY 28 / XMIN YMIN Z2 / XMAX YMAX ZMAX 
dump/tmp_hex01.inp/ cmohex

# Vertical spacing from middle 1 to middle 2 
createpts/brick/xyz/NX NY 11 / XMIN YMIN Z1 / XMAX YMAX Z2 
dump/tmp_hex0102.inp/ cmohex

# Vertical spacing bottom to middle
createpts/brick/xyz/NX NY 6 / XMIN YMIN ZMIN / XMAX YMAX Z1 
dump/tmp_hex010203.inp/ cmohex

cmo/printatt/cmohex/ -xyz- minmax
quality


#-------------------------------------------------
# Create wave surfaces for regions 
# translate to make multiple surfaces
# make sure surfaces extend beyond boundaries of hex mesh 
# z(i)=sin(x(i))
# x coordinates 1-10 gives 2 complete 2pi cycles (0 to 6.28)
# x coordinates 1-20 gives 4 cycles 

# Use x coordinates to compute z based on sin()
# y coordinate detirmines length of surface
define / XMAX / 20. 
define / NX / 40
define / YMAX / 10. 
define / NY / 10
define SCALE_VAL 2.

cmo / create / mosurf / / / quad
cmo select mosurf

quadxy /NX NY/0. 0. 0./ XMAX 0. 0./XMAX YMAX 0./0. YMAX 0./ 1 1 1
createpts/brick/xyz/NX,NY,1/1 0 0 / connect
math/sin/mosurf/zic/1,0,0/mosurf/xic
scale / 1 0 0 / relative / xyz / SCALE_VAL SCALE_VAL SCALE_VAL 
cmo / setatt / mosurf / imt / 1 0 0 / 1
cmo / setatt / mosurf / itetclr / 1 0 0 / 1
resetpts / itp
cmo / printatt / mosurf / -xyz- / minmax

trans / 1 0 0 / 0. 0. 0. / -1. -1. 8.0
dump / surf_sin_01.inp / mosurf
cmo/printatt/mosurf/ zic minmax

cmo/copy/mosurf2/mosurf
cmo/select/mosurf2
trans / 1 0 0 / 0. 0. 0. / -1. 0. 6.8
dump / surf_sin_02.inp / mosurf2
cmo/printatt/mosurf2/ zic minmax

#-------------------------------------------------
# Define Geometry
# Color Hex Mesh with regions defined by surfaces
# Save this hex mesh for coloring tets into hex shapes 

cmo select cmohex
surface / s1 / intrface / sheet  / mosurf
surface / s2 / intrface / sheet  / mosurf2

# set interfaces equal to one surface
region / r1 / lt s1 and le s2
region / r2 / gt s1 and lt s2
region / r3 / ge s2
 
# interfaces are not equal to surface
mregion / mr1 / lt s1 and lt s2
mregion / mr2 / gt s1 and lt s2
mregion / mr3 / gt s2

setpts
settets
resetpts/itp

# Write the colored hex mesh
dump/ hex.inp/ cmohex

#-------------------------------------------------
# CREATE TET MESH
# Connect point distribution into tet mesh

cmo/create/cmotet
copypts/cmotet/cmohex

# remove duplicate points, reset attributes
filter/1,0,0
rmpoint/compress
cmo / setatt / cmotet / imt / 1 0 0 / 1
cmo / setatt / cmotet / itp / 1 0 0 / 0 

connect noadd
resetpts / itp
quality


#-------------------------------------------------
# COLOR MATERIALS and TRUNCATE TOP
# Interpolate material colors from hex mesh
# Remove material above top surface

interpolate/map/cmotet itetclr /1,0,0/ cmohex itetclr
interpolate/voronoi/cmotet imt /1,0,0/ cmohex imt
cmo/printatt/cmotet/itetclr minmax
dump/tet.inp/cmotet

# remove cells above material 2
# check for material 3 nodes left behind
eltset/edel/ itetclr / gt 2
rmpoint element eltset,get,edel
rmpoint/compress
resetpts/itp

cmo/printatt/cmotet/itetclr minmax
cmo/printatt/cmotet/imt minmax

# FIX imt of nodes left behind after removal of cells
define CMO cmotet
define MAX_MAT 2
infile reset_truncated_nodes_imt.mlgi

#-------------------------------------------------
# SET OUTSIDE BOUNDARY NODES
# tag outside nodes and save attributes
# we need to use the max node located in each xyz direction
# normal zone_outside will tag stair-steps incorrectly

dump/ zone_outside_minmax/ tet_max /cmotet/ keepatt
dump/ tet_trunc.inp/ cmotet

#-------------------------------------------------
# check quality
# check geometric coefficients
cmo/setatt/cmotet/imt 1	
dump/ stor/ tet_trunc / cmotet
quality

# EXIT
finish