# Creates a 3D 3D 3D ice wedge

# read the outer polygon points (no elements)
# connect points into flat shape of triangles
cmo / create / mo / / / triplane
read / avs / poly40_sim1.inp / mo
  triangulate / counterclockwise
  cmo / setatt / mo / imt / 1,0,0 / 1
  cmo / setatt / mo / itetclr / 1,0,0 / 1
  resetpts / itp
  dump / tmp_poly40_tri1.gmv / mo

# check triangulation
quality

# read the inner polygon points (no elements)
# connect points into flat shape of triangles
cmo / create / mo_inner / / / triplane
read / avs / poly40_sim2.inp / mo_inner
  triangulate / counterclockwise
  cmo / setatt / mo_inner / imt / 1,0,0 / 2
  cmo / setatt / mo_inner / itetclr / 1,0,0 / 2
  resetpts / itp 
  dump / gmv / tmp_poly40_tri2.gmv / mo_inner

############### CHECK tri surfs #############

# the number of triangles is related to number of poly points
cmo select mo
quality
cmo select mo_inner
quality

# Both shapes should be within same area, 2nd poly inside
cmo printatt mo -xyz- minmax
cmo printatt mo_inner -xyz- minmax

# Write merged poly shapes
addmesh / merge / mo_merge / mo / mo_inner
  dump / gmv / tmp_poly40.gmv / mo_merge

# There should be 2 materials
cmo / printatt / mo_merge / imt / minmax
cmo / printatt / mo_merge / itetclr / minmax

#############################################

# Break big triangles into smaller for better distribution
cmo select mo_merge
massage / 5 / 1.e-5 / 1.e-5 / strictmergelength
smooth / position / esug / 1,0,0;recon 0;
smooth / position / esug / 1,0,0;recon 0;
smooth / position / esug / 1,0,0;recon 0;
smooth / position / esug / 1,0,0;recon 0;
smooth / position / esug / 1,0,0;recon 0;
smooth / position / esug / 1,0,0;recon 0;

massage / 1 / 1.e-5 / 1.e-5 / strictmergelength
smooth / position / esug / 1,0,0;recon 0;
smooth / position / esug / 1,0,0;recon 0;
smooth / position / esug / 1,0,0;recon 0;
smooth / position / esug / 1,0,0;recon 0;
smooth / position / esug / 1,0,0;recon 0;
smooth / position / esug / 1,0,0;recon 0;
filter / 1,0,0
rmpoint / compress

# Flatten input elevations and interpolate to poly shape
read / gmv / lidar_site3.gmv / mo_lidar
#read / gmv / poly40_topo.gmv / mo_lidar
cmo / copyatt / mo_lidar / mo_lidar / z_val / zic
cmo / setatt / mo_lidar / zic / 1,0,0 / 0.
cmo / addatt / mo_merge / z_val / vdouble / scalar / nnodes
intrp/continuous/mo_merge z_val/1,0,0/mo_lidar z_val/linear
  cmo / copyatt / mo_merge / mo_merge / zic / z_val
  dump / gmv / tmp_poly_massage.gmv / mo_merge
  cmo / printatt / mo_merge / zic / minmax

# REFINE - Depends on Parameter for horizontal spacing (edge lengths)
cmo / copy / mo_poly_tri / mo_merge
cmo / select / mo_poly_tri
cmo / setatt / mo_poly_tri / imt / 1
cmo / setatt / mo_poly_tri / itetclr / 1
cmo / setatt / mo_poly_tri / itp / 0
resetpts / itp
# each call to refine2d will split edges into smaller triangles
refine2d
refine2d
#refine2d
  dump / tmp_poly40_refine.inp / mo_poly_tri
  resetpts itp
  cmo delete mo_poly_tri

# Extract the edge lines around the polygon
read / avs / tmp_poly40_refine.inp / mo_poly_tri
extract / surfmesh / 1,0,0 / mo_line_work / mo_poly_tri
cmo / status / brief
dump / gmv / tmp_line.gmv / mo_line_work

# Create distance attibute from outer polygon edge
compute / distance_field / mo_merge / mo_line_work / dfield
cmo printatt mo_merge dfield minmax
dump / tmp_poly40_refine.gmv / mo_merge

# Sizes for ice wedge used in math routines to define shape
define / HALF_WIDTH / 1.4
define / DEPTH / 2.8
define / SLOPE / 2.0

define / L1 / 0.02
define / L2 / 0.24
define / L3 / 0.1
define / L4 / 0.05

# Create top and bottom layer of ice wedge mesh object
# top of ice wedge is translated to requested layer
# to put top under a specified layer
# These distances are same as defined for full mesh
define / L1 / 0.02
# this is all material 2 
define / L2 / 0.24
# this is all material 3 
define / L3 / 0.24
# this is all material 4 
define / L4 / 0.5
# this is all material 5
define / L5 / 5

# subtract each previous layer until ice top is positioned correctly
# CHECK Starting elevation is same as top elevation
cmo printatt mo_merge zic minmax

math / sub / mo_merge / zic / 1,0,0 / mo_merge / zic / L1
#dump / avs / tmp_poly40_1.inp / mo_merge

math / sub / mo_merge / zic / 1,0,0 / mo_merge / zic / L2
#dump / avs / tmp_poly40_2.inp / mo_merge

math / sub / mo_merge / zic / 1,0,0 / mo_merge / zic / L3
#dump / avs / tmp_poly40_3.inp / mo_merge

#math / sub / mo_merge / zic / 1,0,0 / mo_merge / zic / L4
#math / sub / mo_merge / zic / 1,0,0 / mo_merge / zic / L5
#dump / avs / tmp_poly40_4.inp / mo_merge

# top file for stack routine
dump / avs / tmp_top.inp / mo_merge

# CHECK top elevation should be located at intended layer
cmo printatt mo_merge zic minmax

# math routines to pinchout ice wedge at the bottom
cmo / addatt / mo_merge / z_ele / vdouble / scalar / nnodes
cmo / copyatt / mo_merge / mo_merge / z_save / zic
cmo / printatt / mo_merge / zic / minmax
cmo / addatt / mo_merge / depth / vdouble / scalar / nnodes
cmo / setatt / mo_merge / depth / 1,0,0 / DEPTH
cmo / printatt / mo_merge / imt / minmax
#eltset / e2 / itetclr / eq / 2
#pset / p2 / eltset / e2
#cmo / setatt / mo_merge / imt / pset get p2 / 2
pset / pdown / attribute / dfield / 1,0,0 / lt / HALF_WIDTH
#pset / pdown / attribute / imt / 1,0,0 / eq / 1
#math / divide / mo_merge / slope / pset get pdown / mo_merge / depth / mo_merge / dfield
math / multiply / mo_merge / z_ele / pset get pdown / mo_merge / dfield / SLOPE 
math / sub / mo_merge / z_ele / pset get pdown / mo_merge / z_ele / DEPTH
math / add / mo_merge / z_ele / 1,0,0 / &
mo_merge / z_ele / mo_merge / z_save
#math / ceiling / mo_merge / z_save / 1,0,0 / mo_merge / z_save / 0.

cmo / printatt / mo_merge / z_ele / minmax
cmo / copyatt / mo_merge / mo_merge / zic / z_ele

# bottom layer for stack routine
dump / avs / tmp_bot.inp / mo_merge

# CHECK bottom elevation should be DEPTH below top layer 
cmo printatt mo_merge zic minmax

# stack top and bottom into a prism mesh ice wedge object
cmo / create / mo_stack
stack / layers / avs / tmp_bot.inp 1 tmp_top.inp 1 / flip
stack / fill / mo_prism / mo_stack
filter / 1,0,0
rmpoint / compress
rmpoint / zero_volume / 1.e-12

# Final output file to be used later 
dump / gmv / tmp_ice_wedge.gmv / mo_prism

############# CHECK FINAL WORK #####################
# write surface triangulation sizes
cmo select mo_poly_tri
quality
quality/edge_min
quality/edge_max

# coordinates for polygon surface (top elevations)
cmo printatt mo_poly_tri -xyz- minmax

# coordinates for stacked ice mesh object
cmo printatt mo_prism -xyz- minmax

# there should be one material for ice mesh object
cmo printatt mo_prism imt minmax
cmo printatt mo_prism itetclr minmax

##################################################

# DONE 
# file tmp_poly40_refine.inp surface for full mesh layers
# file tmp_ice_wedge.gmv 3D ice wedge mesh for interpolation

finish