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

/*******************************************************************************************
*
*                              T W O - 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   S E G M E N T
*      -------------------------------------------------------------------------------
*
*                         N. Sukumar, University of California, Davis
*
*                                   C + +   C O D E
*                                   ---------------
*
* Purpose
* =======
* Partition a 2-dimensional finite element (triangle or quadrilateral) by a segment that
* is given by a point on it and its bi-normal (m) [t x m is a vector in the z-direction,
* where t is the tangent vector of the segment].  The nodal connectivities of the finite
* element are stored in counter-clockwise orientation (ccw). The finite element is 
* partitioned into 3-noded triangular elements on either side of the segment. 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
* 2-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 -DIxTRI3  main.C [triangular element], or
* 2. g++ -Wno-deprecated -DIxQUAD4 main.C [quadrilateral element], 
* with a.out 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 segment (such as a crack segment)
* is that of the class GxElement_c in the X-FEM code.  The number of nodes, faces (edges
* in R^2), maximum node and element IDs, connectivity information for the element and the 
* edges are set for the triangle (IxTRI3) as well as the quadrilateral (IxQUAD4) element 
* to test this code
*
* Output
* ======
* The graphics plotting package TECPLOT [Amtec Inc.] is used to display the results. The
* edges of the finite element are output in Tecplot format to the file "edges.dat",
* the partitioned elements are output in Tecplot format to the file "subelements2d.dat", 
* and the crack segment is output in Tecplot format to the file "crack.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 segment (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 edges 
* of the finite element is carried out, and if there is an intersection for an edge, the 
* point of intersection is computed and added to the node and face 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 edge lists of the 
* finite element. Now, using the edge lists of the finite element, two new edge lists for the 
* surfaces ABOVE and BELOW the (crack) segment are created.  The two edges (upper and lower) 
* that correspond to the (crack) segement are also created such that ccw orientation of the 
* nodes is preserved. The centroid of the two domains is evaluated with orientation = 1000 
* for the ABOVE surface and orientation = -1000 for BELOW surface. A loop over the two edge
* lists is performed and using the centroid for each surface, 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]); 
//
// tests if a point is ON, ABOVE, or BELOW the given segment
// int OnAboveOrBelow(double px, double py, double pz, double p[3], double m[3]);
//
// determines the intersection of the given segment 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;
}

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
// parametric representation of a line: p1 + t*(p2 - p1), where t = [0,1] if 
// there is an intersection
  
  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;
}


int main() {

using namespace std; // using std namespace

// Partitioning a 2D element that is cut by a (crack) segment into sub-elements (triangles) 
// for the purpose of numerical integration 

  int nbNodes;                // no of nodes for the (IxElement_c) element

  int nbFaces;                // no of faces (edges) of the IxElement_c

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

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

  int maxNodeId;              // max node Id

  int maxElemId;              // max element Id

#ifdef IxTRI3 // 3-noded constant strain finite element
  nbNodes = 3; 
  nbFaces = 3; 
  nbNodesFace = 2;
  maxElemId = 1; 
  maxFaceId = 3;
  maxNodeId = 3;
#elif defined(IxQUAD4) // 4-noded bilinear quadrilateral element
  nbNodes = 4; 
  nbFaces = 4; 
  nbNodesFace = 2; 
  maxElemId = 1; 
  maxFaceId = 4;
  maxNodeId = 4;
#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, typedef and const iterators for node ID to nodal coordinates
  std::map<int,std::vector<double>,std::less<int> > Coord;
  std::map<int,std::vector<double>,std::less<int> >::const_iterator itC, itCEnd;
  typedef std::map<int,std::vector<double>,std::less<int> >::value_type coord_add;
// 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;

// INPUT PARAMETERS 
// segment is given by its bi-normal m and a POINT p

  double m[3] = {1.0, 0.3, 0.0};
  double p[3] = {0.6, 0.4, 0.0};

// Initializations

  int i, j, k;
  double xyz[3];
  std::vector<double> xyzvec;

// node Id to nodal coordinates map

#ifdef IxTRI3 // 3-noded constant strain finite element
  cout << "\n Partitioning a IxTRI3 which is cut by a crack segment\n";
  double x[3] = {0.,1.,0.};
  double y[3] = {0.,0.,1.};
  double z[3] = {0.,0.,0.};
  int elemconnect[3] = {1,2,3};

  int v[3][2] = { {0,1}, {1,2}, {2,0} };

#elif defined(IxQUAD4) // 4-noded bilinear quadrilateral element
  cout << "\n Partitioning a IxQUAD4 which is cut by a crack segment\n";
  double x[4] = {0.,1.,1.,0.};
  double y[4] = {0.,0.,1.,1.};
  double z[4] = {0.,0.,0.,0.};
  int elemconnect[4] = {1,2,3,4};

  int v[4][2] = { {0,1}, {1,2}, {2,3}, {3,0} };

#endif

// node Id to nodal coordinates map
  for (i = 0; i < nbNodes; ++i) {
    xyz[0] = x[i]; 
    xyz[1] = y[i]; 
    xyz[2] = z[i]; 
    xyzvec.clear();
    for (int ii = 0; ii < 3; ++ii) xyzvec.push_back(xyz[ii]);
    Coord.insert(coord_add(i+1,xyzvec)); // start from ID 1
  }
 
// element connectivity
  for (i = 0; i < nbNodes; ++i) {
    connect.push_back(elemconnect[i]);
  }

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

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

  double px, py, pz;
  int orient;

  for (i = 0; i < nbNodes; ++i) {
    itC = Coord.find(connect[i]);
    xyzvec = itC->second;
    px = xyzvec[0]; py = xyzvec[1]; pz = xyzvec[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 line does not cut the element or is 
// along an element edge. 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 edge is cut by the (crack) segment; 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 edges to check edges that are cut by the (crack) segment

  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) {
    isCut = false;
    orient_compare = 0;
    facenodes = itF->second;
// two nodes per face/edge in 2D
    itA1 = Abo.find(facenodes[0]);
    orient1 = itA1->second;
    itA2 = Abo.find(facenodes[1]);
    orient2 = itA2->second;
    if (orient1*orient2 == -1) isCut = true;
    cout << "\n Face " << itF->first << " is " << (isCut ? "" : "NOT ") << "CUT by the crack ";   

// for edges cut by the crack, find the intersection point

    if (isCut) {

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

      itC = Coord.find(facenodes[0]);
      xyzvec = itC->second;         
      p1x = xyzvec[0]; p1y = xyzvec[1]; p1z = xyzvec[2]; 
      itC = Coord.find(facenodes[1]);
      xyzvec = itC->second;         
      p2x = xyzvec[0]; p2y = xyzvec[1]; p2z = xyzvec[2];

// 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;

      xyzvec.clear();
      xyzvec.push_back(px); xyzvec.push_back(py); xyzvec.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 segment, 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,xyzvec));

// insert previous node and the point (ID is maxNodeId) in the NEW connectivity list                      
      NodesPoints.push_back(facenodes[0]);
      NodesPoints.push_back(maxNodeId);
      NodesPoints.push_back(facenodes[1]);
    }
    else {
// copy the edge to the new edge (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 connectivity for the nodes/points above and below the
// crack segment

  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 edge is below the (crack) segment
        below.push_back(itN->second[i]);
      } 
      if (itO->second == 1 || itO->second == 0 || itO->second == -10 || itO->second == 10) {
//   node/point of the edge is above the (crack) segment
        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) segment for the top and bottom edges

  std::vector<int> unsorted;
  std::vector<double> vec1, vec2;

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

  for (; itO != itOEnd; ++itO) {
    if (itO->second == 0 || itO->second == -10 || itO->second == 10) {
      unsorted.push_back(itO->first);
    } 

  }

  if (unsorted.size() != 2) assert(0); // *ONLY* two points on the (crack) segement

  itC = Coord.find(unsorted[0]); vec1 = itC->second;
  itC = Coord.find(unsorted[1]); vec2 = itC->second; 
  xyz[0] = - (vec2[1] - vec1[1]);
  xyz[1] = vec2[0] - vec1[0];
  xyz[2] = 0.0;

// sort the nodes on the OnAbove surface of the crack so that the connectivity is 
// correct -> a ccw connectivity results in a +ve sense; for the OnBelow surface too 
// ccw connectivity results in a positive sense

  if (vecdotvec(xyz,m) > 0.0) { // keep the same order of the nodes
    above.push_back(unsorted[0]);
    above.push_back(unsorted[1]);
  }
  else { // switch the order of the nodes
    above.push_back(unsorted[1]);
    above.push_back(unsorted[0]);
  }

// ccw connectivity for the OnBelow surface
  below.push_back(above[1]);
  below.push_back(above[0]);

// create the edge OnAbove using the sorted points from vector above

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

// create the edge OnBelow using the sorted points from vector below

// maxFaceId is updated in the case for OnAbove
  FaceBelowCrack.insert(new_add(maxFaceId,below)); below.clear();

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

  double below_centroid[3] = {0.,0.,0.};
  double above_centroid[3] = {0.,0.,0.};

  itO = Abo.begin();
  itOEnd = Abo.end();
  int nbabove = 0, nbbelow = 0;
  for (; itO != itOEnd; ++itO) {
    itC = Coord.find(itO->first);     
    xyzvec = itC->second;
    if (itO->second == -1 || itO->second == 0 || itO->second == -10 || itO->second == 10) {
      for (i = 0, ++nbbelow; i < 3; ++i) below_centroid[i] += xyzvec[i]; 
    }
    if (itO->second == 1 || itO->second == 0 || itO->second == 10 || itO->second == -10) {
      for (i = 0, ++nbabove; i < 3; ++i) above_centroid[i] += xyzvec[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 area above and below the crack segment,
// respectively)

  orient = 1000; maxNodeId++; 
  int aboveAreaPoint = maxNodeId;
  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));
 
  orient = -1000; maxNodeId++; 
  int belowAreaPoint = 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));

// construct triangles (sub-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[3];
  std::vector<int> connvec;

  itN = FaceAboveCrack.begin();
  itNEnd = FaceAboveCrack.end();
  for (; itN != itNEnd; ++itN) {
    subconn[0] = itN->second[0];
    subconn[1] = itN->second[1];
    subconn[2] = aboveAreaPoint;
    ++maxElemId;
    connvec.clear();
    for (int ii = 0; ii < 3; ++ii) connvec.push_back(subconn[ii]);
    SubElements.insert(sub_add(maxElemId,connvec));
  }

  itN = FaceBelowCrack.begin();
  itNEnd = FaceBelowCrack.end();
  for (; itN != itNEnd; ++itN) {
    subconn[0] = itN->second[0];
    subconn[1] = itN->second[1];
    subconn[2] = belowAreaPoint;
    ++maxElemId;
    connvec.clear();
    for (int ii = 0; ii < 3; ++ii) connvec.push_back(subconn[ii]);
    SubElements.insert(sub_add(maxElemId,connvec));
  }

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

  itN = NewFaceConn.begin();
  itNEnd = NewFaceConn.end();
  for (; itN != itNEnd; ++itN) {
    cout << "\n Edge 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 Edge 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("subelements2d.dat","w");

  fprintf(OUTPUT, " TITLE = \"Example: Triangular Sub-Element Data\" "); fprintf(OUTPUT,"\n");
  fprintf(OUTPUT, " VARIABLES = \"X\" ");
  fprintf(OUTPUT, ", \"Y\" "); fprintf(OUTPUT,"\n");

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

// Output nodal coordinates

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

// Output triangle 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",itS->second[0],itS->second[1],itS->second[2]);
    fprintf(OUTPUT,"\n");
  }

  fclose(OUTPUT);

// output each edge of the (IxElement_c) element in Tecplot format

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

   itF = Face.begin();
   itFEnd = Face.end();
   for (; itF != itFEnd; ++itF) {
     facenodes = itF->second;
     fprintf(OUTPUT2, "GEOMETRY T=LINE, C=GREEN, L =SOLID, LT = %12.5e, CS = GRID\n",0.5);
     fprintf(OUTPUT2,"%d\n",1); fprintf(OUTPUT2,"%d\n",nbNodesFace);
     for (j = 0; j < nbNodesFace; ++j) {
       itC = Coord.find(facenodes[j]);
       xyzvec = itC->second;
       fprintf(OUTPUT2,"%22.15e %22.15e\n",xyzvec[0],xyzvec[1]);
     }
   }
   fclose(OUTPUT2);

// output the crack segment

   FILE *OUTPUT3;
   OUTPUT3 = fopen("crack.dat","w");

// two arbitrary points on the crack: p1[3] and p2[3]
   double p1[3] = { p[0] + m[1], p[1] - m[0], 0.0 };
   double p2[3] = { p[0] - m[1], p[1] + m[0], 0.0 };
   fprintf(OUTPUT3, "GEOMETRY T=LINE, C=BLACK, L =SOLID, LT = %12.5e,CS = GRID\n",0.5);
   fprintf(OUTPUT3,"%d\n%d\n",1,2); 
   fprintf(OUTPUT3,"%22.15e %22.15e\n",p1[0],p1[1]);
   fprintf(OUTPUT3,"%22.15e %22.15e",p2[0],p2[1]);
   fclose(OUTPUT3);

   return 1;
}