# read 2D cylinder point set and create Delaunay 2D triangulation

define WELL_ID 1
define MAT_ID  2
read avs cylinder_rim.inp mo_rim
read avs cylinder_well.inp mo_well
  cmo setatt mo_rim  imt MAT_ID
  cmo setatt mo_well imt WELL_ID

# merge the points
addmesh/merge/ mo / mo_rim / mo_well
dump gmv tmp_merge_points.gmv mo

# add points to 2D mesh object
cmo create cmotri///triplane     
  copypts / cmotri / mo
  cmo setatt cmotri imt MAT_ID

####################################
# CONNECT into triangulation

cmo select cmotri
cmo / setatt / cmotri / imt 1

connect

cmo / setatt / cmotri / itetclr 1
  rmpoint compress
  resetpts itp

# set distance from center for well point selections
compute/distance_field/cmotri mo_well / wdist

dump / gmv / tmp_tri1.gmv /cmotri
dump / avs / tmp_tri1.inp / cmotri

####################################
# Optimize with refine and smooth

# Refine by halves, max edge starts near 10.
define REFINE1 5.0
refine/rivara///edge/1 0 0/REFINE1///inclusive
dump / gmv / tmp_tri2.gmv / cmotri
dump / avs / tmp_tri2.inp / cmotri

define REFINE1 2.0
refine/rivara///edge/1 0 0/REFINE1///inclusive
dump / gmv / tmp_tri3.gmv / cmotri
dump / avs / tmp_tri3.inp / cmotri

# Smooth into nicer point distribution
# the more times this is repeated, the nicer it is
# Select points OUTSIDE the well for movement

pset/psmooth/attribute wdist/1,0,0/ gt .01 
smooth/position/esug / pset get psmooth
recon / 0
rmpoint compress
  resetpts itp
  dump gmv tmp_smooth1.gmv cmotri

pset psmooth delete
pset/psmooth/attribute wdist/1,0,0/ gt .01 
smooth/position/ esug / pset get psmooth
recon / 0
rmpoint compress
  resetpts itp
  dump gmv tmp_smooth2.gmv cmotri

* smooth entire mesh
* do not smooth too close to wells
pset psmooth delete
pset/psmooth/attribute wdist/1,0,0/ gt .01 
smooth / position / esug / pset,get,psmooth
recon / 0
  rmpoint compress
  resetpts itp
  dump gmv tmp_smooth3.gmv cmotri

# merge thin triangles
massage 5. .8 /0.1/ 1,0,0
  rmpoint compress
  resetpts itp
  dump gmv tmp_mass4.gmv cmotri

* do not smooth too close to wells
pset psmooth delete
pset/psmooth/attribute wdist/1,0,0/ gt .01 
smooth / position / esug / pset,get,psmooth
recon / 0
  rmpoint compress
  resetpts itp
  dump gmv tmp_smooth4.gmv cmotri

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

# Final cell and node colors
# Select center cylinder set of nodes and cells for well
cmo setatt cmotri imt MAT_ID 
cmo setatt cmotri itetclr MAT_ID 

# pset select the center node
pset/px/attribute xic/1,0,0/ eq 0.
pset/pcenter/ attribute yic /pset,get,px/ eq 0.
# eltset select cells connected to center node
eltset/ewell/ inclusive pset,get,pcenter
# pset select all node vertices of well cells
pset/pwell/ eltset / ewell

# overwrite material MAT_ID with WELL_ID for well
cmo setatt cmotri imt pset,get,pwell WELL_ID 
cmo setatt cmotri itetclr eltset,get,ewell WELL_ID 

# reset some lagrit book-keeping arrays
resetpts / parent
rmpoint / compress
resetpts / itp

dump / gmv / tmp_smooth_final.gmv / cmotri
dump / avs / tmp_smooth_final.inp / cmotri

# re-connect with final recon call
# to ensure the triangulation is Delaunay
recon / 1

# Write final triangulation
dump / gmv / cyl_2D_delaunay.gmv / cmotri
dump / avs / cyl_2D_delaunay.inp / cmotri

# Write some mesh information
cmo/status
cmo/printatt//-all-/minmax
quality

# end
finish