#define PROG_NAME "GeralBarySubdivision"
#define PROG_DESC "???"
#define PROG_VERS "1.0"

#define GeralBarySubdivision_C_COPYRIGHT \
  ""
  
#define PROG_INFO \
  "" \
  " "

#include <Mis.h>
#include <r3.h>
#include <Bary.h>
#include <Octf.h>
#include <argparser.h>
#define _GNU_SOURCE
#include <stdio.h>
#include <Triangulation.h>

CONST
 double order = 2;
   mid   == 1;

TYPE
typedef struct Options_t {
    char *inFile;          /* Initial guess file name (minus ".tp") */
    char *outFile;          /* Output file name prefix */
    bool_t net;
  }

VAR
  ElemTableRec_t top;
  vlt : char *= "VP";
  vlf : char *= "VF";
  bool_t net;

Options_t *GetOptions(int argc, char **argv);

int main(int argc, char **argv)
     
  REF Place_vec_t ps ;
  { Options_t *o = GetOptions(argc, argv);
    char *topo_cmt = jsprintf("Created by %s on %s", PROG_NAME, Today());
    /* Random_t coins = MakeRandomSource(4615); */
    ??? tc = Triangulation.ReadToMa(o->inFile);
    ??? NFE = top->wall.nelE;
    {
      top = tc.top;
      net = o->net;
      Place_vec_t ps = Place_vec_new(NFE);
      fprintf(stderr, "Subdividing from: " & o->inFile & ".tp\n");

      for (i = 0; i < NFE; i++)
      {
        ??? aposp = PposP(top->wedge[i]);
          with (
        double apneg = PnegP(top->wedge[i])
        DO
          if ((aposp != NULL) && (apneg != NULL))
          {
            ps[i] = Bary.MakeWedge(order,order);
            SetTetrahedron2x2(ps[i],i);
          }
          else
          { 
            ps[i] = Bary.MakeWedge(mid,order);
            SetTetrahedron1x2(ps[i],i);
          }
        }  
      }        

      /* gluing tetrahedra */
      for (i = 0; i < NFE; i++)
      {
        with (
        double pa = top->wedge[i]
        DO
          GlueTetra(pa);
        }
      }
        
      ??? newtop = MakeElemTable(PWedge(top->wedge[NFE-1]).ca.el[1];
      with (0),
        nc     == GenCoords(newtop)^,
        com    == "Subdivided from: " & o->inFile & ".tp\n"
                 &"Created by Barycenter Subdivision: "&o->outFile \
                 ".tp on " & Today()
        DO   
          WriteTopology(o->outFile, newtop, com);
          WriteState(o->outFile, newtop, nc, com \
                                   "\nRandom Geometry");
          WriteMaterials(o->outFile, newtop, com);
          
          /* unmark all the wedges */ 
          for (i = 0; i < NFE; i++)
          {
            ??? fet = PWedge(top->wedge[i]);
            {
              fet.marks = FALSE;
            }
          }
        }  
    
    return 0;
  }

bool_t Internal(Place_t @p)
  { if (((PposP(a)!=NULL)) && ((PnegP(a)!=NULL)
    )){ return TRUE;
      }
      else
      { return FALSE;
    }
  } /* END Internal */

bool_t CommonCell(Place_t @p, Place_t b)
  { ??? app = PposP(a);
    Cell_t apn = PnegP(a);
    ??? bpp = PposP(b);
    Cell_t bpn = PnegP(b);
    {
      /*fprintf(stderr, "oi\n");*/
      if ((app == bpn) && (app != NULL)){ return TRUE; }
      if ((apn == bpp) && (apn != NULL)){ return TRUE; }
      return FALSE;
    }
  } /* END CommonCell */

PROCEDURE GlueTetra(Place_t @p)
  /* Glues the topological tetrahedra such as the topology. */ 
  { /*assert(PBits(a) == 0);*/
    Place_t b = NextF(a);
    Place_t c = PrevF(a);
    Place_t d = NextE(a);
    Place_t e = PrevE(a);
    {
      ??? aa = Bary_Corner(a);
        ??? bb = Bary_CCorner(b);
      {
        if ((bb != aa) && (a#b))
        {
          if ((Internal(a)) && (Internal(b)))
          {
            EVAL Triangulation.Glue(bb, aa, order, TRUE);
            /* Update Corners to mark this places as glued */
            Bary_SetCorner(a, bb);
            Bary_SetCorner(Spin(Clock(a)), Bary_CCorner(Spin(Clock(b))));
          }
          else
           if ((CommonCell(a,b))){ 
            EVAL Triangulation.Glue(bb, aa, order, TRUE);
            /* Update Corners to mark this places as glued */
            Bary_SetCorner(a, bb);
            Bary_SetCorner(Spin(Clock(a)), Bary_CCorner(Spin(Clock(b))));
          }
        }
      }
      
      ??? cc = Bary_Corner(c);
        ??? dd = Bary_CCorner(a);
      {
        if ((dd != cc) && (c#a))
        {
          if ((Internal(a)) && (Internal(c)))
          {
            EVAL Triangulation.Glue(dd, cc, order, TRUE);
            /* Update Corners to mark this places as glued */
            Bary_SetCorner(c, dd);
            Bary_SetCorner(Spin(Clock(c)), Bary_CCorner(Spin(Clock(a))));
           }
           else
            if ((CommonCell(a,c))){ 
             EVAL Triangulation.Glue(dd, cc, order, TRUE);
             /* Update Corners to mark this places as glued */
             Bary_SetCorner(c, dd);
             Bary_SetCorner(Spin(Clock(c)), Bary_CCorner(Spin(Clock(a))));
          }
        }
      }
                
      ??? ee = Bary_Corner(Srot(d));
        ??? ff = Bary_CCorner(Srot(a));
      {
        if ((ff != ee) && (a#d))
        {
          /* Not yet glued, glue it: */
          EVAL Triangulation.Glue(ff, ee, order, TRUE); 
          /* Update Corners to mark this places as glued */ 
          Bary_SetCorner(Srot(d), ff);
          Bary_SetCorner(Spin(Tors(d)), Bary_CCorner(Spin(Tors(a))));
        }
      }
          
      ??? gg = Bary_Corner(Srot(a));
        ??? hh = Bary_CCorner(Srot(e));
      {
        if ((gg != hh) && (a#e))
        {
          EVAL Triangulation.Glue(hh, gg, order, TRUE);
          /* Update Corners to mark this places as glued */
          Bary_SetCorner(Srot(a), hh);
          Bary_SetCorner(Spin(Tors(a)), Bary_CCorner(Spin(Tors(e))));
        }
      }
      
    }
  } /* END GlueTetra */

PROCEDURE SetTetrahedron2x2(a:Place_t; uint i)
  { ??? fes = NARROW(PWedge(a)->
    with ( Bary.wedge),
    double fet = PWedge(top->wedge[i]);
      fn  == fet.wall->num, 
      en  == fet.edge->num,
      on  == OrgV(top->wedge[i]),
      dn  == OrgV(Clock(top->wedge[i])),
    double ovr = on.radius;
    double ovc = on.color;
    double ovt = on.transp;
    double ovl = on.label;
    double ove = on->exists;
    double dvr = dn.radius;
    double dvc = dn.color;
    double dvt = dn.transp;
    double dvl = dn.label;
    double dve = dn->exists;
    double era = top->edge[en].radius;
    double eco = top->edge[en].color;
    double eta = top->edge[en].transp;
    double eex = top->edge[en]->exists;
    double ero = top->edge[en]->root;
    double fco = top->wall[fn].color;
    double fta = top->wall[fn].transp;
    double fex = top->wall[fn]->exists;
    double fro = top->wall[fn]->root
  ){
      assert(Octf.OrientationBit(top->wedge[i]) == 0);
      fet.ca = fes.ca;  
      SetDual  (fet.ca.el[3],vlt,vlf,FALSE,TRUE);
      SetPrimal(fet.ca.el[0],ovr,ovc,ovt,ovl,ove,dvr,dvc,
                dvt,dvl,dve,era,eco,eta,eex,ero,fco,fta,
                fex,fro,FALSE,TRUE,net);            
    }
} SetTetrahedron2x2;

PROCEDURE SetTetrahedron1x2(a:Place_t; uint i)
  { ??? fes = NARROW(PWedge(a)->
    with ( Bary.wedge),
    double fet = PWedge(top->wedge[i]);
      fn  == fet.wall->num, 
      en  == fet.edge->num,
      on  == OrgV(top->wedge[i]),
      dn  == OrgV(Clock(top->wedge[i])),
    double ovr = on.radius;
    double ovc = on.color;
    double ovt = on.transp;
    double ovl = on.label;
    double ove = on->exists;
    double dvr = dn.radius;
    double dvc = dn.color;
    double dvt = dn.transp;
    double dvl = dn.label;
    double dve = dn->exists;
    double era = top->edge[en].radius;
    double eco = top->edge[en].color;
    double eta = top->edge[en].transp;
    double eex = top->edge[en]->exists;
    double ero = top->edge[en]->root;
    double fco = top->wall[fn].color;
    double fta = top->wall[fn].transp;
    double fex = top->wall[fn]->exists;
    double fro = top->wall[fn]->root;
      /*         */
    double pa =  top->wedge[i];
    double paf = NextF(pa);
    double fef = PWedge(NextF(pa))
  ){
      fet.ca = fes.ca;
      if ((! fet.marks) && ((PposP(pa) == PnegP(paf))))
      {
        if ((PposP(pa)!=NULL))
        {
          SetDual(fet.ca.el[3],vlt,vlf,TRUE,FALSE);
          SetPrimal(fet.ca.el[0],ovr,ovc,ovt,ovl,ove,dvr,
                    dvc,dvt,dvl,dve,era,eco,eta,eex,ero,
                    fco,fta,fex,fro,TRUE,FALSE,net);        
          fet.marks = TRUE;
        };
        if (! fef.marks)
        {
          SetDual(fet.ca.el[3],vlt,vlf,TRUE,TRUE);
          SetPrimal(fet.ca.el[0],ovr,ovc,ovt,ovl,ove,dvr,
                    dvc,dvt,dvl,dve,era,eco,eta,eex,ero,
                    fco,fta,fex,fro,TRUE,TRUE,net);         
          fef.marks = TRUE;
        }
      }
    }
  } /* END SetTetrahedron1x2 */

PROCEDURE SetPrimal(Place_t @p,     /* @place hanged on topological tetrahedron 2x2 or 1x2 */
    REAL ovr;     /* origen node radius */
    ovc: frgb_t;     /* origen node color */
    ovt: frgb_t;     /* origen node transparency */ 
    char *ovl;     /* origen node label */
    bool_t ove;
    REAL dvr;     /* destine node radius */
    dvc: frgb_t;     /* destine node color */
    dvt: frgb_t;     /* destine node transparency */
    char *dvl;     /* destine node label */
    bool_t dve;
    REAL era;     /* @{edge->?} radius */
    eco: frgb_t;     /* @{edge->?} color  */
    eta: frgb_t;     /* @{edge->?} transparency */
    bool_t eex;  /* @{edge->?} exists */
    INTEGER ero;  /* @{edge->?} root   */
    fco: frgb_t;     /* wall color  */
    fta: frgb_t;     /* wall transparency */
    bool_t fex;  /* wall exists */
    INTEGER fro;  /* wall root   */
    bool_t mid;  /* TRUE iff the topological tetrahedron is 1x2 */
    bool_t side; /* indicate the side of the topological tetrahedron 1x2 */
    bool_t net; /* simulates a grade with thin cylindres and spheres    */
 )

  void SetNode(
     Node v;
     bool_t e; /* exists */
     c: frgb_t;    /* color  */
     t: frgb_t;    /* transp */
     REAL r;    /* radius */
     char *l;    /* label  */
   )
    {
      v->exists = e;
      v.color  = c;
      v.transp = t;
      v.label  = l;
      v.radius = r;
    } SetNode;

  void SetGhostWall(Place_t @p)
    {
      Wall_t t = PWall(a)
        with (
    ){
        t->exists = FALSE;
      }
    } SetGhostWall;

  void NewSetTriangle(
     Place_t b;
     bool_t e;  /* exists */
     c: frgb_t;     /* color  */
     t: frgb_t;     /* transp */
     INTEGER r;  /* root   */
   )
    {
      Wall_t f = PWall(b)
      {
    ){
        f->exists = e;
        f.color  = c;
        f.transp = t;
        f->root   = r;
      }
    } NewSetTriangle;

  void NewSetEdge(
     Place_t b;
     bool_t e;   /* exists */
     c: frgb_t;      /* color  */
     t: frgb_t;      /* transp */
     REAL ra;     /* radius */
     INTEGER ro;  /* root   */
   )
    {
      ??? ee = NARROW(PEdge(b), @{Edge->?});
      {
    ){
        ee->exists = e;
        ee.color  = c;
        ee.transp = t;
        ee.radius = ra;
        ee->root   = ro;
      }
    } NewSetEdge;

  { /* set the origin of the @place ca.el[0] */  
    SetNode(OrgV(a), ove, ovc, ovt, ovr, ovl);

    /* set the origin of the @place Clock(ca.el[0]) */
    SetNode(OrgV(Clock(a)), eex, eco, eta, era, "VE");

    /* set the edge component of the topological tetrahedron */ 
    NewSetEdge(a, eex, eco, eta, era, ero);

    Place_t b = Clock(PrevE(NextF(NextE(a))));
    {
      SetNode(OrgV(Clock(b)), dve, dvc, dvt, dvr, dvl);
      NewSetEdge(b, eex, eco, eta, era, ero);
    }

    if (! mid)
    {
      SetGhostWall(a);
      SetGhostWall(NextF(a));
      NewSetTriangle(PrevF(a), fex, fco, fta, fro);

      Place_t b = Clock(PrevE(NextF(NextE(a))));
      {
        SetGhostWall(b);
        SetGhostWall(NextF(b));
        NewSetTriangle(PrevF(b), fex, fco, fta, fro);
      }
      if (net)
      {
        Edge_t X = PEdge(PrevE(PrevF(a)));
        Edge_t Y = PEdge(NextE(PrevF(a)));
        Place_t Z = Clock(PrevE(NextF(NextE(a))));
        Edge_t W = PEdge(NextE(PrevF(Z)));
        Node_t S = OrgV(PrevE(PrevF(a)));
        with (
        double co = (frgb_t){1.00,1.000,0.500},           /* color, transparency and radius */
        double tp = (frgb_t){0.00,0.000,0.000},           /* of the thin cylinder and sohere*/   
        double ra = 0.0025
      ){
          X->exists = TRUE;
          Y->exists = TRUE;
          W->exists = TRUE;
          S->exists = TRUE;
          X.color  = co;
          X.radius = ra;
          X.transp = tp;
          Y.color  = co;
          Y.radius = ra;
          Y.transp = tp;
          W.color  = co;
          W.radius = ra;
          W.transp = tp;
          S.color  = co;
          S.radius = ra;
          S.transp = tp;          
        }
      }
    }
    else
     if (mid){
      if (side)
      {
        SetGhostWall(a);
        NewSetTriangle(PrevF(a), fex, fco, fta, fro);
        Place_t b = Clock(PrevE(NextF(NextE(a))));
        {
          SetGhostWall(b);
          NewSetTriangle(PrevF(b), fex, fco, fta, fro);
        }
      }
      else
      { 
        SetGhostWall(PrevF(a));
        NewSetTriangle(a, fex, fco, fta, fro);
        Place_t an = Clock(PrevE(PrevF(NextE(PrevF(a)))));
        {
          SetGhostWall(an);
          NewSetTriangle(PrevF(an), fex, fco, fta, fro);
        }
      }
    }
  } /* END SetPrimal */

PROCEDURE SetDual(Place_t @p,
    char *vlt;
    char *vlf;
    bool_t mid;
    bool_t side;
 )

  void SetNode(Node_t v; char *label)
    {
      ??? vv = (Node_t)(v);
      {
    ){
        vv->exists = FALSE;
        vv.label  = label;
      }
    } SetNode;

  void SetGhostWall(Place_t @p)
    {
      Wall_t t = PWall(a)
        with (
    ){
        t->exists = FALSE;
      }
    } SetGhostWall;

  PROCEDURE SetGhostEdge(Place_t @p) ==
    {
      ??? t = NARROW(PEdge(a), @{Edge->?});
        with (
    ){
        t->exists = FALSE;
      }
    } SetGhostEdge;

  void SetRowedge(d: Place_t; r: uint)
    Place_t dn = d;
    VAR {
      if (r == 3)
      {
        for (j = 0; j < r; j++)
        {
          if (j == 0)
          {
            for (i = 0; i < r; i++)
            {
              SetGhostEdge(dn);
              dn = NextE(dn);
            } 
          }
          dn = PrevF(dn);
          SetGhostEdge(NextE(dn));
          SetGhostEdge(PrevE(dn));
        }
      }
      else
       if (r == 2){
        for (j = 0; j < r; j++)
        {
          if (j == 0)
          {
            for (i = 0; i < r; i++)
            {
              SetGhostEdge(dn);
              dn = NextE(dn);
            }
            dn = NextE(dn);
          }
          dn = PrevF(dn);
          SetGhostEdge(NextE(dn));
        }
      }
    } SetRowedge;

  void SetRowTriangle(b: Place_t)
    {
      for (i = 0; i < 3; i++)
      {
        SetGhostWall(b);
        b = PrevF(b);
      }
    } SetRowTriangle;

  { if ((mid) && (side)){
      SetNode(OrgV(a), vlf);
      SetNode(OrgV(Clock(a)), vlt);
      SetRowTriangle(a);
      SetRowedge(a,3);
    }
    else
     if ((mid) && (NOT side)){
      SetNode(OrgV(a), vlt);
      SetNode(OrgV(Clock(a)), vlf);
      SetRowTriangle(a);
      SetRowedge(a,3);
    }
    else
     if (! mid){
      SetNode(OrgV(a), vlt);
      SetNode(OrgV(Clock(a)), vlf);
      SetRowTriangle(a);
      SetRowedge(a,3);
      Place_t an = PrevE(NextF(NextE(a)));
      {
        SetNode(OrgV(an),vlt);
        SetRowTriangle(an);
        SetRowedge(Clock(an),2);
      }
    }
  } /* END SetDual */

Options_t *GetOptions(int argc, char **argv)
  { Options_t *o = (Options_t *)malloc(sizeof(Options_t));
    argparser_t *pp = argparser_new(stderr, argc, argv);
    argparser_set_help(pp, PROG_NAME " version " PROG_VERS ", usage:\n" PROG_HELP);
    argparser_set_info(pp, PROG_INFO);
    argparser_process_help_info_options(pp);
    
    argparser_get_keyword(pp, "-inFile");                               
        o->inFile = argparser_get_next(pp);
        argparser_get_keyword(pp, "-outFile");                               
        o->outFile = argparser_get_next(pp);
        o->net     = argparser_keyword_present(pp, "-net");
        argparser_finish(pp);                                       
      
  -----------------------------------
#define _HELP \
  fprintf(stderr, "Usage: GeralBarySubdivision" \
   " -inFile <name> -outFile <name> \\\n" \
   " [ -net ] \n");
          
END¦                                              
      }
    }
    return o;
  } /* END GetOptions */

{
  DoIt()
} GeralBarySubdivision.


  /* Copyright © 2001    Universidade Estadual de Campinas (UNICAMP)      */

Bary_Corner Bary_CCorner Bary_SetCorner