mesher.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------

#include <iostream>
#include <fstream>
#include <cstring>
#include <string>
#include <cmath>
#include <stdlib.h>
#include <limits>
#include <algorithm>

#include "vect.h"
#include "mesher_p.h"
#include "mesher_d.h"


using namespace opencarp;

typedef enum {Myocardium=1, Isobath, Anisobath } region_t;

#define BOX_CENTERS_GRID    false
#define NODE_GRID           true

Point p_assign_array( float *p )
{
  Point a;
  a.x = p[0]; a.y=p[1]; a.z=p[2];
  return a;
}

const Point e_circ = {1,0,0};
const Point e_long = {0,1,0};
const Point e_rad  = {0,1,0};



const float CM2UM=1.e4;

class TisAxes {
  public:
    void  set_xi(float xi_) { xi=xi_; }
    void  set_axes(float alpha_, float beta_pr_, float gamma_);
    void  set_bath_axes(bool);
    Point fiber(void)  {return f;}
    Point sheet(void)  {return s;}
  private:
    float xi;
    float alpha;
    float beta_pr;
    float gamma;
    Point f;
    Point s;
    Point sn;
};

void
TisAxes::set_axes(float b, float c, float d)
{
  alpha   = b/180*M_PI;
  beta_pr = c/180*M_PI;
  gamma   = d/180*M_PI;

  Point p(cos(alpha), sin(alpha), sin(gamma));
  f = normalize(p);

  Point sp(0.0, sin(double(beta_pr)), cos(double(beta_pr)));
  float lambda = -dot(f, sp)/dot(f, e_circ);
  s = normalize(sp + scal_X(e_circ,lambda));
}

void
TisAxes::set_bath_axes(bool aniso_bath)
{
        f.y = f.z = 0.0;
  s.x = s.y = s.z = 0.0;

  f.x = aniso_bath?1.0:0.0;
}

class Region {
  public:
    Region(Point p, int b): p0_(p), bath_(b), tag_(-1) {}
    virtual ~Region() {}
    virtual bool inside(Point) = 0;
    bool isbath() { return bath_;  }
    int  tag() { return tag_; }
    void tag(int t) { tag_ = t; }
  protected:
    Point p0_;
    int   bath_;
    int   tag_;
};

class BlockRegion: public Region {
  public:
    BlockRegion(Point p0, Point p, int bth): Region(p0, bth), p1_(p) {}
    virtual bool inside(Point p);
  private:
    Point p1_;
};

class SphericalRegion: public Region {
  public:
    SphericalRegion(Point ctr, float r, int bth):
                    Region(ctr, bth), radius2_(r*r) {}
    virtual bool inside(Point p) { return dist_2(p, p0_)<= radius2_; }
  private:
    float radius2_; //square of radius
};

bool
BlockRegion :: inside( Point p )
{
  return p.x >= p0_.x && p.x <= p1_.x && p.y >= p0_.y && p.y <= p1_.y &&
         p.z >= p0_.z && p.z <= p1_.z;
}


class CylindricalRegion: public Region {
  public:
    CylindricalRegion(Point origin, Point dir, float r, float l, int bth):
                 Region(origin, bth), radius2_(r*r)
                {axis_=normalize(dir);len_=(l==0.)?1.e36:l;}
    virtual bool inside( Point p );
  private:
    Point axis_;
    float radius2_;
    float len_;
};

bool
CylindricalRegion :: inside(Point p)
{
  Point R = p - p0_;
  float d = dot(R, axis_);

  return d>=0 && d<=len_ && mag2(R)-d*d <= radius2_;
}


class Element {
  public:
    Element(int N, const char *t):n_(N),p_(new int[N]),type_(t){}
    ~Element() { delete[] p_; }
    Point centre(Point *ctr);
    int     num(){ return n_; };
    friend std::ostream&  operator<<(std::ostream &, Element& );
    TisAxes  ax;
  protected:
    int      n_;
    int     *p_;
    std::string   type_;
};


Point Element::centre(Point* ctr) {
  Point result = ctr[p_[0]];
  for( int i=1; i<n_; i++ ) result = result + ctr[p_[i]];
  return scal_X( result, 1./(float)n_);
}


class Tetrahedron:public Element {
  public:
  Tetrahedron( int A, int B, int C, int D ): Element(4,"Tt") {
    p_[0]=A; p_[1]=B; p_[2]=C; p_[3]=D; }
  void chkNegVolume(Point *pts);
};


void Tetrahedron::chkNegVolume(Point* pts){
  Point  P01 = pts[p_[1]] - pts[p_[0]];
  Point  P02 = pts[p_[2]] - pts[p_[0]];
  Point  P03 = pts[p_[3]] - pts[p_[0]];
  double vol = det3(P01,P02,P03); // 6x volume
  if( vol < 0. )
    std::swap(p_[0],p_[1]);
}

class Hexahedron:public Element {
  public:
    Hexahedron( int A, int B, int C, int D, int E, int F, int G, int H ): Element(8,"Hx") {
      p_[0]=A; p_[1]=B; p_[2]=C; p_[3]=D; p_[4]=E; p_[5]=F; p_[6]=G; p_[7]=H; }
};


class Quadrilateral:public Element {
  public:
    Quadrilateral( int A, int B, int C, int D ): Element(4,"Qd") {
        p_[0] = A; p_[1] = B; p_[2] = C; p_[3] = D;
    }
};


class Triangle:public Element {
  public:
    Triangle( int A, int B, int C ): Element(3,"Tr") {
        p_[0] = A; p_[1] = B; p_[2] = C;
    }
};


class Line:public Element {
  public:
    Line( int A, int B ): Element(2,"Ln"){p_[0]=A; p_[1]=B;}
};


std::ostream& operator<<( std::ostream& out, Element &e ) {
  out << e.type_ << " ";
  for( int i=0; i<e.num()-1; i++ )
    out << e.p_[i] << " ";
  out << e.p_[e.num()-1];
  return out;
}

std::ostream& operator<<( std::ostream& out, Point p ) {
  out << p.x << " " << p.y << " " << p.z;
  return out;
}


class tmProfile {
  public:
    ~tmProfile() { free(xi); free(ang); }
    int read(char *fname);
    float lookup(float xi);
    int linear(float ang_endo, float ang_epi, float z_endo, float z_epi );
  private:
    int    N;
    float *xi;
    float *ang;
};


int
tmProfile::linear (float ang_endo,float ang_epi, float z_endo, float z_epi)
{
  float dz = 10;         // sample the profile at a 10 um resolution
  if( z_epi==z_endo )    // 2D case
    N = 0;
  else
    N = (z_epi-z_endo)/dz;

  xi  = (float *)malloc(sizeof(float)*(N+1));
  ang = (float *)malloc(sizeof(float)*(N+1));

  if(xi==NULL || ang==NULL)  {
    std::cerr << "Memory allocation failed." << std::endl;
    return -1;
  }
  else if( !N ) {    // 2D
    xi[0]  = 0;
    ang[0] = ang_endo;
  } else  {
    float K = (ang_epi-ang_endo)/(z_epi-z_endo);
    for(int i=0;i<=N;i++)  {
      float z = z_endo+i*dz;
      xi[i] = (z-z_endo)/(z_epi-z_endo)-0.5;
      ang[i] = ang_endo + K*(z-z_endo);
    }
  }
  return 0;
}

float
tmProfile::lookup(float xi_)
{
  if( !N ) return ang[0];

  int i=0;
  while(xi[i]<xi_ && i<N) i++;

  if((i==0)|| (i==N))
    return ang[i];
  else
    return ang[i-1]+(ang[i]-ang[i-1])/(xi[i]-xi[i-1])*(xi_-xi[i-1]);
}

int
tmProfile::read(char *fname)
{
  FILE *profile = fopen(fname,"rt");
  if(profile==NULL)  {
    fprintf( stderr, "Can't open transmural profile data file %s.\n", fname);
    return -1;
  }

  int err = fscanf(profile,"%d",&N);
  xi  = (float *)malloc(sizeof(float)*N);
  ang = (float *)malloc(sizeof(float)*N);
  for(int i=0; i<N; i++)  {
    xi[i] = (float)i/(float)(N-1)-0.5;
    err = fscanf(profile,"%f",ang+i);
  }
  fclose(profile);

  return err;
}

class BBoxDef  {
  public:
    void update( Point p );
//  private:
    float  xd   = 0;
    float  x_mn = std::numeric_limits<float>::max();
    float  x_mx = std::numeric_limits<float>::min();
    float  yd   = 0;
    float  y_mn = std::numeric_limits<float>::max();
    float  y_mx = std::numeric_limits<float>::min();
    float  zd   = 0;
    float  z_mn = std::numeric_limits<float>::max();
    float  z_mx = std::numeric_limits<float>::min();
};

void
BBoxDef::update( Point p)
{
  if(p.x<x_mn) x_mn = p.x;
  if(p.y<y_mn) y_mn = p.y;
  if(p.z<z_mn) z_mn = p.z;
  if(p.x>x_mx) x_mx = p.x;
  if(p.y>y_mx) y_mx = p.y;
  if(p.z>z_mx) z_mx = p.z;
}


class BoundingBox {
  public:
    virtual ~BoundingBox() { delete[] bx; delete[] res; }
    void  init( int *_bx, float *_res, Point p0);
    void  dims(void);
    float z2xi(float z, bool nodeGrid);
    float mn_z(bool nodeGrid){ return nodeGrid?nodes.z_mn:bctrs.z_mn; };
    float mx_z(bool nodeGrid){ return nodeGrid?nodes.z_mx:bctrs.z_mx; };
  public:
    int     *bx;
    int      bx_inds[3][2];
  protected:
    float   *res;
    BBoxDef  bctrs;
    BBoxDef  nodes;
};

float
BoundingBox::z2xi (float z, bool nodeGrid)
{
  float zd   = nodeGrid?nodes.zd:bctrs.zd;
  float z_mn = nodeGrid?nodes.z_mn:bctrs.z_mn;

  if(zd==0.)
    return 0.;
  else
    return (z-z_mn)/zd-0.5;
}

void
BoundingBox::init(int *_bx, float *_res, Point p0)
{
  bx  = new int [3];
  res = new float [3];
  for(int i=0;i<3;i++)  {
    bx [i] = _bx [i];
    res[i] = _res[i];
    bx_inds[i][0] = 0;
    bx_inds[i][1] = bx[i]-1;
  }

  // min corner of nodal bbx
  nodes.update(p0);

  // min corner of center bbx
  Point pc0   = p0;
  pc0.x = p0.x + res[0]/2;
  pc0.y = p0.y + res[1]/2;
  pc0.z = p0.z + res[2]/2;
  bctrs.update(pc0);

  // max corner of nodal bbx
  Point p1 = p0;
  p1.x = p0.x+bx[0]*res[0];
  p1.y = p0.y+bx[1]*res[1];
  p1.z = p0.z+bx[2]*res[2];
  nodes.update(p1);

  // max corner of center bbx
  Point pc1 = p1;
  pc1.x = p1.x - res[0]/2;
  pc1.y = p1.y - res[1]/2;
  pc1.z = p1.z - res[2]/2;
  bctrs.update(pc1);

  dims();
}

void
BoundingBox::dims(void)
{
  // figure out dims in terms of nodes
  nodes.xd = nodes.x_mx - nodes.x_mn;
  nodes.yd = nodes.y_mx - nodes.y_mn;
  nodes.zd = nodes.z_mx - nodes.z_mn;

  // and in terms of box centers
  bctrs.xd = bctrs.x_mx - bctrs.x_mn;
  bctrs.yd = bctrs.y_mx - bctrs.y_mn;
  bctrs.zd = bctrs.z_mx - bctrs.z_mn;
}

class fibDef {
  public:
    void setFiberDefs(char *f_prof, float fEndo, float fEpi, float imbr,
                      char *s_prof, float sEndo, float sEpi);
    bool withSheets(void);
    char *f_name(void)  { return f_Prof; };
    char *s_name(void)  { return s_Prof; };
    float imbrication() { return f_imbr; };
    float rotEndo()     { return f_Endo; };
    float rotEpi()      { return f_Epi;  };
    float sheetEndo()   { return s_Endo; };
    float sheetEpi()    { return s_Epi;  };

  private:
    char *f_Prof;
    float f_Endo;
    float f_Epi;
    float f_imbr;
    char *s_Prof;
    float s_Endo;
    float s_Epi;
};

void fibDef::setFiberDefs(char *f_prof, float fEndo, float fEpi, float imbr,
                      char *s_prof, float sEndo, float sEpi)
{
  f_Prof = strdup(f_prof);
  f_Endo = fEndo;
  f_Epi  = fEpi;
  f_imbr = imbr;
  s_Prof = strdup(s_prof);
  s_Endo = sEndo;
  s_Epi  = sEpi;
}

bool
fibDef::withSheets(void)
{
  if((f_imbr!=0.0) ||  (s_Endo!=0.0) || (s_Epi!=0.0)  || (strcmp(s_Prof,"")))
    return true;
  else
    return false;

}


class Grid {
  public:
    Grid( char *, Region **, int );
    virtual ~Grid() {}
    void unPrMFiberDefs(void);
    virtual void  build_mesh(float*,float*,float*,bool *,float*,float,bool,int)=0;
    virtual void  output_boundary( char * )=0;
    void add_element( Element &, region_t );
    void  set_indx_bounds(bool *sym);
    bool  chk_bath( int i, int j, int k );
    bool  os_good();
  protected:
    Point       *pt;
    int          npt;
    BoundingBox  b_bbx;
    BoundingBox  t_bbx;
    std::ofstream pt_os, elem_os, lon_os, elemc_os, vec_os;
    Region     **region;
    int          num_axes;
    fibDef       f_def;
    tmProfile    f_xi;
    tmProfile    s_xi;
    int          dim;
};

class Grid2D : public Grid {
  public:
    virtual void build_mesh(float*, float*, float*, bool *, float*, float, bool, int);
    virtual void  output_boundary( char *fn ){};
    Grid2D( char *m, Region **r ):Grid(m,r,2){}
  private:
    void add_tri( int, int, int );
};

class Grid3D : public Grid {
  public:
    virtual void  build_mesh(float*, float*, float*, bool *, float*, float, bool, int);
    virtual void  output_boundary( char *fn ){}
    Grid3D( char *m, Region **r ):Grid(m,r,3){}
  private:
    void add_tet( int, int, int, int );
};

class Grid1D : public Grid {
  public:
    virtual void  build_mesh(float*, float*, float*, bool *, float*, float, bool, int);
    virtual void  output_boundary( char *fn ){}
    Grid1D( char *m, Region **r ):Grid(m,r,1){}
  private:
    void add_line( int, int, int, int );
};


bool Grid::os_good()
{
  return pt_os.good() && lon_os.good() && elem_os.good();
}


Grid::Grid( char *msh, Region **r, int d ): region(r),dim(d)
{
  std::string fname(msh);
  fname  = msh;
  fname += ".lon";
  lon_os.open(fname.c_str());

  fname  = msh;
  fname += ".pts";
  pt_os.open(fname.c_str());

  fname  = msh;
  fname += ".elem";
  elem_os.open(fname.c_str());

  fname  = msh;
  fname += ".vpts";
  elemc_os.open(fname.c_str());

  fname  = msh;
  fname += ".vec";
  vec_os.open(fname.c_str());
}


void
Grid::set_indx_bounds(bool *sym)
{
  for(int i=0;i<3;i++)  {
    if(b_bbx.bx[i]>t_bbx.bx[i])  {
      t_bbx.bx_inds[i][0] = 0;
      if(sym[i])
        t_bbx.bx_inds[i][0] = (b_bbx.bx[i]-t_bbx.bx[i])/2;
      t_bbx.bx_inds[i][1] = t_bbx.bx_inds[i][0]+t_bbx.bx[i]-1;
    }
  }
}


bool
Grid::chk_bath(int i, int j, int k)
{
  bool bth[] = { false, false, false };
  int  idx[] = { i, j, k };

  for(int cnt=0; cnt<dim; cnt++ )
    if( idx[cnt]<t_bbx.bx_inds[cnt][0] || idx[cnt]>t_bbx.bx_inds[cnt][1] )
      return true;

  return false;
}


void
Grid::unPrMFiberDefs(void)
{
  f_def.setFiberDefs(param_globals::fibers.tm_fiber_profile,
                     param_globals::fibers.rotEndo,
                     param_globals::fibers.rotEpi,
                     param_globals::fibers.imbrication,
                     param_globals::fibers.tm_sheet_profile,
                     param_globals::fibers.sheetEndo,
                     param_globals::fibers.sheetEpi);
}

void
Grid::add_element( Element &elem, region_t regtype )
{
  elem_os << elem;
  Point c = elem.centre(pt);
  elemc_os << c << std::endl;
  int r=-1;
  int regid = regtype==Myocardium?1:0;
  while( region[++r] != NULL ) {
    if( region[r]->inside( c ) ) {
      if( region[r]->isbath() ) {            // myo or bath becomes bath
        regid   = region[r]->tag()<0 ? 0 : region[r]->tag();
        regtype = region[r]->isbath()==1?Isobath:Anisobath;
      } else if( regtype == Myocardium ) {   // only myo becomes other myo
        regid   = region[r]->tag()<0 ? r+2 : region[r]->tag();
      }
      if(param_globals::first_reg)
        break;
    }
  }

  if(regtype==Myocardium)  {
    float xi = t_bbx.z2xi(c.z,BOX_CENTERS_GRID);
    elem.ax.set_xi(xi);
    elem.ax.set_axes( f_xi.lookup(xi),s_xi.lookup(xi),f_def.imbrication());
  }
  else  {
    elem.ax.set_bath_axes(regtype==Anisobath);
    if(regtype==Anisobath)
      regid *= -1;
  }

  elem_os << " " << regid << std::endl;
  if(num_axes>1)
    lon_os << elem.ax.fiber() << " " << elem.ax.sheet() << std::endl;
  else
    lon_os << elem.ax.fiber() << std::endl;

  // write fiber to vector file
  vec_os << elem.ax.fiber() << std::endl;
}


void
Grid1D::build_mesh(float* x0, float* x, float* tissue, bool *sym, float *res,
                   float pert, bool aniso_bath, int periodic_bc)
{
  bool eletype;
  int  p1, p2, p3, p4;
  int  bnx[3];            // number of cubes in each direction
  int  tnx[3];            // number of tissue cubes in each direction

  // determine number of boxes
  for( int i=0; i<3; i++ )  {
    bnx[i] = (int)(x[i]/res[i]);
    tnx[i] = (int)(tissue[i]/res[i]);
  }

  // determine origins for tissue and bath bounding box
  Point pt0, pb0;
  pt0.x = -tnx[0]*res[0]/2 + x0[0];
  pt0.y = 0;
  pt0.z = 0;

  pb0.x = sym[0]?pt0.x-(bnx[0]-tnx[0])*res[0]/2:pt0.x+x0[0];
  pb0.y = 0.0;
  pb0.z = 0.0;

  // determine bounding box of bath and tissue grids
  b_bbx.init(bnx,res,pb0);
  t_bbx.init(tnx,res,pt0);

  // set tissue index boundaries relative to bath grid
  set_indx_bounds(sym);

  f_xi.linear( f_def.rotEndo(), f_def.rotEpi(), 0, 0 );
  s_xi.linear( f_def.sheetEndo(), f_def.sheetEpi(), 0, 0 );

  // output the points file
  npt = (bnx[0]+1);
  std::cout << "Number of points: " << npt << std::endl;
  pt = new Point[npt];
  pt_os << npt << std::endl;
  npt = 0;
  for( int i=0; i<=bnx[0]; i++ ) {
    pt[npt].assign<double>(x0[0]+i*res[0], 0, 0);

    if( i!=bnx[0] ) {
      pt[npt].x += ((double)random()/(double)RAND_MAX)*pert*res[0];
    }
    pt_os << pt[npt++] << std::endl;
  }
  pt_os.close();

  int num_in_layer = (bnx[0]+1);
  elem_os << bnx[0] << std::endl;
  std::cout << "Number of linear elements: " << bnx[0] << std::endl;

  for( int i=0; i<bnx[0]; i++ ) {
    p1 = i;
    p2 = p1+1;

    region_t regtype = Myocardium;
    if( chk_bath(i,0,0) )
      regtype = aniso_bath?Anisobath:Isobath;

    Line l(p1, p2 );
    add_element( l, regtype );
  }
  elem_os.close();
  lon_os.close();
  delete[] pt;
}


void
Grid2D::build_mesh(float* x0, float* x, float* tissue, bool *sym, float *res,
                   float pert, bool aniso_bath, int periodic_bc)
{
  bool eletype;
  int  p1, p2, p3, p4;
  int  bnx[3];            // number of cubes in each direction
  int  tnx[3];            // number of tissue cubes in each direction

  // determine number of boxes
  for( int i=0; i<3; i++ )  {
    bnx[i] = (int)(x[i]/res[i]);
    tnx[i] = (int)(tissue[i]/res[i]);
  }

  // determine origins for tissue and bath bounding box
  Point pt0, pb0;
  pt0.x = -tnx[0]*res[0]/2 + x0[0];
  pt0.y = -tnx[1]*res[1]/2 + x0[1];
  pt0.z = 0.0;

  pb0.x = sym[0]?pt0.x-(bnx[0]-tnx[0])*res[0]/2:pt0.x+ x0[0];
  pb0.y = sym[0]?pt0.y-(bnx[0]-tnx[0])*res[0]/2:pt0.y+ x0[1];
  pb0.z = 0.0;

  // determine bounding box of bath and tissue grids
  b_bbx.init(bnx,res,pb0);
  t_bbx.init(tnx,res,pt0);

  // set tissue index boundaries relative to bath grid
  set_indx_bounds(sym);

  // output the points file
  npt = (bnx[0]+1)*(bnx[1]+1);
  std::cout << "Number of points: " << npt << std::endl;
  pt = new Point[npt];
  pt_os << npt << std::endl;
  npt = 0;
  for( int j=0; j<=bnx[1]; j++ )
    for( int i=0; i<=bnx[0]; i++ ) {
      pt[npt].assign<double>(x0[0]+i*res[0], x0[1]+j*res[1], 0);

      if( j && j!=bnx[1] && i && i!=bnx[0] ) {
        pt[npt].x += ((double)random()/(double)RAND_MAX)*pert*res[0];
        pt[npt].y += ((double)random()/(double)RAND_MAX)*pert*res[1];
      }
      pt_os << pt[npt++] << std::endl;
    }
  pt_os.close();

  // determine periodic b.c. connections
  int nper_cnnx = 0;
  if( (periodic_bc&1) == 1 )
    nper_cnnx += bnx[1]+1;
  if( (periodic_bc&2) == 2 )
    nper_cnnx += bnx[0]+1;

  int num_in_layer = (bnx[0]+1)*(bnx[1]+1);

  if(param_globals::tri2D) {
    elem_os  << 2*bnx[0]*bnx[1]+nper_cnnx << std::endl;
    elemc_os << 2*bnx[0]*bnx[1]+nper_cnnx << std::endl;
    std::cout << "Number of triangles: " << 2*bnx[0]*bnx[1] << std::endl;
  } else {
    elem_os  << bnx[0]*bnx[1]+nper_cnnx << std::endl;
    elemc_os << bnx[0]*bnx[1]+nper_cnnx << std::endl;
    std::cout << "Number of quadrilaterals: " << bnx[0]*bnx[1] << std::endl;
  }
  if( nper_cnnx )
    std::cout << "Number of periodic connections: " << nper_cnnx << std::endl;

  f_xi.linear( f_def.rotEndo(), f_def.rotEpi(), 0, 0 );
  s_xi.linear( f_def.sheetEndo(), f_def.sheetEpi(), 0, 0 );

  num_axes = f_def.withSheets()?2:1;
  lon_os << num_axes << std::endl;
  if(num_axes>1)
    std::cout << "Using orthotropic fiber setup." << std::endl;
  else
    std::cout << "Using transversely isotropic fiber setup." << std::endl;
    
  for( int j=0; j<bnx[1]; j++ ) {
    eletype = j%2;
    for( int i=0; i<bnx[0]; i++ ) {
      p1 = j*(bnx[0]+1) + i;
      p2 = p1+1;
      p3 = (j+1)*(bnx[0]+1) + i;
      p4 = p3+1;

      region_t regtype = Myocardium;
      if( chk_bath(i,j,0) )
        regtype = aniso_bath?Anisobath:Isobath;

      if( param_globals::tri2D ) {
        if( eletype ) {
          Triangle t1(p1, p2, p3), t2(p2, p4, p3);
          add_element( t1, regtype );
          add_element( t2, regtype );
        } else {
          Triangle t1(p1, p2, p4), t2(p1, p4, p3);
          add_element( t1, regtype );
          add_element( t2, regtype );
        }
        eletype = !eletype;
      } else {
        Quadrilateral q(p1,p2,p4,p3);
        add_element( q, regtype );
      }
    }
  }

  if( (periodic_bc&1) == 1 ) {
    for( int i=0; i<=bnx[1]; i++ ) {
      Line l( i*(bnx[0]+1), (i+1)*(bnx[0]+1)-1 );
      elem_os << l << " " << param_globals::periodic_tag<<std::endl;
      lon_os << "1 0 0" << std::endl;
    }
  }
  if( (periodic_bc&2) == 2 ) {
    for( int i=0; i<=bnx[0]; i++ ) {
      Line l( i, (bnx[0]+1)*(bnx[1])+i );
      elem_os << l << " " << param_globals::periodic_tag+1 << std::endl;
      lon_os << "0 1 0" << std::endl;
    }
  }

  elem_os.close();
  lon_os.close();
  delete[] pt;
}


void
Grid3D::build_mesh( float* x0, float* x, float* tissue, bool *sym, float *res,
                    float pert, bool aniso_bath, int periodic_bc)
{
  bool eletype;
  int  p1, p2, p3, p4, p5, p6, p7,p8;
  int  bnx[3];            // number of cubes in each direction
  int  tnx[3];            // number of tissue cubes in each direction

  // determine number of boxes
  for( int i=0; i<3; i++ )  {
    bnx[i] = (int)(x[i]/res[i]);
    tnx[i] = (int)(tissue[i]/res[i]);
  }

  // determine origins for tissue and bath bounding box
  Point pt0, pb0;
  pt0 = {x0[0],x0[1],x0[2]};

  pb0.x = sym[0]?pt0.x-(bnx[0]-tnx[0])*res[0]/2:pt0.x + x0[0];
  pb0.y = sym[1]?pt0.y-(bnx[1]-tnx[1])*res[1]/2:pt0.y + x0[1];
  pb0.z = sym[2]?pt0.z-(bnx[2]-tnx[2])*res[2]/2:pt0.z + x0[2];

  // determine bounding box of bath and tissue grids
  b_bbx.init(bnx,res,pb0);
  t_bbx.init(tnx,res,pt0);

  // set tissue index boundaries relative to bath grid
  set_indx_bounds(sym);


  // define fiber and sheet arrangements
  if(strcmp(f_def.f_name(),""))  {
    std::cout << "Reading transmural fiber rotation profile from " << f_def.f_name() << std::endl;
    f_xi.read(f_def.f_name());
  }
  else
    f_xi.linear(f_def.rotEndo(),f_def.rotEpi(),t_bbx.mn_z(BOX_CENTERS_GRID),
                t_bbx.mx_z(BOX_CENTERS_GRID));

  if(strcmp(f_def.s_name(),""))  {
    std::cout << "Reading transmural sheet profile from " << f_def.s_name() << std::endl;
    s_xi.read(f_def.s_name());
  }
  else
    s_xi.linear(f_def.sheetEndo(),f_def.sheetEpi(),t_bbx.mn_z(BOX_CENTERS_GRID),
                t_bbx.mx_z(BOX_CENTERS_GRID));

  // output the points file
  npt = (bnx[0]+1)*(bnx[1]+1)*(bnx[2]+1);
  std::cout << "Number of points: " << npt << std::endl;
  pt = new Point[npt];
  pt_os << npt << std::endl;
  npt = 0;
  for( int k=0; k<=bnx[2]; k++ )
    for( int j=0; j<=bnx[1]; j++ )
      for( int i=0; i<=bnx[0]; i++ ) {
        pt[npt] = {x0[0]+i*res[0], x0[1]+j*res[1], x0[2]+k*res[2]};
        if( k && k!=bnx[2] && j && j!=bnx[1] && i && i!=bnx[0] ) {
          pt[npt].x += ((double)random()/(double)RAND_MAX)*pert*res[0];
          pt[npt].y += ((double)random()/(double)RAND_MAX)*pert*res[1];
          pt[npt].z += ((double)random()/(double)RAND_MAX)*pert*res[2];
        }
        pt_os << pt[npt++] << std::endl;
      }
  pt_os.close();

  int num_in_layer = (bnx[0]+1)*(bnx[1]+1);
  if(!param_globals::Elem3D)  {
    elem_os << 5*bnx[0]*bnx[1]*bnx[2] << std::endl;
    elemc_os << 5*bnx[0]*bnx[1]*bnx[2] << std::endl;
    std::cout << "Number of Tetrahedra: " << 5*bnx[0]*bnx[1]*bnx[2] << std::endl;
  }
  else  {
    elem_os << bnx[0]*bnx[1]*bnx[2] << std::endl;
    elemc_os << bnx[0]*bnx[1]*bnx[2] << std::endl;
    std::cout << "Number of Hexahedra: " << bnx[0]*bnx[1]*bnx[2] << std::endl;
  }

  num_axes = f_def.withSheets()?2:1;
  lon_os << num_axes << std::endl;
  if(num_axes>1)
    std::cout << "Using orthotropic fiber setup." << std::endl;
  else
    std::cout << "Using transversely isotropic fiber setup." << std::endl;


  for( int k=0; k<bnx[2]; k++ ) {
    eletype = k%2;
    for( int j=0; j<bnx[1]; j++ ) {
      // if bnx[0] is even we need to toggle
      // otherwise we break the chessboard pattern
      if(j && !(bnx[0]%2) )
        eletype = !eletype;

      for( int i=0; i<bnx[0]; i++ ) {
        p1 = k*num_in_layer + j*(bnx[0]+1) + i;
        p2 = p1+1;
        p3 = k*num_in_layer + (j+1)*(bnx[0]+1) + i;
        p4 = p3+1;
        p5 = (k+1)*num_in_layer + j*(bnx[0]+1) + i;
        p6 = p5+1;
        p7 = (k+1)*num_in_layer + (j+1)*(bnx[0]+1) + i;
        p8 = p7+1;

        region_t regtype = Myocardium;
        if(chk_bath(i,j,k))
          regtype = aniso_bath?Anisobath:Isobath;

        if(!param_globals::Elem3D) {
          // generate tet mesh
          if( eletype ) {
            Tetrahedron t1(p1, p2, p4, p6),
                        t2(p1, p3, p4, p7),
                        t3(p1, p4, p6, p7),
                        t4(p1, p5, p6, p7),
                        t5(p4, p6, p7, p8);
            t1.chkNegVolume( pt ); add_element( t1, regtype );
            t2.chkNegVolume( pt ); add_element( t2, regtype );
            t3.chkNegVolume( pt ); add_element( t3, regtype );
            t4.chkNegVolume( pt ); add_element( t4, regtype );
            t5.chkNegVolume( pt ); add_element( t5, regtype );
          } else {
            Tetrahedron t1(p1, p2, p3, p5),
                        t2(p2, p3, p4, p8),
                        t3(p2, p5, p6, p8),
                        t4(p2, p3, p5, p8),
                        t5(p3, p5, p7, p8);
            t1.chkNegVolume( pt ); add_element( t1, regtype );
            t2.chkNegVolume( pt ); add_element( t2, regtype );
            t3.chkNegVolume( pt ); add_element( t3, regtype );
            t4.chkNegVolume( pt ); add_element( t4, regtype );
            t5.chkNegVolume( pt ); add_element( t5, regtype );
          }
          eletype = !eletype;
        }
        else  {
          // Hex

          // this is the original hex ordering. not compatible with the SF ansatzfunc scheme..
          // Hexahedron h1(p4, p3, p1, p2, p8, p6, p5, p7);

          Hexahedron h1(p5, p7, p8, p6, p1, p2, p4, p3);
          add_element( h1, regtype );
        }
      }
    }
  }
  elem_os.close();
  lon_os.close();
  elemc_os.close();
  vec_os.close();
  delete[] pt;
}


int main( int argc, char *argv[] )
{
  int   status;
  float tsize[3], x0[3], t0[3];
  bool  symbath[3];

  do{
    status = param(PARAMETERS, &argc, argv );
    if( status==PrMERROR||status==PrMFATAL )
      fprintf( stderr, "\n*** Error reading parameters\n\n");
    else if( status==PrMQUIT )
      fprintf( stderr, "\n*** Quitting by user's request\n\n");
  }while (status==PrMERROR);
  if(status&&PrMERROR) exit( status );


  for( int i=0; i<3; i++ ) {
    param_globals::size[i]   *= CM2UM;
    param_globals::bath[i]   *= CM2UM;
    param_globals::center[i] *= CM2UM;
    // in 2D case we don't allow a bath attached in the direction
    // perpendicular to the 2D surface
    if(param_globals::size[i]==0.) param_globals::bath[i] = 0.0;
    if(param_globals::bath[i]>=0)  {
      tsize[i]   = param_globals::size[i]+param_globals::bath[i];
      symbath[i] = false;
      x0[i]      = -0.5*param_globals::size[i]+param_globals::center[i];
    }
    else  {
      tsize[i]   = param_globals::size[i]-2*param_globals::bath[i];
      symbath[i] = true;
      x0[i]      = -0.5*tsize[i]+param_globals::center[i];
    }
  }

  Region **region = (Region**)calloc(param_globals::numRegions+1,sizeof(Region*));
  Point ctr;
  ctr.x = param_globals::center[0];
  ctr.y = param_globals::center[1];
  ctr.z = param_globals::center[2];

  RegionDef* regdef = param_globals::regdef;

  for (int r = 0; r < param_globals::numRegions; r++ ) {
    Point p0 = scal_X( p_assign_array(regdef[r].p0), CM2UM );
    if( !param_globals::size[2]) p0.z = 0;
    Point p1;
    switch(regdef[r].type) {
      case 0:
        if( !param_globals::size[2] ) p1.z = 0;
        p1 = scal_X( p_assign_array( regdef[r].p1 ), CM2UM );
        p1 = p1 + ctr;
        if( p0.x>p1.x )  std::swap( p0.x, p1.x );
        if( p0.y>p1.y )  std::swap( p0.y, p1.y );
        if( p0.z>p1.z )  std::swap( p0.z, p1.z );
        region[r] = new BlockRegion( p0, p1, regdef[r].bath );
        break;
      case 1:
        regdef[r].rad *= CM2UM;
        region[r] = new SphericalRegion( p0,regdef[r].rad, regdef[r].bath );
        break;
      case 2:
        Point axis;
        if( !param_globals::size[2] )
          axis = {0,0,1};
        else
          axis = p_assign_array( regdef[r].p1 );
        regdef[r].rad    *= CM2UM;
        regdef[r].cyllen *= CM2UM;
        region[r] = new CylindricalRegion( p0, axis, regdef[r].rad,
                                            regdef[r].cyllen, regdef[r].bath );
         break;
    }
    region[r]->tag(regdef[r].tag);
  }
  region[param_globals::numRegions] = NULL;


  Grid *grid;
  if( !param_globals::size[1] && !param_globals::size[2] )
    grid = new Grid1D( param_globals::mesh, region );
  else if( !param_globals::size[2] )
    grid = new Grid2D( param_globals::mesh, region );
  else  {
    grid = new Grid3D( param_globals::mesh, region );
  }

  if( grid->os_good() ) {
    grid->unPrMFiberDefs();
    grid->build_mesh(x0, tsize, param_globals::size, symbath, param_globals::resolution,
                     param_globals::perturb, param_globals::anisoBath, param_globals::periodic);
    delete grid;
  }
  else {
    std::cerr << "ERROR:" << std::endl;
    std::cerr << "Could not open mesh files " << param_globals::mesh << ".* for writing!" << std::endl;
    std::cerr << "Aborting!" << std::endl;
    delete grid;
    return 1;
  }

  return 0;
}