# Tutorial LaGriT Introduction
# Step 3. Assign Materials to the Tet Mesh 
# Using commands pset, regions, and interpolation

define IN_FILE 02_tet_mesh.inp
define OUT_FILE tet_colors.inp

define IN_FILE 01_hex_mesh.inp 
define OUT_FILE hex_colors.inp 

# Read the IN_FILE to color materials 

read/avs/ IN_FILE / mo_mat
cmo/status/mo_mat/ brief
cmo/select/mo_mat


# -----------------------------------------------------------
# Assign materials by Selected sets

# define some elevations for layers between Z= 0 and 80

define MAT1_Ztop 40.
define MAT2_Ztop 62.


# SELECT NODE SETS based on attribute zic (Z coordinate) 
pset/pmat1/attribute zic/1,0,0/ lt MAT1_Ztop
pset/pmat2/attribute zic/1,0,0/ lt MAT2_Ztop
pset/pmat3/attribute zic/1,0,0/ ge MAT2_Ztop

# SET NODE IMT in order such that last overwrites previous
# layer 2 overwrites layer 3, layer 1 overwrites layer 2
cmo/setatt/mo_mat/imt/ pset,get,pmat3 / 3
cmo/setatt/mo_mat/imt/ pset,get,pmat2 / 2
cmo/setatt/mo_mat/imt/ pset,get,pmat1 / 1

# SELECT ELEMENT SETS based on point sets
# Use inclusive for any node of element is in set

eltset/emat1/inclusive/pset,get,pmat1
eltset/emat2/inclusive/pset,get,pmat2
eltset/emat3/inclusive/pset,get,pmat3

# SET ELEMENT ITETCLR in order
cmo/setatt/mo_mat/itetclr/eltset,get,emat3/ 3
cmo/setatt/mo_mat/itetclr/eltset,get,emat2/ 2
cmo/setatt/mo_mat/itetclr/eltset,get,emat1/ 1

# write a temporary file to view colors so far
dump/ tmp_layers.inp / mo_mat

# -----------------------------------------------------------
# Assign materials by Surfaces and Regions 

# Create surfaces to define a fourth material

cmo / create / motmp
cmo / select / motmp

# Create a quad mesh based on corner coordinates
# The surface must be equal or larger than the mesh domain 

quadxy / 10 5 /-1. -1. 10.   /100.1 -1.  100./ &
               100.1 51. 100./   -1. 51. 10.
createpts/brick/xyz/ 10 5 1 /1,0,0/connect
hextotet/4/ mosurf1 / motmp
cmo / printatt / mosurf1 / -xyz- / minmax
dump/ tmp_surf1.inp / mosurf1

# Create a second surface offset from first
offsetsurf/ mosurf2 / mosurf1 / 15.
dump/ tmp_surf2.inp / mosurf2

# Define mesh objects as surfaces for geometry
# Make mesh current so geometry is assigned to it

cmo / select / mo_mat
surface / s_mosurf1 / reflect / sheet / mosurf1
surface / s_mosurf2 / reflect / sheet / mosurf2

# Define region between the surfaces
region/ r_slant / ge s_mosurf1 and le s_mosurf2

# Select node and element sets in the slanted region
pset/ pslant / region / r_slant
eltset/ eslant / region / r_slant

# SET element and node materials
cmo/ setatt / mo_mat / imt / pset,get,pslant / 4
cmo/ setatt / mo_mat / itetclr / eltset,get,eslant / 4

# set boundary and interface nodes
resetpts/itp

dump / OUT_FILE / mo_mat

cmo / status / mo_mat / brief

# uncomment Early finish to check results
# finish

# -----------------------------------------------------------
# Assign materials by Interpolation
# This will preserve the stair-step interfaces of materials

# remove unneeded mesh objects
cmo/delete/motmp
cmo/delete/mosurf1
cmo/delete/mosurf2
cmo/list

# Read tet mesh from Step 2 
read/avs/ 02_tet_mesh.inp / mo_tet
cmo/select/mo_tet

interpolate/map/mo_tet/ itetclr /1,0,0/ mo_mat itetclr
interpolate/voronoi/mo_tet/ imt /1,0,0/ mo_mat imt

# check interpolated values
cmo/printatt/mo_tet/imt minmax
cmo/printatt/mo_tet/itetclr minmax

# Set boundary nodes and write view file
resetpts/itp
dump/avs/tet_interp_materials.inp mo_tet

 
# -----------------------------------------------------------
# Report Material Quantities

# Report minmax for node and element materials
cmo/select/mo_tet
cmo/printatt/mo_tet/ imt minmax
cmo/printatt/mo_tet/ itetclr minmax

# Report element volumes by material
quality/ volume / material 

# Report node imt quantities
# These are what FEHM will use
dump/zone_imt/ tet / mo_tet

# Report node voronoi volume totals by node material
# First add attribute vorvol, then report sums
cmo addatt/mo_tet/vor_volume/vorvol
pset/p1/attribute/imt/1,0,0/ eq 1
math/sum/mo_tet/vol_tot/pset,get,p1/mo_tet/vorvol
pset/p1/attribute/imt/1,0,0/ eq 2
math/sum/mo_tet/vol_tot/pset,get,p1/mo_tet/vorvol
pset/p1/attribute/imt/1,0,0/ eq 3
math/sum/mo_tet/vol_tot/pset,get,p1/mo_tet/vorvol
pset/p1/attribute/imt/1,0,0/ eq 4
math/sum/mo_tet/vol_tot/pset,get,p1/mo_tet/vorvol


finish