****************************** 
* 
* Carl W. Gable
* gable@lanl.gov
* LaGriT Control File
* http://lagrit.lanl.gov
* http://meshing.lanl.gov
* 
****************************** 
*
* Build tet meshes 
* and insert a fault plane.
*
* Carl Gable
* gable@lanl.gov
*
*
infile build_background_mesh.mlgi
*
*
* This will build a set of triangulated faults and return
* with the faults in a mesh object named:
*
define / FLT_SET / cmo_tri1
*
infile build_fault_set.mlgi
*
*
* Find where the tet mesh is close to the fault plane.
*
compute / distance_field / CMO_NAME / cmo_tri1 / dfield
*
dump / gmv / dfield.gmv / CMO_NAME

cmo / select / CMO_NAME
pset / prefine / attribute / dfield / 1 0 0 / 3000. / lt
eltset / erefine / inclusive / pset get prefine

refine / eltset / eltset get erefine
rmpoint / compress
eltset / erefine / delete
pset   / prefine / delete

dump / gmv / refine1.gmv / CMO_NAME

cmo / DELATT / CMO_NAME / dfield
compute / distance_field / CMO_NAME / cmo_tri1 / dfield
pset / prefine / attribute / dfield / 1 0 0 / 1500. / lt
eltset / erefine / inclusive / pset get prefine
refine / eltset / eltset get erefine
rmpoint / compress

dump / gmv / refine2.gmv / CMO_NAME
*
* Remove nodes that are close to the fault planes.
*
* Some degenerate triangles exist so intersect_elements will have
* a lot of warning messages, but that is ok.
*
intersect_elements / CMO_NAME / cmo_tri1 / in_flt
eltset / e_in_flt / in_flt / gt / 0
rmpoint / element / eltset get e_in_flt
rmpoint / compress

dump / gmv / rmpoint1.gmv / CMO_NAME
*
* Intersect can leave behind some nodes that are very very close
* to the fault plane. Use an alternative method to delete the
* nodes close to the plane.
*
define / RM_DISTANCE / 1500.
*
cmo / copy / cmo_tmp_tri / cmo_tri1
cmo / select / cmo_tmp_tri
refine2d
cmo / DELATT / CMO_NAME / dfield
compute / distance_field / CMO_NAME / cmo_tmp_tri / dfield
cmo / delete / cmo_tmp_tri
cmo / status / brief
cmo / select / CMO_NAME
pset / pdel / attribute / dfield / 1 0 0 / RM_DISTANCE / lt
eltset / edel / inclusive / pset get pdel
rmpoint / element / eltset get edel
rmpoint / compress
eltset / edel / delete

dump / gmv / rmpoint2.gmv / CMO_NAME
*
* Now put the two point distributions together and connect to form tets.
*
cmo / create / cmo_tet2 / / / tet
copypts / cmo_tet2 / cmo_tri1

filter / 1 0 0 / 50.
rmpoint / compress

copypts / cmo_tet2 / CMO_NAME

filter / 1 0 0 / 40.
rmpoint / compress

cmo / printatt / cmo_tet2 / -xyz- / minmax
trans / 1 0 0 / zero
cmo / printatt / cmo_tet2 / -xyz- / minmax

cmo / setatt / CMO_NAME / imt / 1 0 0 / 1
cmo / setatt / CMO_NAME / itp / 1 0 0 / 0
connect / noadd
rmpoint / compress

trans / 1 0 0 / original
cmo / printatt / cmo_tet2 / -xyz- / minmax

dump / gmv / connect.gmv / cmo_tet2

cmo / DELATT / cmo_tet2 / dfield
compute / distance_field / cmo_tet2 / cmo_tri1 / dfield
pset / pfault / attribute / dfield / 1 0 0 / lt / 3.
cmo / addatt / cmo_tet2 / if_flt / vint / scalar / nnodes
cmo / setatt / cmo_tet2 / if_flt / 1 0 0 / 1
cmo / setatt / cmo_tet2 / if_flt / pset get pfault / 2

dump / gmv / tets_w_fault.gmv / cmo_tet2
*
* Identify elements ABOVE and BELOW the fault plane.
* The algorithm tests 
*
surface / s_flt1 / reflect / sheet / cmo_tri1
cmo / select / cmo_tet2
region / r_above / ge s_flt1
region / r_below / lt s_flt1
eltset / e_below / region / r_below

cmo / addatt / cmo_tet2 / id_block / vint / scalar / nelements
cmo / setatt / cmo_tet2 / id_block / 1 0 0 / 1
cmo / setatt / cmo_tet2 / id_block / eltset get e_below / 2

quality / aspect / y
dump / gmv / tets_id_flt.gmv / cmo_tet2
*
* Put the tet mesh and the fault triangles in the same mesh object
* and output a GMV file.
*
cmo/copyatt/cmo_tet2/cmo_tet2/itetclr id_block
addmesh append cmotmp2 cmo_tet2 cmo_tri1
dump/gmv/tets_w_triflts.gmv/cmotmp2
cmo / delete / cmotmp2

finish