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 <cctype>
#include <cstdio>
#include <cstring>
#include <string>
#include <cmath>
#include <stdlib.h>
#include <limits>
#include <algorithm>
#include <sys/wait.h>
#include <unistd.h>
#include <vector>
 
#include "vect.h"
#include <mesher_p.h>
#include <mesher_d.h>
#include "mesher_schema.hpp"
#include "runtime.hpp"
#include "snapshot_file_io.hpp"
 
 
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;
}
 
namespace {
 
enum class ParserCompareMode {
  Off,
  Warn,
  Strict,
};
 
enum class ParserFallbackMode {
  Off,
  Legacy,
};
 
struct RuntimeCompatOptions {
  ParserFallbackMode fallback_mode = ParserFallbackMode::Off;
};
 
struct LegacyCompareInput {
  bool                     available = false;
  std::vector<std::string> runtime_args;
  std::string              unavailable_reason;
};
 
std::string trim_copy(const std::string& value)
{
  std::string::size_type first = 0;
  while (first < value.size() && std::isspace(static_cast<unsigned char>(value[first]))) {
    ++first;
  }
 
  std::string::size_type last = value.size();
  while (last > first && std::isspace(static_cast<unsigned char>(value[last - 1]))) {
    --last;
  }
 
  return value.substr(first, last - first);
}
 
std::string to_lower_ascii(std::string value)
{
  for (std::string::size_type i = 0; i < value.size(); ++i) {
    value[i] = static_cast<char>(std::tolower(static_cast<unsigned char>(value[i])));
  }
  return value;
}
 
bool parse_option_argument(int*               index,
                           int                argc,
                           char**             argv,
                           const std::string& attached_value,
                           const char*        option_name,
                           std::string*       value,
                           std::string*       error)
{
  if (!attached_value.empty()) {
    *value = attached_value;
    return true;
  }
  if (*index + 1 >= argc) {
    *error = "Missing argument after " + std::string(option_name);
    return false;
  }
  *value = argv[++(*index)];
  return true;
}
 
bool filename_has_suffix(const std::string& filename, const char* suffix)
{
  const std::string normalized_filename = to_lower_ascii(trim_copy(filename));
  const std::string normalized_suffix = to_lower_ascii(std::string(suffix == NULL ? "" : suffix));
  return normalized_filename.size() >= normalized_suffix.size() &&
         normalized_filename.compare(normalized_filename.size() - normalized_suffix.size(),
                                     normalized_suffix.size(),
                                     normalized_suffix) == 0;
}
 
bool match_long_option(const std::string& token, const char* option, std::string* attached_value)
{
  *attached_value = std::string();
  if (token == option) {
    return true;
  }
 
  const std::string prefix = std::string(option) + "=";
  if (token.size() > prefix.size() && token.compare(0, prefix.size(), prefix) == 0) {
    *attached_value = token.substr(prefix.size());
    return true;
  }
 
  return false;
}
 
bool is_help_topic_candidate(const char* token)
{
  return token != NULL && token[0] != '\0' && token[0] != '-' && token[0] != '+';
}
 
bool is_long_option_argument_error(const std::string& token,
                                   const char*        option,
                                   const std::string& attached_value,
                                   std::string*       error)
{
  if (attached_value.empty()) {
    return false;
  }
 
  *error = "Unexpected argument for " + std::string(option) + " in '" + token + "'";
  return true;
}
 
bool normalize_runtime_args(int                        argc,
                            char**                     argv,
                            std::vector<std::string>* normalized,
                            RuntimeCompatOptions*      compat,
                            std::string*               error)
{
  normalized->clear();
  compat->fallback_mode = ParserFallbackMode::Off;
 
  if (argc <= 0 || argv == NULL || argv[0] == NULL) {
    *error = "Missing program name";
    return false;
  }
 
  normalized->push_back(argv[0]);
 
  for (int i = 1; i < argc; ++i) {
    const std::string token = argv[i];
    if (token == "+") {
      break;
    }
 
    std::string attached_value;
 
    if (token == "+Help" || match_long_option(token, "--help", &attached_value)) {
      std::string topic = "PrM";
      if (!attached_value.empty()) {
        topic = attached_value;
      } else if (i + 1 < argc && is_help_topic_candidate(argv[i + 1])) {
        topic = argv[++i];
      }
      normalized->push_back("+Help");
      normalized->push_back(topic);
      break;
    }
 
    if (token == "+Doc" || match_long_option(token, "--doc", &attached_value)) {
      std::string topic = "ALL";
      if (!attached_value.empty()) {
        topic = attached_value;
      } else if (i + 1 < argc && is_help_topic_candidate(argv[i + 1])) {
        topic = argv[++i];
      }
      normalized->push_back("+Doc");
      normalized->push_back(topic);
      break;
    }
 
    if (token == "+Default" || match_long_option(token, "--default", &attached_value)) {
      if (token != "+Default" && is_long_option_argument_error(token, "--default", attached_value, error)) {
        return false;
      }
      normalized->push_back("+Default");
      break;
    }
 
    if (token == "+Run" || match_long_option(token, "--run", &attached_value)) {
      if (token != "+Run" && is_long_option_argument_error(token, "--run", attached_value, error)) {
        return false;
      }
      normalized->push_back("+Run");
      continue;
    }
 
    if (token == "+I" || match_long_option(token, "--interactive", &attached_value)) {
      if (!attached_value.empty()) {
        *error = "Unexpected argument for --interactive in '" + token + "'";
      } else {
        *error = "Unsupported option " + token + " (interactive mode is not available)";
      }
      return false;
    }
 
    if (match_long_option(token, "--param-fallback", &attached_value)) {
      std::string mode = attached_value;
      if (mode.empty()) {
        if (i + 1 >= argc) {
          *error = "Missing argument after --param-fallback";
          return false;
        }
        mode = argv[++i];
      }
 
      if (to_lower_ascii(trim_copy(mode)) != "legacy") {
        *error = "Unsupported value '" + mode + "' for --param-fallback (expected legacy)";
        return false;
      }
 
      compat->fallback_mode = ParserFallbackMode::Legacy;
      continue;
    }
 
    if (token == "+F" || match_long_option(token, "--file", &attached_value)) {
      std::string filename = attached_value;
      if (filename.empty()) {
        if (i + 1 >= argc) {
          *error = "Missing filename after " + token;
          return false;
        }
        filename = argv[++i];
      }
      normalized->push_back("+F");
      normalized->push_back(filename);
      continue;
    }
 
    if (token == "+Save" || match_long_option(token, "--save", &attached_value)) {
      std::string filename = attached_value;
      if (filename.empty()) {
        if (i + 1 >= argc) {
          *error = "Missing argument after " + token;
          return false;
        }
        filename = argv[++i];
      }
      normalized->push_back("+Save");
      normalized->push_back(filename);
      continue;
    }
 
    normalized->push_back(token);
  }
 
  return true;
}
 
bool build_legacy_compare_input(int argc, char** argv, LegacyCompareInput* input, std::string* error)
{
  input->available = false;
  input->runtime_args.clear();
  input->unavailable_reason.clear();
 
  if (argc <= 0 || argv == NULL || argv[0] == NULL) {
    *error = "Missing program name";
    return false;
  }
 
  input->runtime_args.push_back(argv[0]);
  bool saw_passthrough_token_before_file = false;
 
  for (int i = 1; i < argc; ++i) {
    const std::string token = argv[i];
    if (token == "+") {
      break;
    }
 
    std::string attached_value;
 
    if (match_long_option(token, "--param-fallback", &attached_value)) {
      std::string ignored;
      if (!parse_option_argument(&i, argc, argv, attached_value, "--param-fallback", &ignored, error)) {
        return false;
      }
      continue;
    }
 
    if (token == "+Save" || match_long_option(token, "--save", &attached_value)) {
      std::string ignored;
      if (!parse_option_argument(&i, argc, argv, attached_value, token.c_str(), &ignored, error)) {
        return false;
      }
      continue;
    }
 
    if (token == "+Help" || match_long_option(token, "--help", &attached_value)) {
      std::string topic = "PrM";
      if (!attached_value.empty()) {
        topic = attached_value;
      } else if (i + 1 < argc && is_help_topic_candidate(argv[i + 1])) {
        topic = argv[++i];
      }
      input->runtime_args.push_back("+Help");
      input->runtime_args.push_back(topic);
      input->available = true;
      return true;
    }
 
    if (token == "+Doc" || match_long_option(token, "--doc", &attached_value)) {
      std::string topic = "ALL";
      if (!attached_value.empty()) {
        topic = attached_value;
      } else if (i + 1 < argc && is_help_topic_candidate(argv[i + 1])) {
        topic = argv[++i];
      }
      input->runtime_args.push_back("+Doc");
      input->runtime_args.push_back(topic);
      input->available = true;
      return true;
    }
 
    if (token == "+Default" || match_long_option(token, "--default", &attached_value)) {
      if (token != "+Default" && is_long_option_argument_error(token, "--default", attached_value, error)) {
        return false;
      }
      input->runtime_args.push_back("+Default");
      continue;
    }
 
    if (token == "+Run" || match_long_option(token, "--run", &attached_value)) {
      if (token != "+Run" && is_long_option_argument_error(token, "--run", attached_value, error)) {
        return false;
      }
      input->runtime_args.push_back("+Run");
      continue;
    }
 
    if (token == "+F" || match_long_option(token, "--file", &attached_value)) {
      std::string filename;
      if (!parse_option_argument(&i, argc, argv, attached_value, token.c_str(), &filename, error)) {
        return false;
      }
 
      if (saw_passthrough_token_before_file) {
        input->unavailable_reason =
            "legacy compare is unavailable when direct parameter arguments precede a .par input file";
        input->runtime_args.clear();
        input->runtime_args.push_back(argv[0]);
        return true;
      }
 
      if (!filename_has_suffix(filename, ".par")) {
        input->unavailable_reason =
            "legacy compare is unavailable for non-.par input file '" + filename + "'";
        input->runtime_args.clear();
        input->runtime_args.push_back(argv[0]);
        return true;
      }
 
      input->runtime_args.push_back("+F");
      input->runtime_args.push_back(filename);
      continue;
    }
 
    input->runtime_args.push_back(token);
    saw_passthrough_token_before_file = true;
  }
 
  input->available = input->runtime_args.size() > 1;
  if (!input->available && input->unavailable_reason.empty()) {
    input->unavailable_reason = "unable to reconstruct a legacy-compatible parameter input";
  }
  return true;
}
 
void populate_arg_pointers(const std::vector<std::string>& values, std::vector<char*>* argv)
{
  argv->assign(values.size(), NULL);
  for (std::size_t i = 0; i < values.size(); ++i) {
    (*argv)[i] = const_cast<char*>(values[i].c_str());
  }
}
 
void print_lines(FILE* stream, const char* label, const std::vector<std::string>& lines)
{
  for (std::size_t i = 0; i < lines.size(); ++i) {
    std::fprintf(stream, "%s%s\n", label, lines[i].c_str());
  }
}
 
ParserCompareMode parser_compare_mode()
{
  const char* raw = std::getenv("OPENCARP_PARAM_COMPARE");
  if (raw == NULL) {
    return ParserCompareMode::Strict;
  }
 
  const std::string normalized = to_lower_ascii(trim_copy(raw));
  if (normalized.empty() || normalized == "1" || normalized == "on" || normalized == "true" ||
      normalized == "yes" || normalized == "strict" || normalized == "fail" || normalized == "error") {
    return ParserCompareMode::Strict;
  }
  if (normalized == "warn") {
    return ParserCompareMode::Warn;
  }
  if (normalized == "0" || normalized == "off" || normalized == "false" || normalized == "no") {
    return ParserCompareMode::Off;
  }
 
  return ParserCompareMode::Strict;
}
 
ParserFallbackMode parser_fallback_mode()
{
  const char* raw = std::getenv("OPENCARP_PARAM_FALLBACK");
  if (raw == NULL) {
    return ParserFallbackMode::Off;
  }
 
  const std::string normalized = to_lower_ascii(trim_copy(raw));
  if (normalized.empty() || normalized == "0" || normalized == "off" || normalized == "false" ||
      normalized == "no") {
    return ParserFallbackMode::Off;
  }
  if (normalized == "legacy") {
    return ParserFallbackMode::Legacy;
  }
 
  return ParserFallbackMode::Off;
}
 
paramschema::ExecutionResult execute_parser_runtime_args(const std::vector<std::string>& runtime_args,
                                                         const bool                      allow_save)
{
  std::vector<char*> argv;
  populate_arg_pointers(runtime_args, &argv);
 
  paramschema::ExecutionOptions options;
  options.allow_save = allow_save;
  return paramschema::execute_legacy_cli(paramschema::mesher_schema(),
                                         static_cast<int>(argv.size()),
                                         argv.data(),
                                         options);
}
 
bool apply_parser_runtime_args(const std::vector<std::string>& runtime_args,
                               const bool                      allow_save,
                               paramschema::ExecutionResult*   executed_out)
{
  const paramschema::ExecutionResult executed = execute_parser_runtime_args(runtime_args, allow_save);
  if (executed_out != NULL) {
    *executed_out = executed;
  }
  print_lines(stderr, "parameter parser warning: ", executed.warnings);
  if (!executed.rendered_output.empty()) {
    std::fputs(executed.rendered_output.c_str(), stdout);
  }
  if (!executed.errors.empty() || executed.status == paramschema::ExecutionStatus::Fatal) {
    print_lines(stderr, "parameter parser error: ", executed.errors);
    return false;
  }
 
  if (executed.status == paramschema::ExecutionStatus::Help) {
    exit(EXIT_SUCCESS);
  }
 
  return true;
}
 
std::string parent_directory(const std::string& path)
{
  const std::string::size_type slash = path.rfind('/');
  if (slash == std::string::npos) {
    return ".";
  }
  if (slash == 0) {
    return "/";
  }
  return path.substr(0, slash);
}
 
std::string join_path(const std::string& left, const std::string& right)
{
  if (left.empty() || left == ".") {
    return right;
  }
  if (!left.empty() && left[left.size() - 1] == '/') {
    return left + right;
  }
  return left + "/" + right;
}
 
bool find_executable_on_path(const std::string& name, std::string* resolved_path)
{
  if (name.empty() || name.find('/') != std::string::npos) {
    return false;
  }
 
  const char* path_env = std::getenv("PATH");
  if (path_env == NULL || path_env[0] == '\0') {
    return false;
  }
 
  const std::string path_list = path_env;
  std::string::size_type start = 0;
  while (start <= path_list.size()) {
    std::string::size_type end = path_list.find(':', start);
    if (end == std::string::npos) {
      end = path_list.size();
    }
 
    const std::string directory = path_list.substr(start, end - start);
    const std::string candidate = join_path(directory.empty() ? "." : directory, name);
    if (access(candidate.c_str(), X_OK) == 0) {
      *resolved_path = candidate;
      return true;
    }
 
    if (end == path_list.size()) {
      break;
    }
    start = end + 1;
  }
 
  return false;
}
 
bool resolve_legacy_snapshot_helper(const std::string& program_path, std::string* helper_path)
{
  std::vector<std::string> resolved_program_paths;
  if (!program_path.empty()) {
    resolved_program_paths.push_back(program_path);
  }
 
  if (program_path.find('/') == std::string::npos) {
    std::string resolved_program_path;
    if (find_executable_on_path(program_path, &resolved_program_path)) {
      resolved_program_paths.push_back(resolved_program_path);
    }
  }
 
  for (std::size_t i = 0; i < resolved_program_paths.size(); ++i) {
    const std::string program_dir = parent_directory(resolved_program_paths[i]);
    const std::string parent_dir = parent_directory(program_dir);
 
    const std::vector<std::string> candidates = {
        join_path(program_dir, "mesher-param-parser-legacy-snapshot"),
        join_path(join_path(parent_dir, "tools/mesher"), "mesher-param-parser-legacy-snapshot"),
    };
 
    for (std::size_t j = 0; j < candidates.size(); ++j) {
      if (access(candidates[j].c_str(), X_OK) == 0) {
        *helper_path = candidates[j];
        return true;
      }
    }
  }
 
  return find_executable_on_path("mesher-param-parser-legacy-snapshot", helper_path);
}
 
bool create_temp_output_path(const char* suffix, std::string* path)
{
  char temp_path[128];
  if (suffix != NULL && suffix[0] != '\0') {
    std::snprintf(temp_path, sizeof temp_path, "/tmp/mesher-parser-compare-XXXXXX%s", suffix);
  } else {
    std::snprintf(temp_path, sizeof temp_path, "/tmp/mesher-parser-compare-XXXXXX");
  }
 
  const int fd = suffix != NULL && suffix[0] != '\0' ? mkstemps(temp_path, static_cast<int>(std::strlen(suffix))) :
                 mkstemp(temp_path);
  if (fd < 0) {
    return false;
  }
  close(fd);
  *path = temp_path;
  return true;
}
 
bool capture_current_parser_snapshot(paramschema::SnapshotResult* snapshot)
{
  *snapshot = paramschema::snapshot_schema_state(paramschema::mesher_schema());
  print_lines(stderr, "parameter compare warning: ", snapshot->warnings);
  if (!snapshot->errors.empty()) {
    print_lines(stderr, "parameter compare error: ", snapshot->errors);
    return false;
  }
  return true;
}
 
void maybe_force_test_mismatch(paramschema::SnapshotResult* snapshot)
{
  const char* raw = std::getenv("OPENCARP_PARAM_TEST_FORCE_MISMATCH");
  if (raw == NULL) {
    return;
  }
 
  const std::string normalized = to_lower_ascii(trim_copy(raw));
  if (normalized.empty() || normalized == "0" || normalized == "off" || normalized == "false" ||
      normalized == "no") {
    return;
  }
 
  if (!snapshot->entries.empty()) {
    snapshot->entries[0].value += "__forced_parser_compare_mismatch__";
    return;
  }
 
  paramschema::SnapshotEntry entry;
  entry.path = "mesh";
  entry.value = "__forced_parser_compare_mismatch__";
  snapshot->entries.push_back(entry);
}
 
bool run_legacy_snapshot_helper(const std::string&              helper_path,
                                const std::vector<std::string>& runtime_args,
                                const std::string&              output_path)
{
  std::vector<std::string> child_values;
  child_values.reserve(runtime_args.size() + 4);
  child_values.push_back(helper_path);
  child_values.push_back("--snapshot-out");
  child_values.push_back(output_path);
  child_values.push_back("--");
  child_values.insert(child_values.end(), runtime_args.begin(), runtime_args.end());
 
  std::vector<char*> child_argv;
  child_argv.reserve(child_values.size() + 1);
  for (std::size_t i = 0; i < child_values.size(); ++i) {
    child_argv.push_back(const_cast<char*>(child_values[i].c_str()));
  }
  child_argv.push_back(NULL);
 
  const pid_t pid = fork();
  if (pid < 0) {
    std::perror("fork");
    return false;
  }
 
  if (pid == 0) {
    execv(helper_path.c_str(), child_argv.data());
    std::perror("execv");
    _exit(127);
  }
 
  int status = 0;
  if (waitpid(pid, &status, 0) < 0) {
    std::perror("waitpid");
    return false;
  }
 
  if (!WIFEXITED(status) || WEXITSTATUS(status) != 0) {
    if (WIFSIGNALED(status)) {
      std::fprintf(stderr, "parameter compare error: legacy snapshot helper terminated with signal %d\n",
                   WTERMSIG(status));
    } else {
      std::fprintf(stderr, "parameter compare error: legacy snapshot helper exited with status %d\n",
                   WEXITSTATUS(status));
    }
    return false;
  }
 
  return true;
}
 
bool legacy_fallback_requested(const RuntimeCompatOptions& compat)
{
  return compat.fallback_mode == ParserFallbackMode::Legacy ||
         parser_fallback_mode() == ParserFallbackMode::Legacy;
}
 
void print_legacy_fallback_workaround()
{
  std::fprintf(stderr,
               "parameter compare error: rerun with OPENCARP_PARAM_FALLBACK=legacy to continue with the legacy parameter state\n");
  std::fprintf(stderr,
               "parameter compare error: or add --param-fallback=legacy to the mesher command line\n");
}
 
bool restore_legacy_snapshot_state(const paramschema::SnapshotResult& legacy_snapshot)
{
  const paramschema::SnapshotRestoreResult restored =
      paramschema::restore_snapshot_state(paramschema::mesher_schema(), legacy_snapshot);
  print_lines(stderr, "parameter compare warning: ", restored.warnings);
  if (!restored.errors.empty()) {
    print_lines(stderr, "parameter compare error: ", restored.errors);
    return false;
  }
  return true;
}
 
bool run_parser_legacy_compare(const std::vector<std::string>& runtime_args,
                               const LegacyCompareInput&       legacy_input,
                               const RuntimeCompatOptions&     compat)
{
  const ParserCompareMode mode = parser_compare_mode();
  if (mode == ParserCompareMode::Off) {
    return true;
  }
  const bool fallback_to_legacy = legacy_fallback_requested(compat);
 
  if (runtime_args.empty()) {
    std::fprintf(stderr, "parameter compare error: unable to reconstruct mesher argv\n");
    if (fallback_to_legacy) {
      print_legacy_fallback_workaround();
    }
    return mode != ParserCompareMode::Strict;
  }
 
  paramschema::SnapshotResult parser_snapshot;
  if (!capture_current_parser_snapshot(&parser_snapshot)) {
    std::fprintf(stderr, "parameter compare error: unable to snapshot the parser runtime state\n");
    if (fallback_to_legacy) {
      print_legacy_fallback_workaround();
    }
    return mode != ParserCompareMode::Strict;
  }
  maybe_force_test_mismatch(&parser_snapshot);
 
  std::string helper_path;
  if (!resolve_legacy_snapshot_helper(runtime_args[0], &helper_path)) {
    std::fprintf(stderr, "parameter compare error: unable to locate mesher-param-parser-legacy-snapshot\n");
    if (fallback_to_legacy) {
      print_legacy_fallback_workaround();
    } else {
      std::fprintf(stderr,
                   "parameter compare error: rerun with OPENCARP_PARAM_COMPARE=0 to bypass this temporary compatibility gate\n");
    }
    return mode != ParserCompareMode::Strict;
  }
 
  if (!legacy_input.available) {
    std::fprintf(stderr, "parameter compare warning: %s\n", legacy_input.unavailable_reason.c_str());
    if (fallback_to_legacy) {
      std::fprintf(stderr,
                   "parameter compare warning: legacy fallback is unavailable because no legacy baseline exists for this input\n");
    }
    return true;
  }
 
  std::string snapshot_path;
  if (!create_temp_output_path("", &snapshot_path)) {
    std::fprintf(stderr, "parameter compare error: unable to create temporary snapshot file\n");
    if (fallback_to_legacy) {
      print_legacy_fallback_workaround();
    } else {
      std::fprintf(stderr,
                   "parameter compare error: rerun with OPENCARP_PARAM_COMPARE=0 to bypass this temporary compatibility gate\n");
    }
    return mode != ParserCompareMode::Strict;
  }
 
  const bool helper_ok = run_legacy_snapshot_helper(helper_path, legacy_input.runtime_args, snapshot_path);
  paramschema::SnapshotResult legacy_snapshot;
  std::string                 read_error;
  const bool loaded = helper_ok &&
                      paramschema::snapshotio::read_snapshot_file(snapshot_path, &legacy_snapshot, &read_error);
  unlink(snapshot_path.c_str());
 
  if (!loaded) {
    if (!read_error.empty()) {
      std::fprintf(stderr, "parameter compare error: %s\n", read_error.c_str());
    } else {
      std::fprintf(stderr, "parameter compare error: unable to capture the legacy parameter state\n");
    }
    if (fallback_to_legacy) {
      print_legacy_fallback_workaround();
    } else {
      std::fprintf(stderr,
                   "parameter compare error: rerun with OPENCARP_PARAM_COMPARE=0 to bypass this temporary compatibility gate\n");
    }
    return mode != ParserCompareMode::Strict;
  }
 
  const paramschema::SnapshotComparisonResult comparison =
      paramschema::compare_snapshot_results(paramschema::mesher_schema(), parser_snapshot, legacy_snapshot);
  if (!comparison.errors.empty() || !comparison.mismatches.empty()) {
    print_lines(stderr, "", paramschema::format_snapshot_comparison_report(comparison, "parser compare"));
    std::fprintf(stderr,
                 "parameter compare error: the parser runtime and legacy param() produced different parameter states\n");
    std::fprintf(stderr,
                 "parameter compare error: please open an issue and include the triggering command line and parameter files\n");
    if (fallback_to_legacy) {
      if (!restore_legacy_snapshot_state(legacy_snapshot)) {
        print_legacy_fallback_workaround();
        return false;
      }
      std::fprintf(stderr, "parameter compare warning: continuing with the legacy parameter state\n");
      return true;
    }
 
    print_legacy_fallback_workaround();
    return mode != ParserCompareMode::Strict;
  }
 
  return true;
}
 
}  // namespace
 
 
int main( int argc, char *argv[] )
{
  LegacyCompareInput        legacy_input;
  std::vector<std::string> normalized_args;
  RuntimeCompatOptions     compat;
  std::string              normalize_error;
  if (!build_legacy_compare_input(argc, argv, &legacy_input, &normalize_error) ||
      !normalize_runtime_args(argc, argv, &normalized_args, &compat, &normalize_error)) {
    std::fprintf(stderr, "\n*** %s\n\n", normalize_error.c_str());
    exit(EXIT_FAILURE);
  }
 
  paramschema::ExecutionResult executed;
  if (!apply_parser_runtime_args(normalized_args, true, &executed)) {
    exit(EXIT_FAILURE);
  }
 
  if (legacy_input.available) {
    for (std::size_t i = 0; i < executed.validation.assignments.size(); ++i) {
      if (!executed.validation.assignments[i].synthesized) {
        continue;
      }
      legacy_input.available = false;
      legacy_input.runtime_args.clear();
      legacy_input.runtime_args.push_back(normalized_args[0]);
      legacy_input.unavailable_reason =
          "legacy compare is unavailable because the parser inferred optional controller counts from the original input";
      break;
    }
  }
 
  if (!run_parser_legacy_compare(normalized_args, legacy_input, compat)) {
    exit(EXIT_FAILURE);
  }
 
  float tsize[3], x0[3];
  bool  symbath[3];
 
 
  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;
}