# Interpolate ice wedge mesh on to full stacked mesh
# This will overwrite material layers with ice 30

# Input the ice wedge mesh object
read / avs / ice_wedge.inp / mo_ice
cmo/printatt mo_ice -xyz- minmax
cmo/printatt mo_ice itetclr minmax
cmo / addatt / mo_ice / mat_orig / VINT / scalar / nelements
cmo / setatt / mo_ice / itetclr / 30
cmo / setatt / mo_ice / mat_orig / 30

# Input the full stacked mesh
read / avs / full_mesh.inp / mo
cmo/printatt mo -xyz- minmax
cmo / addatt / mo / node_map / VINT / scalar / nnodes
cmo / addatt / mo / node_ice / VINT / scalar / nnodes
cmo / addatt / mo / elem_map / VINT / scalar / nelements
cmo / addatt / mo / mat_orig / VINT / scalar / nelements
cmo / copyatt / mo / mo / mat_orig / itetclr 

# Interpolate ice wedge mesh on to full mesh node attribute node_map
# Save itetclr to mat_orig before overwriting
cmo/printatt mo itetclr minmax
interpolate/map/ mo node_map/1,0,0/mo_ice itetclr
interpolate/map/ mo elem_map/1,0,0/mo_ice itetclr

# CHECK nodes should be values 30 and 31 (outside ice)
cmo / printatt / mo / node_map / minmax

cmo/select/mo

# Find nodes and elems set with ice value 30 
# note inclusion or exclusion options select more or less
# eltset / e_ice / exclusive / pset get p_ice

pset / p_ice / attribute / node_map / 1,0,0 / eq / 30
eltset / e_ice1 / inclusive / pset get p_ice
eltset / e_ice2 / elem_map / eq 30
eltset / e_ice / union e_ice1 e_ice2

# Find elements in layers where ice exists
# Note Material order 1 at top to bottom n
define MAX_ICELAYER 4
eltset e_all itetclr ge 1
eltset / e_maxl / itetclr / gt / MAX_ICELAYER

# CHECK bounds of ice over all materials
eltset e_all_e30 inter e_ice e_all

# Find elements of ice in permitted layers only
eltset / e30 / inter / e_ice e_maxl

# set itetclr to ice material
cmo / setatt / mo / itetclr / eltset,get,e30 / 30

# set node_ice to ice material
pset/pts_ice2 / eltset / e30
cmo setatt mo node_ice 0
cmo setatt mo node_ice pset,get,pts_ice2/ 1

# make sure all nodes on outside faces are captured
pset/pout/attribute itp/pset,get,pts_ice2/ ge 10
pset/pfix/attribute node_map/pset,get,pout/ lt 10
cmo setatt mo node_map pset,get,pfix 3

dump / proto_full_poly.gmv / mo
# remove extra attributes
cmo/DELATT/mo/ ikey_utr
cmo/DELATT/mo/ node_map
cmo/DELATT/mo/ elem_map

# Final Mesh with ice:
dump / proto_full_poly.exo / mo 
dump / proto_full_poly.inp / mo

############# CHECK FINAL MESH #######################                          

# merge for viewing in same space
cmo copy mo_merge mo
cmo setatt mo_merge itetclr 1
cmo setatt mo_ice itetclr 2
addmesh/merge/mo_merge / mo_merge / mo_ice
dump gmv ck_merge_full_ice.gmv mo_merge
cmo delete mo_merge

# ice wedge mesh object
cmo/printatt mo_ice itetclr minmax
cmo/printatt mo_ice mat_orig minmax
cmo/printatt mo_ice -xyz- minmax

# Full stacked mesh
cmo/printatt mo itetclr minmax
cmo/printatt mo mat_orig minmax
cmo/printatt mo -xyz- minmax
cmo/printatt mo -all- minmax

finish