#include <iostream.h>
#include <stdio.h>
#include <iomanip.h>
#include <math.h>
#include <vector>
#include <map>
#include <deque>

/*******************************************************************************************
*
*                             T H R E E - D I M E N S I O N A L
*                             ---------------------------------
*        F I N I T E   E L E M E N T   P A R T I T I O N I N G   B Y   A   P L A N E
*        ---------------------------------------------------------------------------
*
*                         N. Sukumar, University of California, Davis
*
*                                   C + +   C O D E
*                                   ---------------
*
* Purpose
* =======
* Partition a 3-dimensional finite element (tetrahedron or hexahedron) by a plane that is
* given by a point on it and its bi-normal (m) [m is a vector that is normal to the (crack) 
* plane].  The nodal connectivities of the faces of the finite element are stored in 
* counter-clockwise orientation (ccw), and the face connectivities are also stored such that 
* an edge that is common to two faces is traversed in opposite directions. The nodal and face 
* connectivity of a tetrahedron are such that if the nodes on a face are traversed in a 
* certain sense, then the thumb, using the right-hand ruls, points in the direction of the 
* fouth node. The finite element is partitioned into 4-noded tetrahedral elements on either 
* side of the (crack) plane. Degenerate cases are handled using a tolerance in the partitioning 
* algorithm. This function is used in the application of the eXtended Finite Element Method 
* (X-FEM) to crack modeling in 3-dimensions. A paper, which describes the partitioning 
* algorithm can be found at http://dilbert.engr.ucdavis.edu/~suku/xfem/papers/xfem_cracks3d.pdf

* Author
* ======
* N. Sukumar, University of California, Davis
* E-mail: n-sukumar@ucdavis.edu
* First version: June 1999; Updated: August 2003
*
* Reference
* =========
* N. Sukumar, `Element partitioning code in 2-d and 3-d for the extended finite
* element method,' Available from http://dilbert.engr.ucdavis.edu/~suku/xfem, 2000
*
* Compilation
* ===========
* The code has been tested on Red Hat Linux 8.0 with gcc versions 2.95.3
* 20010315 release as well as gcc version 3.2 20020903 release. For
* compiling, use
*
* 1. g++ -Wno-deprecated -DVERTWO -DIxTETRA4  main.C [tetrahedral element], or
* 2. g++ -Wno-deprecated -DVERTHREE -DIxHEX8  main.C [hexahedral element], 
* where the macro VERTWO or VERTHREE is used to indicate the gcc version and
* a.out is the executable in either case (main.C is this file)
*
* Inputs
* ======
* The input parameters and settings are indicated in the code itself. The finite element
* is an instance of the class IxElement_c and the given plane (such as a crack plane)
* is that of the class GxElement_c in the X-FEM code.  The number of nodes, faces, maximum
* node and element IDs, connectivity information for the element and the faces are set for the
* tetrahedral (IxTETRA4) as well as the hexahedral (IxHEX8) element to test this code
*
* Output
* ======
* The graphics plotting package TECPLOT [Amtec Inc.] is used to display the results. The
* faces of the finite element are output in Tecplot format to the file "faces.dat" and
* the partitioned elements are output in Tecplot format to the file "subelements3d.dat"
*
* Algorithm
* =========
* A short description of the algorithm follows. The orientation of the nodes of the finite
* element are first set (ABOVE = +1, BELOW = -1, or ON = 0) using the plane (GxElement_c)
* that partitions the element. Partitioning is carried out if and only if there exists two 
* nodes with orientation = 1 and orientation = -1 in the connectivity. A loop over the faces 
* of the finite element is carried out, and if there is an intersection for an edge belonging 
* to a face, the point of intersection is computed and added to the node and face-edge
* list. The orientation of the point of interesection is +10 if the orientation of the first 
* point of the edge is 1, and -10 if the orientation of the first point is -1. This updates the 
* face lists of the finite element. Now, using the face list of the finite element, two new 
* face lists for the domains ABOVE and BELOW the (crack) plane are created.  The two faces
* (upper and lower) that correspond to the (crack) plane are also created such the orientation
* of the nodes in the face is per requirement - ccw for nodes in the ABOVE surface and cw
* orientation for the nodes in the BELOW surface preserve the correct orientation of the
* partitioned tetrahedrons. Now, the centroid of the faces is evaluated and the orientation of 
* the new point (node) is set to 100 for the faces in the ABOVE list and -100 for the faces in 
* BELOW list. The centroid of the two volumes is evaluated with orientation = 1000 for the 
* ABOVE volume and orientation = -1000 for BELOW volume. A loop over the two face lists is 
* carried out and using the centroid for each face and the centroid of the two volume domains, 
* the partitioned elements are obtained
*
* Disclaimer/Questions
* ====================
* This code is in the public-domain and is provided as is with no guarantees whatsoever. If 
* you do have any comments and/or spot any bugs/improvements, please do contact me. If you 
* find this code of use in your application, I'd appreciate an e-mail with references to 
* your work(s)
* 
*******************************************************************************************/

#define FALSE 0

#define TRUE 1

const double TOLERANCE = 1.e-6;

// Functions
// =========
// dot product of two vectors
// double vecdotvec(double vec1[3], double vec2[3]); 
//
// L_2-norm of a vector
// double norm(double vec[3]);
//
// cross product of two vectors
// void vecXvec(double a[3], double b[3], double res[3]);
//
// tests if a point is ON, ABOVE, or BELOW the given (crack) plane
// int OnAboveOrBelow(double px, double py, double pz, double p[3], double m[3]);
//
// determines the intersection of the given (crack) plane with the line joining p1 to p2
// int LineIntersection(double p1x, double p1y, double p1z, 
//                      double p2x, double p2y, double p2z,
//                      double& px, double& py, double& pz, double p0[3], double m[3]);

double vecdotvec(double vec1[3], double vec2[3]) {
  double res= 0.0;
  for(int k=0;k<3;k++){
    res += vec1[k] * vec2[k];
  }
  return res;
}

double norm(double vec[3]) {
  double mod = sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
  return mod;
}

void vecXvec(double a[3], double b[3], double res[3]) {
  res[0] = a[1] * b[2] - a[2] * b[1];
  res[1] = a[2] * b[0] - a[0] * b[2];
  res[2] = a[0] * b[1] - a[1] * b[0];
  return;
}

inline int OnAboveOrBelow(double px, double py, double pz, double p[3], double m[3]) {
// line equation in 2-D and that of the plane in 3-D 
  double planeEq = m[0]*(px - p[0]) + m[1]*(py - p[1]) + m[2]*(pz - p[2]);
  
  if (fabs(planeEq) <= TOLERANCE) return 0; 
  else if (planeEq > TOLERANCE)   return 1;
  else                            return -1; 
}

int LineIntersection(double p1x, double p1y, double p1z, 
                     double p2x, double p2y, double p2z,
                     double& px, double& py, double& pz, double p0[3], double m[3]) {

// point of intersection of a line (p1,p2) and a GxElement_c in 3-D;
// equation of the plane: m . (p - p0) = 0, where p0 is a given point on the plane
// and m is the normal to the plane;
// parametric representation of a point on a line (pa,pb) is: pa + t*(pb - pa);
// if there is an intersection, then the point lines on the line as well as on the
// plane, whence t = [0,1]

  double denominator = m[0]*(p1x - p2x) + m[1]*(p1y - p2y) + m[2]*(p1z - p2z);

  if (fabs(denominator) < TOLERANCE) {
    cout << "\n Error: Line is either parallel to the plane or lies in the plane\n";
    assert (0);
  }

  double t = (m[0]*(p1x - p0[0]) + m[1]*(p1y - p0[1]) + m[2]*(p1z - p0[2]))/denominator;
 
  if (t < 0.0 || t > 1.0) 
    return FALSE;
  else {
    px = p1x + t*(p2x - p1x);
    py = p1y + t*(p2y - p1y);
    pz = p1z + t*(p2z - p1z);
  }

  return TRUE;
}

// using namepspace std;

int main() {

using namespace std; // std namespace

// Partitioning a 3D element that is cut by a (crack) plane into sub-elements (tetrahedrons) 
// for the purpose of numerical integration 
 
  int nbNodes;                // no of nodes for the (IxElement_c) element

  int nbFaces;                // no of faces of the IxElement_c

  int nbNodesFace;            // no of nodes for a face of the IxElement_c

  int maxFaceId;              // max face Id in the IxElement_c (local max.). IDs of
                              // the faces are locally numbered as 0,1,...,Nbfaces-1

  int maxNodeId;             // max node Id

  int maxElemId;             // max element Id

#ifdef IxTETRA4
  nbNodes = 4; 
  nbFaces = 4; 
  nbNodesFace = 3; 
  maxElemId = 1; 
  maxFaceId = 4;
  maxNodeId = 4;
#elif defined(IxHEX8)
  nbNodes = 8; 
  nbFaces = 6; 
  nbNodesFace = 4; 
  maxElemId = 1; 
  maxFaceId = 6;
  maxNodeId = 8;
#endif

// vector to store the element connectivity
  std::vector<int> connect;

// map definition, typedef, and const iterators for face ID to face connectivity
  std::map<int,std::vector<int>,std::less<int> >  Face;
  typedef std::map<int,std::vector<int>,std::less<int> >::value_type face_add;
  std::map<int,std::vector<int>,std::less<int> >::const_iterator itF, itFEnd;
  std::vector<int> facenodes; // nodal connectivity for the face

// map definition, typedef and const iterators for node ID to nodal coordinates
  std::map<int,std::vector<double>,std::less<int> > Coord;
  typedef std::map<int,std::vector<double>,std::less<int> >::value_type coord_add;
  std::map<int,std::vector<double>,std::less<int> >::const_iterator itC, itCEnd;

// map definition, const iterators, and typedef for insertion for node ID to orientation
  std::map<int,int,std::less<int> >  Abo; 
  typedef std::map<int,int,std::less<int> >::value_type abo_add;
  std::map<int,int,std::less<int> >::const_iterator itO, itOEnd;
  std::map<int,int,std::less<int> >::const_iterator itA;
  std::map<int,int,std::less<int> >::const_iterator itA1;
  std::map<int,int,std::less<int> >::const_iterator itA2;

// Initializations

  int i, j, k;
  std::vector<double> xyz;

// node Id to nodal coordinates map

#ifdef IxTETRA4
  cout << "\n Partitioning a IxTETRA4 which is cut by a (crack) plane\n";
  double x[4] = {0.,1.,0.,0.};
  double y[4] = {0.,0.,1.,0.};
  double z[4] = {0.,0.,0.,1.};
  int elemconnect[4] = {1,2,3,4};

  int v[4][3] = { {0,2,1}, {0,3,2}, {1,2,3}, {0,1,3} };
#elif defined(IxHEX8)
  cout << "\n Partitioning a IxHEX8 which is cut by a (crack) plane\n";
  double x[8] = {0.,1.,1.,0.,0.,1.,1.,0.};
  double y[8] = {0.,0.,1.,1.,0.,0.,1.,1.};
  double z[8] = {0.,0.,0.,0.,1.,1.,1.,1.};
  int elemconnect[8] = {1,2,3,4,5,6,7,8};

  int v[6][4] = { {0,1,2,3}, {4,7,6,5}, {0,4,5,1}, {1,5,6,2}, {3,2,6,7}, {0,3,7,4} };
#endif
  
// node Id to nodal coordinates map
  for (i = 0; i < nbNodes; ++i) {
    xyz.clear();
    xyz.push_back(x[i]); 
    xyz.push_back(y[i]); 
    xyz.push_back(z[i]); 
    Coord.insert(coord_add(i+1,xyz)); // start IDs from 1
  }
 
// element connectivity
  for (i = 0; i < nbNodes; ++i) {
    connect.push_back(elemconnect[i]);
  }

// face connectivity
  for (i = 0; i < nbFaces; ++i) {
    for (j = 0; j < nbNodesFace; ++j) {
      facenodes.push_back(elemconnect[ v[i][j] ]);
    }
    Face.insert(face_add(i+1,facenodes)); // start IDs from 1
    facenodes.clear();
  }

// plane is given by its bi-normal m and a point p

   double m[3] = {0.0, -0.5, 0.5};
//  double m[3] = {-1., 1., 0.};
//  double m[3] = {1.,0.,0.};

   double p[3] = {0.0, 0.0, 0.0};
//  double p[3] = {0., 0., 0.};
//  double p[3] = {0.999999,0.,0.};

// set the node Id to orientation map: if the node is above the crack (GxElement_c) element,
// orientation = 1, if the node is below the crack (GxElement_c) element, orientation = -1,
// and if the node is on the crack (GxElement_c) element, orientation = 0

  double px, py, pz;
  int orient;

  for (i = 0; i < nbNodes; ++i) {
    itC = Coord.find(connect[i]);
    xyz = itC->second;
    px = xyz[0]; py = xyz[1]; pz = xyz[2];
    orient = OnAboveOrBelow(px,py,pz,p,m);
    Abo.insert(abo_add(connect[i],orient));
    cout << "\n Node no. " << connect[i] << " and coords " << px << " " <<  py << " " << pz;
    cout << "\n Node is " << (orient==-1 ? "BELOW" : ((orient==1) ? "ABOVE" : "ON")) << " the crack " << endl;
  }

// check that the nodal orientation of at least two nodes is +1 and -1 in the
// connectivity list of the finite element; if not, then the crack plane does
// not cut the element or is along an element facet. In either case, no
// partitioning of the elements into sub-elements is required

  bool partition = false, nodeAbove = false, nodeBelow = false;
  itO = Abo.begin();
  itOEnd = Abo.end(); 
  for (; itO != itOEnd; ++itO) {
    if (itO->second == -1) nodeBelow = true;
    if (itO->second == 1) nodeAbove = true;
  }
  if (nodeBelow && nodeAbove) partition = true;

// if partition = true, at least one face is cut by the crack plane; if not, then
// no partitioning of the element is required and hence the function returns

  if (!partition) {
     cout << "\nELEMENT PARTIONING NOT REQUIRED\n";
     return 1;
  }

// loop over the faces to check faces that are cut by the crack

  int intersect;
  bool isCut;
  int orient1, orient2;
  int orient_compare;
  double p1x, p1y, p1z, p2x, p2y, p2z;

  std::map<int,std::vector<int>,std::less<int> > NewFaceConn;
  std::vector<int> NodesPoints;
  typedef std::map<int,std::vector<int>,std::less<int> >::value_type new_add;
  std::map<int,std::vector<int>,std::less<int> >::const_iterator itN, itNEnd;

  itF = Face.begin();
  itFEnd = Face.end();
  for (; itF != itFEnd; ++itF) {
    facenodes = itF->second;
    isCut = false;
    orient_compare = 0;
    nbNodesFace = facenodes.size();
    for (j = 0; j < nbNodesFace; ++j) {
      itA = Abo.find(facenodes[j]);
      orient = itA->second;
      if ((orient != 0) && (!orient_compare)) orient_compare = orient;
      if (orient*orient_compare == -1) isCut = true;
    }
    cout << "\n Face " << itF->first << " is " << (isCut ? "" : "NOT ") << "CUT by the crack ";   

// for faces cut by the crack, check if contiguous nodes in the nodal connectivity of the
// face have +1 and -1 (or vice-versa) for the orientations 
   
    if (isCut) {

      cout << "\n ****************************************\n" ;

// at least one face is cut by the crack plane
      for (j = 0; j < nbNodesFace; ++j) {
        k = (j == nbNodesFace-1) ? 0 : j+1;
        itA1 = Abo.find(facenodes[j]);  orient1 = itA1->second;
        itA2 = Abo.find(facenodes[k]); orient2 = itA2->second;

        if (orient1*orient2 == -1) {
// find the intersection and insert new points into the connectivity of the face
          itC = Coord.find(facenodes[j]);
          xyz = itC->second;         
          p1x = xyz[0]; p1y = xyz[1]; p1z = xyz[2]; 
          itC = Coord.find(facenodes[k]);
          xyz = itC->second;         
          p2x = xyz[0]; p2y = xyz[1]; p2z = xyz[2];

          cout << "\n j is: " << j << " : " << p1x << setw(5) << p1y << setw(5) << p1z;
          cout << " and " << p2x << setw(5) << p2y << setw(5) << p2z;

// px, py and pz are the coordinates of the intersection point
          intersect = LineIntersection(p1x,p1y,p1z,p2x,p2y,p2z,px,py,pz,p,m);
          if (!intersect) assert(0);
          cout << "\n Intersection point is: " << px << " " << py << " " << pz;

          xyz.clear();
     	    xyz.push_back(px); xyz.push_back(py); xyz.push_back(pz); 

          maxNodeId++;

// insert the new point in the node/point orientation map: if the first node in the
// connectivity for the point is above the crack plane, insert +10 for the orientation
// and if the first node in the connectivity is below the crack plane, insert
// -10 for the orientation
          orient = (orient1 == -1) ? -10 : 10;
          Abo.insert(abo_add(maxNodeId,orient));

// insert the new point coordinates in the node-coordinates map
          Coord.insert(coord_add(maxNodeId,xyz));

// insert previous node and the point (ID is maxNodeId) in the NEW connectivity list                      
          NodesPoints.push_back(facenodes[j]);
          NodesPoints.push_back(maxNodeId);
        }
        else {
          NodesPoints.push_back(facenodes[j]);
        }
      }

    }
     
    if (!isCut) {
// copy the face to the new face (no intersections)
      for (j = 0; j < nbNodesFace; ++j) NodesPoints.push_back(facenodes[j]);
    }

    NewFaceConn.insert(new_add(itF->first,NodesPoints)); // itF->face ID

    NodesPoints.clear();

    cout << "\n\n **********************************************************\n" ;
  }

// create the face (above and below the crack plane) connectivity for the nodes/points 
// along the edges of the face

  std::map<int,std::vector<int>,std::less<int> > FaceAboveCrack;
  std::map<int,std::vector<int>,std::less<int> > FaceBelowCrack;

  itN = NewFaceConn.begin();
  itNEnd = NewFaceConn.end();

  std::vector<int> above;
  std::vector<int> below;

  for (; itN != itNEnd; ++itN) {
    for (i = 0; i < itN->second.size(); ++i) {
      itO = Abo.find(itN->second[i]);
      if (itO->second == -1 || itO->second == 0 || itO->second == -10 || itO->second == 10) {
//   node/point of the face is below the crack plane
        below.push_back(itN->second[i]);
      } 
      if (itO->second == 1 || itO->second == 0 || itO->second == -10 || itO->second == 10) {
//   node/point of the face is above the crack plane
        above.push_back(itN->second[i]);
      }
    }

    if (below.size() > 2) FaceBelowCrack.insert(new_add(itN->first,below));
    if (above.size() > 2) FaceAboveCrack.insert(new_add(itN->first,above));

    below.clear();
    above.clear();
  }

// creating the face ID and the connectivity map for the nodes/points on the
// crack plane for the top and bottom surfaces

#if VERTWO
  std::vector<int> unsorted;
  std::vector<int>::iterator pos;
#elif defined(VERTHREE)
  std::deque<int> unsorted;
  std::deque<int>::iterator pos;
#endif

  std::vector<int> AbovePointId;
  std::vector<int> BelowPointId;

  double vec1[3], vec2[3], res[3];

  double sign, tmin, tmax;

  double cost, theta;

  int index;

  double above_centroid[3] = {0.0, 0.0, 0.0};

// face OnAbove the crack plane --- BEGIN

  itO = Abo.begin();
  itOEnd = Abo.end(); 
  j = 0;

  for (; itO != itOEnd; ++itO) {
    itC = Coord.find(itO->first);     
    xyz = itC->second;
    if ((itO->second == 0) || (itO->second == 10)) {
      unsorted.push_back(itO->first);
      itC = Coord.find(itO->first);
      xyz = itC->second;
      for (i = 0; i < 3; ++i) above_centroid[i] += xyz[i];
      ++j;
    }
  }

  for (k = 0; k < 3; ++k) above_centroid[k] /= (double) j;
  ++maxNodeId;
  AbovePointId.push_back(maxNodeId);
  orient = 100;
  Abo.insert(abo_add(maxNodeId,orient)); 

  std::vector<double> centroid;
  for (int ii = 0; ii < 3; ++ii) centroid.push_back(above_centroid[ii]);
  Coord.insert(coord_add(maxNodeId,centroid));

// sort the nodes on the OnAbove surface of the crack so that
// the connectivity is correct -> a ccw connectivity of the three 
// nodes in the plane (triangle) results in a +ve sense (right-hand 
// screw opens towards the the fourth node): the four nodes define 
// a tetrahedron in R^3

  index = 0;
  while (!unsorted.empty()) {
    tmin = acos(-1.0);
    itC = Coord.find(unsorted[index]);
    above.push_back(unsorted[index]);
    xyz = itC->second; 
#if VERTWO
    pos = &unsorted[index];
#elif defined(VERTHREE)
    pos = find(unsorted.begin(), unsorted.end(), unsorted[index]);
#endif
    unsorted.erase(pos);

    for (j = 0; j < 3; ++j) { vec1[j] = xyz[j] - above_centroid[j]; }

    for (i = 0; i < unsorted.size(); ++i) {
      itC = Coord.find(unsorted[i]);
      xyz = itC->second;
      for (j = 0; j < 3; ++j) vec2[j] = xyz[j] - above_centroid[j];
      vecXvec(vec1,vec2,res);
      sign = ( vecdotvec(res,m) >= 0.0 ) ? 1.0 : -1.0;
      cost = vecdotvec(vec1,vec2)/(norm(vec1)*norm(vec2));
      theta = sign*(acos((float)cost));
      if (theta > 0 && theta < tmin) { tmin = theta; index = i; }
    }
  }

// create the face OnAbove the crack plane using the sorted points from 
// vector above

  FaceAboveCrack.insert(new_add(++maxFaceId,above)); above.clear();

// face OnAbove the crack plane --- END

// face OnBelow the crack plane --- BEGIN

  double below_centroid[3] = {0.0, 0.0, 0.0};

  itO = Abo.begin();
  itOEnd = Abo.end(); 
  j = 0;
   for (; itO != itOEnd; ++itO) {
    itC = Coord.find(itO->first);     
    xyz = itC->second;
    if ((itO->second == 0) || (itO->second == -10)) {
      unsorted.push_back(itO->first);
      itC = Coord.find(itO->first);
      xyz = itC->second;
      for (i = 0; i < 3; ++i) below_centroid[i] += xyz[i];
      ++j;
    }
  }

  for (k = 0; k < 3; ++k) below_centroid[k] /= (double) j;
  ++maxNodeId;
  BelowPointId.push_back(maxNodeId);
  orient = -100;
  Abo.insert(abo_add(maxNodeId,orient)); 

  centroid.clear();
  for (int ii = 0; ii < 3; ++ii) centroid.push_back(below_centroid[ii]);
  Coord.insert(coord_add(maxNodeId,centroid));

// sort the nodes on the OnBelow surface of the crack so that
// the connectivity is correct -> a cw connectivity of the three 
// nodes in the plane (triangle) results in a +ve sense (right-hand 
// screw opens towards the the fourth node): the four nodes define 
// a tetrahedron in R^3

  index = 0;
  while (!unsorted.empty()) {
    tmax = -acos(-1.0);
    itC = Coord.find(unsorted[index]);
    below.push_back(unsorted[index]);
    xyz = itC->second; 
#if VERTWO
    pos = &unsorted[index];
#elif defined(VERTHREE)
    pos = find(unsorted.begin(), unsorted.end(), unsorted[index]);
#endif
    unsorted.erase(pos);

    for (j = 0; j < 3; ++j) { vec1[j] = xyz[j] - below_centroid[j]; }

    for (i = 0; i < unsorted.size(); ++i) {
      itC = Coord.find(unsorted[i]);
      xyz = itC->second;
      for (j = 0; j < 3; ++j) vec2[j] = xyz[j] - below_centroid[j];
      vecXvec(vec1,vec2,res);
      sign = ( vecdotvec(res,m) >= 0.0 ) ? 1.0 : -1.0;
      cost = vecdotvec(vec1,vec2)/(norm(vec1)*norm(vec2));
      theta = sign*(acos((float)cost));
      if (theta < 0 && theta > tmax) { tmax = theta; index = i; }
    }
  }

// create the face OnBelow the crack plane using the sorted points from 
// vector below

// maxFaceId has been updated in the case for OnAbove and the same Id
// is used here
  FaceBelowCrack.insert(new_add(maxFaceId,below)); below.clear();

// face OnBelow the crack plane --- END

// create the centroid of the coordinates of the faces for facet domains
// below and above the crack plane. Enter the centroid point in the coord and 
// orientation maps and also in the connectivity of the faces

  std::map<int,std::vector<int>,std::less<int> >::const_iterator itFace, itFaceEnd;

  itFace = FaceAboveCrack.begin();
  itFaceEnd = FaceAboveCrack.end(); 
  --itFaceEnd; // - 1 because the centroid of the face OnAbovePlane has 
               //  already been set

  for (; itFace != itFaceEnd; ++itFace) {
    above_centroid[0] = above_centroid[1] = above_centroid[2] = 0.0;
    for (i = 0; i < itFace->second.size(); ++i) {
      itC = Coord.find(itFace->second[i]);
      xyz = itC->second;
      for (j = 0; j < 3; ++j) { above_centroid[j] += xyz[j]; }
    }

// update coordinates and orientation maps (orientation +100 for point above
// the crack plane)
    for (j = 0; j < 3; ++j) { above_centroid[j] /= double (itFace->second.size()); }
    orient = 100; maxNodeId++; 
    AbovePointId.push_back(maxNodeId);
    Abo.insert(abo_add(maxNodeId,orient));
    centroid.clear();
    for (int ii = 0; ii < 3; ++ii) centroid.push_back(above_centroid[ii]);
    Coord.insert(coord_add(maxNodeId,centroid));
  }
 
  itFace = FaceBelowCrack.begin();
  itFaceEnd = FaceBelowCrack.end();
  --itFaceEnd;  // - 1 because the centroid of the face OnBelowPlane has 
                // already been set

  for (; itFace != itFaceEnd; ++itFace) {
    below_centroid[0] = below_centroid[1] = below_centroid[2] = 0.0;
    for (i = 0; i < itFace->second.size(); ++i) {
      itC = Coord.find(itFace->second[i]);
      xyz = itC->second;
      for (j = 0; j < 3; ++j) { below_centroid[j] += xyz[j]; }
    }

// update coordinates and orientation maps (orientation -100 for point below 
// the crack plane)
    for (j = 0; j < 3; ++j) { below_centroid[j] /= double (itFace->second.size()); }
    orient = -100; maxNodeId++; 
    BelowPointId.push_back(maxNodeId);
    Abo.insert(abo_add(maxNodeId,orient)); 
    centroid.clear();
    for (int ii = 0; ii < 3; ++ii) centroid.push_back(below_centroid[ii]);
    Coord.insert(coord_add(maxNodeId,centroid));
  }

// create the centroid of the coordinates of the points/nodes below and above
// the crack plane - the element is partitioned into tetrahedrons above and 
// below the crack plane

  below_centroid[0] = below_centroid[1] = below_centroid[2] = 0.;
  above_centroid[0] = above_centroid[1] = above_centroid[2] = 0.;

  itO = Abo.begin();
  itOEnd = Abo.end();
  int nbabove = 0, nbbelow = 0;
  for (; itO != itOEnd; ++itO) {
    itC = Coord.find(itO->first);     
    xyz = itC->second;
    if ((itO->second == -1) || (itO->second == 0) || (itO->second == -10)) {
      for (i = 0, ++nbbelow; i < 3; ++i) below_centroid[i] += xyz[i]; 
    }
    if ((itO->second == 1) || (itO->second == 0) || (itO->second == 10)) {
      for (i = 0, ++nbabove; i < 3; ++i) above_centroid[i] += xyz[i]; 
    }
  }

  for (i = 0; i < 3; ++i) { 
    below_centroid[i] /= (double) nbbelow;
    above_centroid[i] /= (double) nbabove;
  }

// insert centroid into the coordinates and orientation maps (orientation +1000 and -1000
// for points associated to the element volume above and below the crack plane,
// respectively)

  orient = 1000; maxNodeId++; 
  int aboveVolPoint = maxNodeId;
  Abo.insert(abo_add(maxNodeId,orient)); 
  centroid.clear();
  for (int ii = 0; ii < 3; ++ii) centroid.push_back(above_centroid[ii]);
  Coord.insert(coord_add(maxNodeId,centroid));

  orient = -1000; maxNodeId++; 
  int belowVolPoint = maxNodeId;
  Abo.insert(abo_add(maxNodeId,orient)); 
  centroid.clear();
  for (int ii = 0; ii < 3; ++ii) centroid.push_back(below_centroid[ii]);
  Coord.insert(coord_add(maxNodeId,centroid));

// tetrahedral elements 
// map SubElements contains the element IDs and their connectivities; all sub-elements
// that are a partition of the IxElement_c are contained in the map

  std::map<int,std::vector<int>,std::less<int> > SubElements;
  typedef std::map<int,std::vector<int>,std::less<int> >::value_type sub_add;
  int subconn[4];
  std::vector<int> connvec;


  itN = FaceAboveCrack.begin();
  itNEnd = FaceAboveCrack.end();
  for (k=0; itN != itNEnd; ++itN,++k) {
    index = (k == (FaceAboveCrack.size()-1)) ? 0 : k+1;
    for (i = 0; i < itN->second.size(); ++i) {
      j = (i == (itN->second.size()-1)) ? 0 : i+1; 
      subconn[0] = itN->second[i];
      subconn[1] = itN->second[j];
      subconn[2] = AbovePointId[index];
      subconn[3] = aboveVolPoint;
      ++maxElemId;
      connvec.clear();
      for (int ii = 0; ii < 4; ++ii) connvec.push_back(subconn[ii]);
      SubElements.insert(sub_add(maxElemId,connvec));
    }
  }
  itN = FaceBelowCrack.begin();
  itNEnd = FaceBelowCrack.end();
  for (k=0; itN != itNEnd; ++itN,++k) {
    index = (k == (FaceBelowCrack.size()-1)) ? 0 : k+1;
    for (i = 0; i < itN->second.size(); ++i) {
      j = (i == (itN->second.size()-1)) ? 0 : i+1; 
      subconn[0] = itN->second[i];
      subconn[1] = itN->second[j];
      subconn[2] = BelowPointId[index];
      subconn[3] = belowVolPoint;
      ++maxElemId;
      connvec.clear();
      for (int ii = 0; ii < 4; ++ii) connvec.push_back(subconn[ii]);
      SubElements.insert(sub_add(maxElemId,connvec));
    }
  }

  cout << "\n Size of new face to connectivity (nodes/points) map is: " << NewFaceConn.size() << endl;

  itN = NewFaceConn.begin();
  itNEnd = NewFaceConn.end();
  for (; itN != itNEnd; ++itN) {
    cout << "\n Face is: " << itN->first ; 
    for (k = 0; k < itN->second.size(); ++k) 
      cout << "\n Connectivity k is " << k << " and node/point is: " << itN->second[k];
  }

  cout << "\n\n **********************************************************\n" ;

  cout << "\n Size of FaceAboveCrack to connectivity (nodes/points) map is: " << FaceAboveCrack.size() << endl;
 
  itN = FaceAboveCrack.begin();
  itNEnd = FaceAboveCrack.end();
  for (; itN != itNEnd; ++itN) {
    cout << "\n Face is: " << itN->first ; 
    for (k = 0; k < itN->second.size(); ++k) 
      cout << "\n Connectivity k is " << k << " and node/point is: " << itN->second[k];
  }

  cout << "\n\n **********************************************************\n" ;

  cout << "\n Size of FaceBelowCrack to connectivity (nodes/points) map is: " << FaceBelowCrack.size() << endl;

  itN = FaceBelowCrack.begin();
  itNEnd = FaceBelowCrack.end();
  for (; itN != itNEnd; ++itN) {
    cout << "\n Face is: " << itN->first ; 
    for (k = 0; k < itN->second.size(); ++k) 
      cout << "\n Connectivity k is " << k << " and node/point is: " << itN->second[k];
  }

  cout << "\n\n **********************************************************\n" ;

  cout << "\n Size of Node/Point to Coordinates map is: " << Coord.size() << endl;
  itC = Coord.begin();
  itCEnd = Coord.end();
  for (; itC != itCEnd; ++itC) {
    cout << "\n Node/Point is: " << itC->first << " and coords are: " << itC->second[0] << " " <<
         itC->second[1] << " " << itC->second[2];
  }

  cout << "\n\n **********************************************************\n" ;
 
  cout << "\n Size of Node/Point to Orientation map is: " << Abo.size() << endl;
  itO = Abo.begin();
  itOEnd = Abo.end();
  for (; itO != itOEnd; ++itO) {
    cout << "\n Node/Point is: " << itO->first << " and its orientation is: " << itO->second;
  }

  cout << "\n\n **********************************************************\n" ;


//TECPLOT OUTPUT

  FILE* OUTPUT;
  OUTPUT = fopen("subelements3d.dat","w");
  fprintf(OUTPUT, " TITLE = \"Example: Tetrahedral Sub-Element Data\" "); fprintf(OUTPUT,"\n");
  fprintf(OUTPUT, " VARIABLES = \"X\" ");
  fprintf(OUTPUT, ", \"Y\" ");
  fprintf(OUTPUT, ", \"Z\" "); fprintf(OUTPUT,"\n");

  fprintf(OUTPUT, " ZONE N=%d",Coord.size()); 
  fprintf(OUTPUT, ", E=%d",SubElements.size());
  fprintf(OUTPUT, ", F=FEPOINT");
  fprintf(OUTPUT, ", ET=TETRAHEDRON");

//Output nodal coordinates

  fprintf(OUTPUT,"\n");
  itC = Coord.begin();
  itCEnd = Coord.end();
  for (; itC != itCEnd; ++itC) {
    fprintf(OUTPUT,"%22.15e %22.15e %22.15e",itC->second[0],itC->second[1],itC->second[2]);
    fprintf(OUTPUT,"\n");
  }

//Output tetrahedron connectivity
   std::map<int,std::vector<int>,std::less<int> >::const_iterator itS = SubElements.begin();
   std::map<int,std::vector<int>,std::less<int> >::const_iterator itSEnd = SubElements.end();

  for (; itS != itSEnd; ++itS) {
    fprintf(OUTPUT,"%d %d %d %d",itS->second[0],itS->second[1],itS->second[2],itS->second[3]);
    fprintf(OUTPUT,"\n");
  }

  fclose(OUTPUT);

//output the lines of each face of the (IxElement_c) element 

  FILE *OUTPUT2;
  OUTPUT2 = fopen("faces.dat","w");

  itF = Face.begin();
  itFEnd = Face.end();
  for (; itF != itFEnd; ++itF) {
    facenodes = itF->second;
    fprintf(OUTPUT2, "\nGEOMETRY T=LINE3D, C=GREEN, L =SOLID, LT = %12.5e, CS = GRID",0.5);
    fprintf(OUTPUT2,"\n%d",1); fprintf(OUTPUT2,"\n%d",facenodes.size()+1);
    for (j = 0; j < nbNodesFace; ++j) {
      itC = Coord.find(facenodes[j]);
      xyz = itC->second;
      fprintf(OUTPUT2,"\n%22.15e %22.15e %22.15e ",xyz[0],xyz[1],xyz[2]);
    }
//repeat the first node of the face
    facenodes = itF->second;
    itC = Coord.find(facenodes[0]); xyz = itC->second;
    fprintf(OUTPUT2,"\n%22.15e %22.15e %22.15e ",xyz[0],xyz[1],xyz[2]);
  }
  fclose(OUTPUT2);

  return 1;
}