# Read final tet mesh created for FEHM modeling
# Define mesh cell materials and outside zones for exodus file
# Conventions were setup by ealier meshing so there is no material 1 or 2
#
# FEHM zone id numbers are changed to fit Exodus conventions
# Exodus is based on element (cell) colors in itetclr
# and has boundaries defined as outside surfaces (tri face of tet)
# Note the node imt zone values are unchanged since Exodus does not use them. 
# Only element itetclr and facesets are used.
#
# For Exodus, we have 3 mesh materials plus top facesets (side sets)
#   1 = background (fehm zone 5) 
#   2 = fault (fehm zone 4)
#   3 = Repo, same as fehm zone 
#
# Faceset integer id is assigned by order they are defined:
#   1 = Top Left and middle (fehm zone 1)
#   2 = Top Right (fehm zone 2)
#   3 = Top Fault (fehm zone overwrites top fault to 4)
# 

read avs tet_final_zones.inp mo
cmo printatt mo -all- minmax

### CELL MATERIAL COLORS ##################
# for Exodus set cells (tets) using node imt colors
# Node colors will not be used. 
# Repository should be about 20 vertical

# set everything to 1, then overwrite with fault and repo
cmo select mo
cmo setatt mo itetclr 1

cmo printatt mo imt minmax

# repository stays 3
# use exclusive 
# inclusive selects too many elements making
# the repository more than 20 m high
pset/p3/attribute imt/1,0,0/ eq 3
eltset/e3/ exclusive pset,get,p3
cmo setatt mo itetclr eltset,get,e3 3

# fault find imt 4, set to 2
# watch that it does not get too wide
pset/p4/attribute imt/1,0,0/ eq 4
eltset/e4/ exclusive pset,get,p4
cmo setatt mo itetclr eltset,get,e4 2

# write element sets for each material
eltset / e_mat3 / itetclr / eq / 3
eltset / e_mat2 / itetclr / eq / 2
eltset / e_mat1 / itetclr / eq / 1
eltset/e_mat3/write / output_mat_3.cellset / ascii
eltset/e_mat2/write / output_mat_2.cellset / ascii
eltset/e_mat1/write / output_mat_1.cellset / ascii
dump gmv tet_colors.gmv mo

###### SORT CELLS for Exodus #################

# use sort to order element blocks as exodus will order
# if this is not done, lagrit faceset numbers will not
# correlate to exodus faceset numbers
# itetclr must be ordered correctly
# additional sort keys can be based on user preference

# sort based on element itetclr number and median location
createpts / median
sort / mo / index / ascending / ikey / itetclr xmed ymed zmed
reorder / mo / ikey
  dump / gmv / out_tmp_sort1.gmv / mo
  cmo / DELATT / mo / ikey

# sort nodes based on mesh convention z, y, then x
sort / mo / index / ascending / ikey / zic yic xic
reorder / mo / ikey
  dump / gmv / out_tmp_sort2.gmv / mo
  cmo / DELATT / mo / xmed
  cmo / DELATT / mo / ymed
  cmo / DELATT / mo / zmed
  cmo / DELATT / mo / ikey
  cmo / DELATT / mo / matclr

###### TOP FACESETS ###########################
# define file names here so when dump/exo is called
# it can use the defines to get the files correct later
define FS_1 exo_1_top_left.faceset
define FS_2 exo_2_top_right.faceset

# extract the outside surface of the full mesh
# to get outside face to element relationships
extract / surfmesh / 1 0 0 / mo_surf / mo / external
cmo / printatt / mo_surf / -all- / minmax

# ## Set surface elements to direction they face based on normals
#  These values will be written to itetclr, copy to id_side attribute
#  1 = bottom   2 = top, 4-6 are sides

cmo / addatt / mo_surf / id_side / vint / scalar / nelements
cmo / select mo_surf
settets / normal
cmo / copyatt / mo_surf mo_surf / id_side itetclr
cmo / printatt / mo_surf / id_side / minmax
dump gmv outside_surface.gmv mo_surf

# make sure to remove all attributes except idelem1 and idface1
cmo / DELATT / mo_surf / itetclr0
cmo / DELATT / mo_surf / idnode0
cmo / DELATT / mo_surf / idelem0
cmo / DELATT / mo_surf / facecol
cmo / DELATT / mo_surf / itetclr1
cmo / DELATT / mo_surf / idface0

############### SET FACES - top left ####################

# TOP LEFT = 1
cmo / copy / mo_tmp / mo_surf
cmo / select / mo_tmp
# changed to lt so face sets do not overlap
# pset / psel / attribute / xic / 1,0,0 / le / 700.0
pset / psel / attribute / xic / 1,0,0 / lt / 700.0
eltset esel inclusive pset,get,psel
eltset /e_top / itetclr eq 2
  eltset / e_keep / inter e_top esel
  eltset / e_delete / not / e_keep
  rmpoint / element / eltset get e_delete
  rmpoint / compress
    cmo / DELATT / mo_tmp / id_side
    dump / avs2 / FS_1 / mo_tmp / 0 0 0 2
    cmo / delete / mo_tmp

# TOP RIGHT = 2
cmo / copy / mo_tmp / mo_surf
cmo / select / mo_tmp
pset / psel / attribute / xic / 1,0,0 / gt / 700.0
eltset esel inclusive pset,get,psel
eltset /e_top / itetclr eq 2
  eltset / e_keep / inter e_top esel
  eltset / e_delete / not / e_keep
  rmpoint / element / eltset get e_delete
  rmpoint / compress
    cmo / DELATT / mo_tmp / id_side
    dump / avs2 / FS_2 / mo_tmp / 0 0 0 2
    cmo / delete / mo_tmp

############################################################
#         FINAL FILES

# write final output mesh and with faceset file
# order of files detirmines integer value of the set
# This exo file can be read into GMV, side sets are under Display Surfaces

dump / exo / tet_exodus.exo / mo / / / facesets &
             FS_1 FS_2 

dump gmv tet_exodus.gmv mo

# DONE
finish