mesh.cpp

Go to the documentation of this file.

#include "mesh.h"
#define DEBUG_REFINE
 
//****************************************************
#pragma mark -
#pragma mark * local ( static ) function prototypes *
//----------------------------------------------------
static double CurrentTime(void)
{
  //     static double scale = 0.0;
        
  //     if (0.0 == scale) {
  //         mach_timebase_info_data_t info;
  //         mach_timebase_info(&info);
  //         scale = info.numer / info.denom * 1e-9;
  //     }
        
  //     return mach_absolute_time() * scale;
  return 0.0;
}       // CurrentTime
//****************************************************
#pragma mark -
#pragma mark * exported function implementations *
//----------------------------------------------------
 
#ifndef REAL
#define REAL double
#define REAL_LOCAL
#endif
#include PROTEUS_TRIANGLE_H
#ifdef REAL_LOCAL
#undef REAL
#endif
 
#include "meshio.h"
#include <set>
#include <valarray>
#include <string>
#include <sstream>
#include <iostream>
#include <fstream>
#include <cstring>
 
//mwftodo decide where to put mesh type tags
const int DEFAULT_ELEMENT_MATERIAL=0;
const int DEFAULT_NODE_MATERIAL=-1;
const int INTERIOR_NODE_MATERIAL=0;
const int EXTERIOR_NODE_MATERIAL=1;
const int INTERIOR_ELEMENT_BOUNDARY_MATERIAL=0;
const int EXTERIOR_ELEMENT_BOUNDARY_MATERIAL=1;
 
//todo compute geometric info, node star
extern "C"
{
 
  int edgeMeshElements(const int& nx, Mesh& mesh)
  {
    //todo put in check to see if memory was allocated, do something about material types
    mesh.nNodes_element = 2;
    mesh.nNodes_global = nx;
    mesh.nElements_global = nx-1;
    mesh.elementNodesArray = new int[mesh.nElements_global*mesh.nNodes_element];
    mesh.elementMaterialTypes = new int[mesh.nElements_global];
    memset(mesh.elementMaterialTypes,DEFAULT_ELEMENT_MATERIAL,mesh.nElements_global*sizeof(int));
    for(int i=0,eN=0;i<nx-1;i++)
      eN=newEdge(eN,mesh.elementNodesArray,i,i+1);
    return 0;
  }
  
  int regularEdgeMeshNodes(const int& nx, const double& Lx, Mesh& mesh)
  {
    const double hx=Lx/(nx-1.0);
    mesh.nodeArray = new double[mesh.nNodes_global*3];
    memset(mesh.nodeArray,0,nx*3*sizeof(double));
    mesh.nodeMaterialTypes = new int[mesh.nNodes_global];
    //set interior and exterior node material flags after get boundary info in
    //constructElementBoundaryElementsArray_edge 
    //if nodeMaterialTypes left as DEFAULT_NODE_MATERIAL
    memset(mesh.nodeMaterialTypes,DEFAULT_NODE_MATERIAL,mesh.nNodes_global*sizeof(int));
    for(int i=0;i<nx;i++)
      mesh.nodeArray[i*3+0]=i*hx;
    return 0;
  }
  
  int regularRectangularToTriangularMeshElements(const int& nx, const int& ny, Mesh& mesh,int triangleFlag)
  {
    mesh.nNodes_element = 3;
    mesh.nElements_global = 2*(nx-1)*(ny-1);
    mesh.nNodes_global = nx*ny;
    mesh.elementNodesArray = new int[mesh.nElements_global*mesh.nNodes_element];
    mesh.elementMaterialTypes = new int[mesh.nElements_global];
    memset(mesh.elementMaterialTypes,DEFAULT_ELEMENT_MATERIAL,mesh.nElements_global*sizeof(int));
    //for local refinement
    mesh.newestNodeBases = new int[mesh.nElements_global];
 
    for(int i=0,eN=0;i<ny-1;i++)
      for(int j=0;j<nx-1;j++)
        {
          int 
            n0 =(j+0) + (i+0)*nx,
            n1 =(j+1) + (i+0)*nx,
            n2 =(j+0) + (i+1)*nx,
            n3 =(j+1) + (i+1)*nx;
          if (triangleFlag == 2)//left leaning
            {
              eN=newTriangle(eN,mesh.elementNodesArray,n0,n1,n2);
              eN=newTriangle(eN,mesh.elementNodesArray,n2,n1,n3);
 
            }
          else if (triangleFlag == 1)//union jack
            {
              if (i%2 + j%2 == 0 || i%2 + j%2 == 2)
                {
                  eN=newTriangle(eN,mesh.elementNodesArray,n0,n1,n3);
                  eN=newTriangle(eN,mesh.elementNodesArray,n0,n3,n2);
                }
              else
                {
                  eN=newTriangle(eN,mesh.elementNodesArray,n0,n1,n2);
                  eN=newTriangle(eN,mesh.elementNodesArray,n2,n1,n3);
                }
            }
          else //right leaning
            {
              eN=newTriangle(eN,mesh.elementNodesArray,n0,n1,n3);
              eN=newTriangle(eN,mesh.elementNodesArray,n0,n3,n2);
            }
          //NW element is (n0,n3,n1), SE element is (n0,n2,n3)
          //eN was incremented twice just now
          mesh.newestNodeBases[eN-2] = 2; //SE across from node n1
          mesh.newestNodeBases[eN-1] = 1; //NW is across from node n2
        }
    return 0;
  }
  
  int regularRectangularToTriangularElementBoundaryMaterials(const double& Lx, const double& Ly, Mesh& mesh)
  {
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN,nN_0,nN_1;
        double x_0,y_0,x_1,y_1,epsilon=1.0e-8;
        ebN = mesh.exteriorElementBoundariesArray[ebNE];
        nN_0 = mesh.elementBoundaryNodesArray[ebN*2 + 0];
        nN_1 = mesh.elementBoundaryNodesArray[ebN*2 + 1];
        x_0 = mesh.nodeArray[nN_0*3+0];
        y_0 = mesh.nodeArray[nN_0*3+1];
        x_1 = mesh.nodeArray[nN_1*3+0];
        y_1 = mesh.nodeArray[nN_1*3+1];
        if (y_0 <= epsilon && y_1 <= epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 1;
        else if (y_0 >= Ly - epsilon && y_1 >= Ly -  epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 3;
        else if (x_0 <= epsilon && x_1 <= epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 4;
        else if (x_0 >= Lx - epsilon && x_1 >= Lx - epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 2;
        else
          assert(false);
      }
    return 0;
  }
 
  int regularRectangularToTriangularMeshNodes(const int& nx, const int& ny, const double& Lx, const double& Ly, Mesh& mesh)
  {
    const double hx=Lx/(nx-1.0),hy=Ly/(ny-1.0);
    mesh.nodeArray = new double[mesh.nNodes_global*3];
    memset(mesh.nodeArray,0,mesh.nNodes_global*3*sizeof(double));
    mesh.nodeMaterialTypes = new int[mesh.nNodes_global];
    //set interior and exterior node material flags after get boundary info in
    //constructElementBoundaryElementsArray_edge 
    //if nodeMaterialTypes left as DEFAULT_NODE_MATERIAL
    memset(mesh.nodeMaterialTypes,DEFAULT_NODE_MATERIAL,mesh.nNodes_global*sizeof(int));
    int nN;
    for(int i=0;i<ny;i++)
      for(int j=0;j<nx;j++)
        {
          nN = i*nx+j;
          mesh.nodeArray[3*nN+0]=j*hx;
          mesh.nodeArray[3*nN+1]=i*hy;
          if (i==0)
            mesh.nodeMaterialTypes[nN] = 1;
          else if(i==ny-1)
            mesh.nodeMaterialTypes[nN] = 3;
          else if (j==0)
            mesh.nodeMaterialTypes[nN] = 4;
          else if(j==nx-1)
            mesh.nodeMaterialTypes[nN] = 2;
          else
            mesh.nodeMaterialTypes[nN] = 0;
        }
    return 0;
  }
  
  inline void reorientNodes_tet(double* nodeArray, int* nodes)
  {
    int &n0(nodes[0]),&n1(nodes[1]),&n2(nodes[2]),&n3(nodes[3]);
    double t[3][3];
    t[0][0] = nodeArray[n1*3+0] - nodeArray[n0*3+0];
    t[0][1] = nodeArray[n1*3+1] - nodeArray[n0*3+1];
    t[0][2] = nodeArray[n1*3+2] - nodeArray[n0*3+2];
    
    t[1][0] = nodeArray[n2*3+0] - nodeArray[n0*3+0];
    t[1][1] = nodeArray[n2*3+1] - nodeArray[n0*3+1];
    t[1][2] = nodeArray[n2*3+2] - nodeArray[n0*3+2];
    
    t[2][0] = nodeArray[n3*3+0] - nodeArray[n0*3+0];
    t[2][1] = nodeArray[n3*3+1] - nodeArray[n0*3+1];
    t[2][2] = nodeArray[n3*3+2] - nodeArray[n0*3+2];
    
    register double det = t[0][0]*(t[1][1]*t[2][2] - t[1][2]*t[2][1]) -   
      t[0][1]*(t[1][0]*t[2][2] - t[1][2]*t[2][0]) +       
      t[0][2]*(t[1][0]*t[2][1] - t[1][1]*t[2][0]);
    
    if (det < 0)
      {
        int tmp = n2;
        n2 = n3;
        n3 = tmp;
      }
  }
  
  int reorientTetrahedralMesh(Mesh& mesh)
  {
    for (int eN=0;eN<mesh.nElements_global;eN++)
      reorientNodes_tet(mesh.nodeArray,&mesh.elementNodesArray[eN*mesh.nNodes_element]);
    return 0;
  }
  
  int regularHexahedralToTetrahedralMeshElements(const int& nx, 
                                                 const int& ny, 
                                                 const int& nz, 
                                                 Mesh& mesh)
  {
    mesh.nNodes_element = 4;
    mesh.nNodes_global=nx*ny*nz;
    mesh.nElements_global = 6*(nx-1)*(ny-1)*(nz-1);
    int nxy=nx*ny;
    //std::cout<<"element nodes array"<<std::endl;
    mesh.elementNodesArray = new int[mesh.nElements_global*mesh.nNodes_element];
    mesh.elementMaterialTypes = new int[mesh.nElements_global];
    memset(mesh.elementMaterialTypes,DEFAULT_ELEMENT_MATERIAL,mesh.nElements_global*sizeof(int));
    for(int i=0,eN=0;i<nz-1;i++)
      for(int j=0;j<ny-1;j++)
        for(int k=0;k<nx-1;k++)
          {
            int 
              n0 = (k+0) + (j+0)*nx + (i+0)*nxy,
              n1 = (k+1) + (j+0)*nx + (i+0)*nxy,
              n2 = (k+0) + (j+1)*nx + (i+0)*nxy,
              n3 = (k+1) + (j+1)*nx + (i+0)*nxy,
              n4 = (k+0) + (j+0)*nx + (i+1)*nxy,
              n5 = (k+1) + (j+0)*nx + (i+1)*nxy,
              n6 = (k+0) + (j+1)*nx + (i+1)*nxy,
              n7 = (k+1) + (j+1)*nx + (i+1)*nxy;
            eN=newTetrahedron(eN,mesh.elementNodesArray,n0,n1,n3,n7);
            eN=newTetrahedron(eN,mesh.elementNodesArray,n0,n1,n5,n7);
            eN=newTetrahedron(eN,mesh.elementNodesArray,n0,n2,n3,n7);
            eN=newTetrahedron(eN,mesh.elementNodesArray,n0,n2,n6,n7);
            eN=newTetrahedron(eN,mesh.elementNodesArray,n0,n4,n5,n7);
            eN=newTetrahedron(eN,mesh.elementNodesArray,n0,n4,n6,n7);
          }
    return 0;
  }
  
  int regularHexahedralToTetrahedralElementBoundaryMaterials(const double& Lx, const double& Ly, const double& Lz, Mesh& mesh)
  {
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN,nN_0,nN_1,nN_2;
        double x_0,y_0,z_0,
          x_1,y_1,z_1,
          x_2,y_2,z_2,
          epsilon=1.0e-8;
        ebN = mesh.exteriorElementBoundariesArray[ebNE];
        nN_0 = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary + 0];
        nN_1 = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary + 1];
        nN_2 = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary + 2];
 
        x_0 = mesh.nodeArray[nN_0*3+0];
        y_0 = mesh.nodeArray[nN_0*3+1];
        z_0 = mesh.nodeArray[nN_0*3+2];
 
        x_1 = mesh.nodeArray[nN_1*3+0];
        y_1 = mesh.nodeArray[nN_1*3+1]; 
        z_1 = mesh.nodeArray[nN_1*3+2];
 
        x_2 = mesh.nodeArray[nN_2*3+0];
        y_2 = mesh.nodeArray[nN_2*3+1]; 
        z_2 = mesh.nodeArray[nN_2*3+2];
 
        if (z_0 <= epsilon && z_1 <= epsilon && z_2 <= epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 1;
        else if (z_0 >= Lz - epsilon && z_1 >= Lz -  epsilon && z_2 >= Lz -  epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 4;
        else if (y_0 <= epsilon && y_1 <= epsilon && y_2 <= epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 2;
        else if (y_0 >= Ly - epsilon && y_1 >= Ly -  epsilon && y_2 >= Ly -  epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 6;
        else if (x_0 <= epsilon && x_1 <= epsilon && x_2 <= epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 3;
        else if (x_0 >= Lx - epsilon && x_1 >= Lx - epsilon && x_2 >= Lx - epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 5;
        else
          assert(false);
      }
    return 0;
  }
  
  int regularHexahedralMeshElementBoundaryMaterials(const double& Lx, const double& Ly, const double& Lz, Mesh& mesh)
  {
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN,nN_0,nN_1,nN_2, nN_3;
        double x_0,y_0,z_0,
          x_1,y_1,z_1,
          x_2,y_2,z_2,
          x_3,y_3,z_3,
          epsilon=1.0e-8;
        ebN = mesh.exteriorElementBoundariesArray[ebNE];
        nN_0 = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary + 0];
        nN_1 = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary + 1];
        nN_2 = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary + 2];
        nN_3 = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary + 3];
 
        x_0 = mesh.nodeArray[nN_0*3+0];
        y_0 = mesh.nodeArray[nN_0*3+1];
        z_0 = mesh.nodeArray[nN_0*3+2];
 
        x_1 = mesh.nodeArray[nN_1*3+0];
        y_1 = mesh.nodeArray[nN_1*3+1]; 
        z_1 = mesh.nodeArray[nN_1*3+2];
 
        x_2 = mesh.nodeArray[nN_2*3+0];
        y_2 = mesh.nodeArray[nN_2*3+1]; 
        z_2 = mesh.nodeArray[nN_2*3+2];
        
        x_3 = mesh.nodeArray[nN_3*3+0];
        y_3 = mesh.nodeArray[nN_3*3+1]; 
        z_3 = mesh.nodeArray[nN_3*3+2];
 
        if (z_0 <= epsilon && z_1 <= epsilon && z_2 <= epsilon && z_3 <= epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 1;
        else if (z_0 >= Lz - epsilon && z_1 >= Lz -  epsilon && z_2 >= Lz - epsilon && z_3 >= Lz -  epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 4;
        else if (y_0 <= epsilon && y_1 <= epsilon && y_2 <= epsilon && y_3 <= epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 2;
        else if (y_0 >= Ly - epsilon && y_1 >= Ly -  epsilon && y_2 >= Ly -  epsilon && y_3 >= Ly -  epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 6;
        else if (x_0 <= epsilon && x_1 <= epsilon && x_2 <= epsilon && x_3 <= epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 3;
        else if (x_0 >= Lx - epsilon && x_1 >= Lx - epsilon && x_2 >= Lx - epsilon && x_3 >= Lx - epsilon)
          mesh.elementBoundaryMaterialTypes[ebN] = 5;
        else
          assert(false);
      }
    return 0;
  }
 
  int regularQuadrilateralMeshElementBoundaryMaterials(const double& Lx, const double& Ly, Mesh& mesh)
  {
    int result = regularRectangularToTriangularElementBoundaryMaterials(Lx, Ly, mesh);
    return result;
  }
  
  int regularMeshNodes(const int& nx, 
                       const int& ny, 
                       const int& nz,
                       const double& Lx, 
                       const double& Ly, 
                       const double& Lz,
                       Mesh& mesh)
  {
    const int nxy=nx*ny;
    const double hx=Lx/(nx-1.0),
      hy=Ly/(ny-1.0),
      hz=Lz/(nz-1.0);
    mesh.nNodes_global=nx*ny*nz;
    //std::cout<<"node array length "<<mesh.nNodes_global*3<<std::endl;
    mesh.nodeArray = new double[mesh.nNodes_global*3];
    //std::cout<<"memset "<<std::endl;
    memset(mesh.nodeArray,0,mesh.nNodes_global*3*sizeof(double));
    //std::cout<<"mat types"<<std::endl;
    mesh.nodeMaterialTypes = new int[mesh.nNodes_global];
    //set interior and exterior node material flags after get boundary info in
    //constructElementBoundaryElementsArray_edge 
    //if nodeMaterialTypes left as DEFAULT_NODE_MATERIAL
    //std::cout<<"memset"<<std::endl;
    memset(mesh.nodeMaterialTypes,DEFAULT_NODE_MATERIAL,mesh.nNodes_global*sizeof(int));
    int nN;
    for(int i=0;i<nz;i++)
      for(int j=0;j<ny;j++)
        for(int k=0;k<nx;k++)
          {
            nN = k + j*nx + i*nxy;
            mesh.nodeArray[3*nN+0]=k*hx;
            mesh.nodeArray[3*nN+1]=j*hy;
            mesh.nodeArray[3*nN+2]=i*hz;
            if (i==0)
              mesh.nodeMaterialTypes[nN] = 1;
            else if(i==nz-1)
              mesh.nodeMaterialTypes[nN] = 4;
            else if (j==0)
              mesh.nodeMaterialTypes[nN] = 2;
            else if(j==ny-1)
              mesh.nodeMaterialTypes[nN] = 6;
            else if (k==0)
              mesh.nodeMaterialTypes[nN] = 3;
            else if(k==nx-1)
              mesh.nodeMaterialTypes[nN] = 5;
            else
              mesh.nodeMaterialTypes[nN] = 0;
          }
    return 0;
  }
 
  int regularMeshNodes2D(const int& nx, 
                         const int& ny,
                         const double& Lx, 
                         const double& Ly, 
                         Mesh& mesh)
  {
    const double hx=Lx/(nx-1.0),
      hy=Ly/(ny-1.0);
    mesh.nNodes_global=nx*ny;   
    mesh.nodeArray = new double[mesh.nNodes_global*3];
    memset(mesh.nodeArray,0,mesh.nNodes_global*3*sizeof(double));
    mesh.nodeMaterialTypes = new int[mesh.nNodes_global];
    //set interior and exterior node material flags after get boundary info in
    //constructElementBoundaryElementsArray_edge 
    //if nodeMaterialTypes left as DEFAULT_NODE_MATERIAL
    memset(mesh.nodeMaterialTypes,DEFAULT_NODE_MATERIAL,mesh.nNodes_global*sizeof(int));
    int nN;
    for(int j=0;j<ny;j++)
      for(int k=0;k<nx;k++)
        {
          nN = k*ny + j;
          mesh.nodeArray[3*nN+0]=k*hx;
          mesh.nodeArray[3*nN+1]=j*hy;
          mesh.nodeArray[3*nN+2]=0.0;
          if (j==0)
            mesh.nodeMaterialTypes[nN] = 1;
          else if(j==ny-1)
            mesh.nodeMaterialTypes[nN] = 3;
          else if (k==0)
            mesh.nodeMaterialTypes[nN] = 4;
          else if(k==nx-1)
            mesh.nodeMaterialTypes[nN] = 2;
          else
            mesh.nodeMaterialTypes[nN] = 0;
        }
    return 0;
  }
 
  int regularHexahedralToTetrahedralMeshNodes(const int& nx, 
                                              const int& ny, 
                                              const int& nz,
                                              const double& Lx, 
                                              const double& Ly, 
                                              const double& Lz,
                                              Mesh& mesh)
  {
    regularMeshNodes(nx,ny,nz,Lx,Ly,Lz,mesh);      
    //std::cout<<"regularHexahedralToTetrahedralMeshNodes is Deprecated\n";    
    return 0;                                     
  }                                            
    
  int regularHexahedralMeshElements(const int& nx, 
                                    const int& ny, 
                                    const int& nz, 
                                    const int& px, 
                                    const int& py, 
                                    const int& pz, 
                                    Mesh& mesh)
  {
    std::cout<<"mesh elements"<<std::endl;
    mesh.nNodes_element = 8;
    mesh.nNodes_global=nx*ny*nz;
    mesh.nElements_global = (nx-1)*(ny-1)*(nz-1);
    int nxy=nx*ny;
 
    mesh.nx=nx;
    mesh.ny=ny;
    mesh.nz=nz;
 
    //mesh.px=px;
    //mesh.py=py;
    //mesh.pz=pz;
 
    mesh.elementNodesArray = new int[mesh.nElements_global*mesh.nNodes_element];
    mesh.elementMaterialTypes = new int[mesh.nElements_global];
    memset(mesh.elementMaterialTypes,DEFAULT_ELEMENT_MATERIAL,mesh.nElements_global*sizeof(int));
    int eN=0;
    for(int i=0;i<nz-1;i++)
      for(int j=0;j<ny-1;j++)
        for(int k=0;k<nx-1;k++)
          {
            eN=newHexahedron(eN,mesh.elementNodesArray,
                             (k+0) + (j+0)*nx + (i+0)*nxy,
                             (k+1) + (j+0)*nx + (i+0)*nxy,
                             (k+1) + (j+1)*nx + (i+0)*nxy,
                             (k+0) + (j+1)*nx + (i+0)*nxy,
                             (k+0) + (j+0)*nx + (i+1)*nxy,
                             (k+1) + (j+0)*nx + (i+1)*nxy,
                             (k+1) + (j+1)*nx + (i+1)*nxy,
                             (k+0) + (j+1)*nx + (i+1)*nxy);
          }
    return 0;
  }
  
  int regularQuadrilateralMeshElements(const int& nx, 
                                       const int& ny,
                                       Mesh& mesh)
  {
    mesh.nNodes_element = 4;
    mesh.nNodes_global=nx*ny;
    mesh.nElements_global = (nx-1)*(ny-1);
 
    mesh.nx=nx;
    mesh.ny=ny;
    mesh.nz=1;
 
    mesh.elementNodesArray = new int[mesh.nElements_global*mesh.nNodes_element];
    mesh.elementMaterialTypes = new int[mesh.nElements_global];
    memset(mesh.elementMaterialTypes,DEFAULT_ELEMENT_MATERIAL,mesh.nElements_global*sizeof(int));
    int eN=0;
    for(int j=0;j<ny-1;j++)
      for(int k=0;k<nx-1;k++)
        {
          eN=newQuadrilateral(eN,mesh.elementNodesArray,
                              (k+0)*ny + (j+0),
                              (k+0)*ny + (j+1),
                              (k+1)*ny + (j+1),
                              (k+1)*ny + (j+0));
        }
    return 0;
  }
  
  int regularNURBSMeshElements(const int& nx, 
                               const int& ny, 
                               const int& nz, 
                               const int& px, 
                               const int& py, 
                               const int& pz, 
                               Mesh& mesh)
  {
 
    mesh.nNodes_element = (px+1)*(py+1)*(pz+1);
    mesh.nNodes_global=nx*ny*nz;
    mesh.nElements_global = (nx-px)*(ny-py)*(nz-pz);
    
    mesh.nx=nx;
    mesh.ny=ny;
    mesh.nz=nz;
    mesh.px=px;
    mesh.py=py;
    mesh.pz=pz;
    
    int nxy=nx*ny;
    mesh.elementNodesArray = new int[mesh.nElements_global*mesh.nNodes_element];
    mesh.elementMaterialTypes = new int[mesh.nElements_global];
    //memset(mesh.elementMaterialTypes,DEFAULT_ELEMENT_MATERIAL,mesh.nElements_global*sizeof(int));
    mesh.elementIJK = new int[mesh.nElements_global*3];
 
    int eN=0;
    for(int i=0;i<nz-px;i++)
      for(int j=0;j<ny-py;j++)
        for(int k=0;k<nx-pz;k++)
          {
            mesh.elementIJK[eN*3+0] = i;
            mesh.elementIJK[eN*3+1] = j;
            mesh.elementIJK[eN*3+2] = k;
 
            int sN = 0;
            for(int ii=0;ii<px+1;ii++)
              for(int jj=0;jj<py+1;jj++)
                for(int kk=0;kk<pz+1;kk++)
                  {             
                    mesh.elementNodesArray[eN*mesh.nNodes_element+sN] = (k+kk) + (j+jj)*nx + (i+ii)*nxy;                    
                    sN++;
                  }    
            eN++;               
 
          } 
                      
    mesh.U_KNOT = new double[nx+px+1];  
    mesh.V_KNOT = new double[ny+py+1];    
    mesh.W_KNOT = new double[nz+pz+1];
 
    for(int i=0;i<px+1;i++)
      mesh.U_KNOT[i] = 0.0;
    for(int i=px+1;i<nx;i++)
      mesh.U_KNOT[i] = double(i-px-1);       
    for(int i=nx;i<nx+px+1;i++)
      mesh.U_KNOT[i] = double(nx);
 
    for(int i=0;i<py+1;i++)
      mesh.V_KNOT[i] = 0.0;
    for(int i=py+1;i<ny;i++)
      mesh.V_KNOT[i] = double(i-py-1);       
    for(int i=ny;i<ny+py+1;i++)
      mesh.V_KNOT[i] = double(ny);
       
    for(int i=0;i<pz+1;i++)
      mesh.W_KNOT[i] = 0.0;
    for(int i=pz+1;i<pz;i++)
      mesh.W_KNOT[i] = double(i-pz-1);       
    for(int i=nz;i<nz+pz+1;i++)
      mesh.W_KNOT[i] = double(nz);           
 
 
    mesh.weights = new double[mesh.nNodes_global];
 
    for(int i=0;i<mesh.nNodes_global;i++)
      mesh.weights[i] = 1.0;
    //std::cout<<"NURBS MESH BUILD"<<std::endl;             
    return 0;
  }
 
  int constructElementBoundaryElementsArray_edge(Mesh& mesh)
  {
    mesh.nNodes_elementBoundary = 1;
    mesh.nElementBoundaries_element = 2;
    using namespace std;
    double start,stop;
    //cout<<"Constructing element boundary map"<<endl;
    map<NodeTuple<1>,
      ElementNeighbors> elementBoundaryElements;
    start=CurrentTime();
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<2;ebN++)
        {
          register int nodes[1];
          nodes[0] = mesh.elementNodesArray[eN*2+(ebN+1)%2];
          NodeTuple<1> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
 
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    //mwf added
    mesh.elementBoundariesArray= new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Populating Arrays"<<endl;
    start = CurrentTime();
    int ebN=0;
    for(map<NodeTuple<1>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++,ebN++)
      {
        mesh.elementBoundaryNodesArray[ebN] = eb->first.nodes[0];
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left;
          }
        else
          exteriorElementBoundaries.insert(ebN);
        //mwf added
        mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = ebN;
        if (eb->second.right != -1)
          {
            mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = ebN;
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    //edges are elements in 1D so we just copy over but could be slicker
    mesh.nEdges_global = mesh.nElements_global;
    mesh.edgeNodesArray = new int[mesh.nElements_global*2];
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        mesh.edgeNodesArray[eN*2+0] = mesh.elementNodesArray[eN*2+0];
        mesh.edgeNodesArray[eN*2+1] = mesh.elementNodesArray[eN*2+1];
      }
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        nodeStar[mesh.edgeNodesArray[eN*2+0]].insert(mesh.edgeNodesArray[eN*2+1]);
        nodeStar[mesh.edgeNodesArray[eN*2+1]].insert(mesh.edgeNodesArray[eN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1]+nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    stop = CurrentTime();
    //repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
 
 
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
  int constructElementBoundaryElementsArray_triangle(Mesh& mesh)
  {
    using namespace std;
    mesh.nNodes_elementBoundary = 2;
    mesh.nElementBoundaries_element = 3;
    double start,stop;
    //cout<<"Constructing element boundary map"<<endl;
    map<NodeTuple<2>,
      ElementNeighbors> elementBoundaryElements;
    start=CurrentTime();
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<3;ebN++)
        {
          register int nodes[2];
          nodes[0] = mesh.elementNodesArray[eN*3+(ebN+1)%3];
          nodes[1] = mesh.elementNodesArray[eN*3+(ebN+2)%3];
          NodeTuple<2> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
    
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    mesh.elementBoundariesArray= new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Populating arrays"<<endl;
    start = CurrentTime();
    int ebN=0;
    for(map<NodeTuple<2>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++,ebN++)
      {
        mesh.elementBoundaryNodesArray[ebN*2 + 0] = eb->first.nodes[0];
        mesh.elementBoundaryNodesArray[ebN*2 + 1] = eb->first.nodes[1];
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left;
          }
        else
          exteriorElementBoundaries.insert(ebN);          
        mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = ebN;
        if (eb->second.right != -1)
          {
            mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = ebN;
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    set<NodeTuple<2> > edges;
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        int nodes[2];
        for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
          for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
            {
              nodes[0] = mesh.elementNodesArray[eN*3+nN_L];
              nodes[1] = mesh.elementNodesArray[eN*3+nN_R];
              edges.insert(NodeTuple<2>(nodes));
            }
      }
    mesh.nEdges_global = edges.size();
    mesh.edgeNodesArray = new int[mesh.nEdges_global*2];
    int edgeN=0;
    for (set<NodeTuple<2> >::iterator edgeTuple_p=edges.begin();edgeTuple_p != edges.end();edgeTuple_p++,edgeN++)
      {
        mesh.edgeNodesArray[edgeN*2+0] = edgeTuple_p->nodes[0];
        mesh.edgeNodesArray[edgeN*2+1] = edgeTuple_p->nodes[1];
      }
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1]+nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    stop = CurrentTime();
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    //repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    //set actual values elsewhere
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
  int constructElementBoundaryElementsArray_quadrilateral(Mesh& mesh)
  {
    using namespace std;
    mesh.nNodes_elementBoundary = 2;
    mesh.nElementBoundaries_element = 4;
    double start,stop;
    //cout<<"Constructing element boundary map"<<endl;
    map<NodeTuple<2>,
      ElementNeighbors> elementBoundaryElements;
    start=CurrentTime();
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<4;ebN++)
        {
          register int nodes[2];
          nodes[0] = mesh.elementNodesArray[eN*4+ebN];
          nodes[1] = mesh.elementNodesArray[eN*4+(ebN+1)%4];
          NodeTuple<2> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
    
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    mesh.elementBoundariesArray= new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Populating arrays"<<endl;
    start = CurrentTime();
    int ebN=0;
    for(map<NodeTuple<2>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++,ebN++)
      {
        mesh.elementBoundaryNodesArray[ebN*2 + 0] = eb->first.nodes[0];
        mesh.elementBoundaryNodesArray[ebN*2 + 1] = eb->first.nodes[1];
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left;
          }
        else
          exteriorElementBoundaries.insert(ebN);          
        mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = ebN;
        if (eb->second.right != -1)
          {
            mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = ebN;
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    set<NodeTuple<2> > edges;
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        int nodes[2];
        for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
          for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
            {
              nodes[0] = mesh.elementNodesArray[eN*4+nN_L];
              nodes[1] = mesh.elementNodesArray[eN*4+nN_R];
              edges.insert(NodeTuple<2>(nodes));
            }
      }
    mesh.nEdges_global = edges.size();
    mesh.edgeNodesArray = new int[mesh.nEdges_global*2];
    int edgeN=0;
    for (set<NodeTuple<2> >::iterator edgeTuple_p=edges.begin();edgeTuple_p != edges.end();edgeTuple_p++,edgeN++)
      {
        mesh.edgeNodesArray[edgeN*2+0] = edgeTuple_p->nodes[0];
        mesh.edgeNodesArray[edgeN*2+1] = edgeTuple_p->nodes[1];
      }
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1]+nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    stop = CurrentTime();
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    //repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    //set actual values elsewhere
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
  int constructElementBoundaryElementsArray_tetrahedron(Mesh& mesh)
  {
    mesh.nNodes_elementBoundary = 3;
    mesh.nElementBoundaries_element = 4;
    using namespace std;
    double start,stop;
    typedef map<NodeTuple<3>, ElementNeighbors> ElementBoundaryElementsType;
    {
    ElementBoundaryElementsType elementBoundaryElements;
    start=CurrentTime();
    //std::cout<<"Extracting boundary elements"<<endl;
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<mesh.nElementBoundaries_element;ebN++)
        {
          register int nodes[3];
          nodes[0] = mesh.elementNodesArray[eN*4+((ebN+1)%4)];
          nodes[1] = mesh.elementNodesArray[eN*4+((ebN+2)%4)];
          nodes[2] = mesh.elementNodesArray[eN*4+((ebN+3)%4)];
          NodeTuple<3> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
            }
        }
    stop = CurrentTime();
    //std::cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //std::cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
 
    //std::cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    typedef set<int> ElementBoundariesType;
    ElementBoundariesType interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    mesh.elementBoundariesArray= new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //std::cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
    //std::cout<<"Generating elementBoundaryElementsArray and elementBoundaryNodesArray"<<endl;
    start = CurrentTime();
    int ebN=0;
    for(map<NodeTuple<3>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++,ebN++)
      {
        mesh.elementBoundaryNodesArray[ebN*3 + 0] = eb->first.nodes[0];
        mesh.elementBoundaryNodesArray[ebN*3 + 1] = eb->first.nodes[1];
        mesh.elementBoundaryNodesArray[ebN*3 + 2] = eb->first.nodes[2];
 
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left; 
          }
        else
          exteriorElementBoundaries.insert(ebN);
        mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = ebN;
        if (eb->second.right != -1)
          {
            mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = ebN;
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    }
    //elementBoundaryElements.~ElementBoundaryElementsType();
    //interiorElementBoundaries.~ElementBoundariesType();
    //exteriorElementBoundaries.~ElementBoundariesType();
    {
      typedef set<NodeTuple<2> > EdgesType;
    EdgesType edges;
    //std::cout<<"extracting edges"<<std::endl;
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        int nodes[2];
        for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
          for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
            {
              nodes[0] = mesh.elementNodesArray[eN*4+nN_L];
              nodes[1] = mesh.elementNodesArray[eN*4+nN_R];
              edges.insert(NodeTuple<2>(nodes));
            }
      }
    mesh.nEdges_global = edges.size();
    mesh.edgeNodesArray = new int[mesh.nEdges_global*2];
    EdgesType::iterator edge_p=edges.begin();
    for (int edgeN=0;edgeN<int(edges.size());edgeN++)
      {
        mesh.edgeNodesArray[edgeN*2+0] = edge_p->nodes[0];
        mesh.edgeNodesArray[edgeN*2+1] = edge_p->nodes[1];
        edge_p++;
      }
    }//edges.~EdgesType();
    {
      //std::cout<<"node star"<<std::endl;
    typedef vector<set<int> > NodeStarType;
    NodeStarType nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1] + nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    }//nodeStar.~NodeStarType();
    {
      //std::cout<<"node elements star"<<std::endl;
    stop = CurrentTime();
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    typedef vector<set<int> > NodeElementsStarType;
    NodeElementsStarType nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    }//nodeElementsStar.~NodeElementsStarType();
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //std::cout<<"nodematerials"<<std::endl;
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //std::cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
  int constructElementBoundaryElementsArray_hexahedron(Mesh& mesh)
  {
    mesh.nNodes_elementBoundary = 4;
    mesh.nElementBoundaries_element = 6;
    using namespace std;
    double start,stop;
    map<NodeTuple<4>,
      ElementNeighbors> elementBoundaryElements;
    start=CurrentTime();
    
    int lface[6][4] = {{0,1,2,3},
                       {0,1,5,4},
                       {1,2,6,5},
                       {2,3,7,6},
                       {3,0,4,7},
                       {4,5,6,7}};
    
    //cout<<"Extracting boundary elements"<<endl;
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<mesh.nElementBoundaries_element;ebN++)
        {
          register int nodes[4];
          nodes[0] = mesh.elementNodesArray[eN*8+lface[ebN][0]];
          nodes[1] = mesh.elementNodesArray[eN*8+lface[ebN][1]];
          nodes[2] = mesh.elementNodesArray[eN*8+lface[ebN][2]];      
          nodes[3] = mesh.elementNodesArray[eN*8+lface[ebN][3]];
                    
          NodeTuple<4> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
            }
        }        
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
 
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    //mwf added
    mesh.elementBoundariesArray= new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Generating elementBoundaryElementsArray and elementBoundaryNodesArray"<<endl;
    start = CurrentTime();
    int ebN=0;
    for(map<NodeTuple<4>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++,ebN++)
      {
        mesh.elementBoundaryNodesArray[ebN*4 + 0] = eb->first.nodes_unsorted[0];
        mesh.elementBoundaryNodesArray[ebN*4 + 1] = eb->first.nodes_unsorted[1];
        mesh.elementBoundaryNodesArray[ebN*4 + 2] = eb->first.nodes_unsorted[2];
        mesh.elementBoundaryNodesArray[ebN*4 + 3] = eb->first.nodes_unsorted[3];
 
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left; 
          }
        else
          exteriorElementBoundaries.insert(ebN);    
          //mwf added
            mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = ebN;
            if (eb->second.right != -1)
            {
               mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = ebN;
            }
      }    
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray  =  new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    set<NodeTuple<2> > edges;
     
    int ledge[12][2] = {{0,1},{1,2},{2,3},{3,0},
                        {0,4},{1,5},{2,6},{3,7},
                        {4,5},{5,6},{6,7},{7,4}};
    
    
    
    
    //std::cout<<"Extracting edges"<<std::endl;
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        int nodes[2];
        for (int e=0;e<12;e++)
          {
        
            nodes[0] = mesh.elementNodesArray[eN*8+ledge[e][0]];
            nodes[1] = mesh.elementNodesArray[eN*8+ledge[e][1]];
                  
            edges.insert(NodeTuple<2>(nodes));
          }   
      }
      
      
    mesh.nEdges_global = edges.size();
    mesh.edgeNodesArray = new int[mesh.nEdges_global*2];
    set<NodeTuple<2> >::iterator edge_p=edges.begin();
    for (int edgeN=0;edgeN<int(edges.size());edgeN++)
      {
        mesh.edgeNodesArray[edgeN*2+0] = edge_p->nodes[0];
        mesh.edgeNodesArray[edgeN*2+1] = edge_p->nodes[1];
        edge_p++;
      }
 
    //std::cout<<"Extracting nodeStar"<<std::endl;      
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
 
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      {
        mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1] + nodeStar[nN-1].size();
      }  
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    stop = CurrentTime();
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    //mwf repeat for node-->elements arrays
 
    //std::cout<<"Extracting nodeElementsStar"<<std::endl;   
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {         
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
 
    //std::cout<<"Set material types"<<std::endl;     
    //mwf end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
        
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
 
  int constructElementBoundaryElementsArray_NURBS(Mesh& mesh)
  {
    using namespace std;
 
    int n0 = 0;
    int n1 = mesh.px ;
    int n2 = (mesh.px+1)*(mesh.py+1)-1 ;
    int n3 = (mesh.px+1)*mesh.py;
 
 
    int n4 = (mesh.px+1)*(mesh.py+1)*mesh.pz + n0;
    int n5 = (mesh.px+1)*(mesh.py+1)*mesh.pz + n1;
    int n6 = (mesh.px+1)*(mesh.py+1)*mesh.pz + n2;
    int n7 = (mesh.px+1)*(mesh.py+1)*mesh.pz + n3;
    
    
    
    mesh.nNodes_elementBoundary = 4;
    mesh.nElementBoundaries_element = 6;
    using namespace std;
    double start,stop;
    map<NodeTuple<4>,
      ElementNeighbors> elementBoundaryElements;
    start=CurrentTime();
    
    int lface[6][4] = {{n0,n1,n2,n3},
                       {n0,n1,n5,n4},
                       {n1,n2,n6,n5},
                       {n2,n3,n7,n6},
                       {n3,n0,n4,n7},
                       {n4,n5,n6,n7}};
    
    
    
    //cout<<"Extracting boundary elements"<<endl;
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<mesh.nElementBoundaries_element;ebN++)
        {
          register int nodes[4];
          nodes[0] = mesh.elementNodesArray[eN*8+lface[ebN][0]];
          nodes[1] = mesh.elementNodesArray[eN*8+lface[ebN][1]];
          nodes[2] = mesh.elementNodesArray[eN*8+lface[ebN][2]];      
          nodes[3] = mesh.elementNodesArray[eN*8+lface[ebN][3]];
                    
          NodeTuple<4> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
 
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    //mwf added
    mesh.elementBoundariesArray= new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Generating elementBoundaryElementsArray and elementBoundaryNodesArray"<<endl;
    start = CurrentTime();
    int ebN=0;
    for(map<NodeTuple<4>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++,ebN++)
      {
        mesh.elementBoundaryNodesArray[ebN*4 + 0] = eb->first.nodes[0];
        mesh.elementBoundaryNodesArray[ebN*4 + 1] = eb->first.nodes[1];
        mesh.elementBoundaryNodesArray[ebN*4 + 2] = eb->first.nodes[2];
        mesh.elementBoundaryNodesArray[ebN*4 + 3] = eb->first.nodes[3];
 
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left; 
          }
        else
          exteriorElementBoundaries.insert(ebN);          
        //mwf added
        mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = ebN;
        if (eb->second.right != -1)
          {
            mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = ebN;
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    set<NodeTuple<2> > edges;    
    
    
 
   
    int ledge[12][2] = {{n0,n1},{n1,n2},{n2,n3},{n3,n0},
                        {n0,n4},{n1,n5},{n2,n6},{n3,n7},
                        {n4,n5},{n5,n6},{n6,n7},{n7,n4}};
    
    
    
    
    //std::cout<<"extracting edges"<<std::endl;
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        int nodes[2];
        for (int e=0;e<12;e++)
          {
                
            nodes[0] = mesh.elementNodesArray[eN*8+ledge[e][0]];
            nodes[1] = mesh.elementNodesArray[eN*8+ledge[e][1]];
                     
            edges.insert(NodeTuple<2>(nodes));
          }   
      }
      
      
    mesh.nEdges_global = edges.size();
    mesh.edgeNodesArray = new int[mesh.nEdges_global*2];
    set<NodeTuple<2> >::iterator edge_p=edges.begin();
    for (int edgeN=0;edgeN<int(edges.size());edgeN++)
      {
        mesh.edgeNodesArray[edgeN*2+0] = edge_p->nodes[0];
        mesh.edgeNodesArray[edgeN*2+1] = edge_p->nodes[1];
        edge_p++;
      }
      
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
 
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1] + nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    stop = CurrentTime();
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    //mwf repeat for node-->elements arrays
  
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {         
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    
    //mwf end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
        
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
  
  int constructElementBoundaryElementsArrayWithGivenElementBoundaryAndEdgeNumbers_NURBS(Mesh& mesh)
  {
 
    int n0 = 0;
    int n1 = mesh.px ;
    int n2 = (mesh.px+1)*(mesh.py+1)-1 ;
    int n3 = (mesh.px+1)*mesh.py;
 
 
    int n4 = (mesh.px+1)*(mesh.py+1)*mesh.pz + n0;
    int n5 = (mesh.px+1)*(mesh.py+1)*mesh.pz + n1;
    int n6 = (mesh.px+1)*(mesh.py+1)*mesh.pz + n2;
    int n7 = (mesh.px+1)*(mesh.py+1)*mesh.pz + n3;
 
    mesh.nNodes_elementBoundary = 4;
    mesh.nElementBoundaries_element = 6;
    assert(mesh.elementBoundariesArray);
    using namespace std;
    double start,stop;
    map<NodeTuple<4>,
      ElementNeighbors> elementBoundaryElements;
    map<NodeTuple<4>,
      int> elementBoundaryIds;
    start=CurrentTime();
    
    int lface[6][4] = {{n0,n1,n2,n3},
                       {n0,n1,n5,n4},
                       {n1,n2,n6,n5},
                       {n2,n3,n7,n6},
                       {n3,n0,n4,n7},
                       {n4,n5,n6,n7}};
    
    //cout<<"Extracting boundary elements"<<endl;
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<mesh.nElementBoundaries_element;ebN++)
        {
          register int ebN_global = mesh.elementBoundariesArray[eN*mesh.nElementBoundaries_element+ebN];
          register int nodes[4];
          nodes[0] = mesh.elementNodesArray[eN*8+lface[ebN][0]];
          nodes[1] = mesh.elementNodesArray[eN*8+lface[ebN][1]];
          nodes[2] = mesh.elementNodesArray[eN*8+lface[ebN][2]];      
          nodes[3] = mesh.elementNodesArray[eN*8+lface[ebN][3]];
          NodeTuple<4> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
 
              // assert(elementBoundaryIds[ebt] == ebN_global);
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
              elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
            }
        }
     stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
 
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Generating elementBoundaryElementsArray and elementBoundaryNodesArray"<<endl;
    start = CurrentTime();
    for(map<NodeTuple<4>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++)
      {
        int ebN = elementBoundaryIds[eb->first];
        mesh.elementBoundaryNodesArray[ebN*4 + 0] = eb->first.nodes[0];
        mesh.elementBoundaryNodesArray[ebN*4 + 1] = eb->first.nodes[1];
        mesh.elementBoundaryNodesArray[ebN*4 + 2] = eb->first.nodes[2];
        mesh.elementBoundaryNodesArray[ebN*4 + 3] = eb->first.nodes[3];
 
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left; 
          }
        else
          exteriorElementBoundaries.insert(ebN);          
        //assert(mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] == ebN);
        if (eb->second.right != -1)
          {
            //  assert(mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] == ebN);
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    //cek/ido todo figure out how to do this or if it's necessary
    //    assert(mesh.edgeNodesArray);
 
    //     set<NodeTuple<2> > edges;
    //     //std::cout<<"extracting edges"<<std::endl;
    //     for (int eN=0;eN<mesh.nElements_global;eN++)
    //       {
    //         int nodes[2];
    //         for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
    //           for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
    //             {
    //               nodes[0] = mesh.elementNodesArray[eN*4+nN_L];
    //               nodes[1] = mesh.elementNodesArray[eN*4+nN_R];
    //               edges.insert(NodeTuple<2>(nodes));
    //             }
    //       }
    //     mesh.nEdges_global = edges.size();
    //     mesh.edgeNodesArray = new int[mesh.nEdges_global*2];
    //     set<NodeTuple<2> >::iterator edge_p=edges.begin();
    //     for (int edgeN=0;edgeN<int(edges.size());edgeN++)
    //       {
    //         mesh.edgeNodesArray[edgeN*2+0] = edge_p->nodes[0];
    //         mesh.edgeNodesArray[edgeN*2+1] = edge_p->nodes[1];
    //         edge_p++;
    //       }
 
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1] + nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    stop = CurrentTime();
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    //mwf repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //mwf end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
 
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
  int constructElementBoundaryElementsArrayWithGivenElementBoundaryNumbers_edge(Mesh& mesh)
  {
    mesh.nNodes_elementBoundary = 1;
    mesh.nElementBoundaries_element = 2;
    assert(mesh.elementBoundariesArray);
    using namespace std;
    double start,stop;
    //cout<<"Constructing element boundary map"<<endl;
    map<NodeTuple<1>,
      ElementNeighbors> elementBoundaryElements;
    map<NodeTuple<1>,
      int> elementBoundaryIds;
    start=CurrentTime();
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<2;ebN++)
        {
          register int ebN_global = mesh.elementBoundariesArray[eN*2+ebN];
          register int nodes[1];
          nodes[0] = mesh.elementNodesArray[eN*2+(ebN+1)%2];
          NodeTuple<1> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
              assert(elementBoundaryIds[ebt] == ebN_global);
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
              elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
 
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Populating Arrays"<<endl;
    start = CurrentTime();
    for(map<NodeTuple<1>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++)
      {
        int ebN = elementBoundaryIds[eb->first];
        mesh.elementBoundaryNodesArray[ebN] = eb->first.nodes[0];
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left;
          }
        else
          exteriorElementBoundaries.insert(ebN);
        assert(mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] == ebN);
        if (eb->second.right != -1)
          {
            assert(mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] == ebN);
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    //edges are elements in 1D so we just copy over but could be slicker
    mesh.nEdges_global = mesh.nElements_global;
    mesh.edgeNodesArray = new int[mesh.nElements_global*2];
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        mesh.edgeNodesArray[eN*2+0] = mesh.elementNodesArray[eN*2+0];
        mesh.edgeNodesArray[eN*2+1] = mesh.elementNodesArray[eN*2+1];
      }
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        nodeStar[mesh.edgeNodesArray[eN*2+0]].insert(mesh.edgeNodesArray[eN*2+1]);
        nodeStar[mesh.edgeNodesArray[eN*2+1]].insert(mesh.edgeNodesArray[eN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1]+nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    stop = CurrentTime();
    //repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
 
 
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
  int constructElementBoundaryElementsArrayWithGivenElementBoundaryNumbers_triangle(Mesh& mesh)
  {
    using namespace std;
    mesh.nNodes_elementBoundary = 2;
    mesh.nElementBoundaries_element = 3;
    assert(mesh.elementBoundariesArray);
    double start,stop;
    //cout<<"Constructing element boundary map"<<endl;
    map<NodeTuple<2>,
      ElementNeighbors> elementBoundaryElements;
    map<NodeTuple<2>,
      int> elementBoundaryIds;
    start=CurrentTime();
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<3;ebN++)
        {
          register int ebN_global = mesh.elementBoundariesArray[eN*3+ebN];
          register int nodes[2];
          nodes[0] = mesh.elementNodesArray[eN*3+(ebN+1)%3];
          nodes[1] = mesh.elementNodesArray[eN*3+(ebN+2)%3];
          NodeTuple<2> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
              assert(elementBoundaryIds[ebt] == ebN_global);
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
              elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
    
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Populating arrays"<<endl;
    start = CurrentTime();
    for(map<NodeTuple<2>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++)
      {
        int ebN = elementBoundaryIds[eb->first];
        mesh.elementBoundaryNodesArray[ebN*2 + 0] = eb->first.nodes[0];
        mesh.elementBoundaryNodesArray[ebN*2 + 1] = eb->first.nodes[1];
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left;
          }
        else
          exteriorElementBoundaries.insert(ebN);          
        assert(mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] == ebN);
        if (eb->second.right != -1)
          {
            assert(mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] == ebN);
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    set<NodeTuple<2> > edges;
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        int nodes[2];
        for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
          for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
            {
              nodes[0] = mesh.elementNodesArray[eN*3+nN_L];
              nodes[1] = mesh.elementNodesArray[eN*3+nN_R];
              edges.insert(NodeTuple<2>(nodes));
            }
      }
    mesh.nEdges_global = edges.size();
    mesh.edgeNodesArray = new int[mesh.nEdges_global*2];
    int edgeN=0;
    for (set<NodeTuple<2> >::iterator edgeTuple_p=edges.begin();edgeTuple_p != edges.end();edgeTuple_p++,edgeN++)
      {
        mesh.edgeNodesArray[edgeN*2+0] = edgeTuple_p->nodes[0];
        mesh.edgeNodesArray[edgeN*2+1] = edgeTuple_p->nodes[1];
      }
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1]+nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    stop = CurrentTime();
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    //repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    //set actual values elsewhere
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
  
  int constructElementBoundaryElementsArrayWithGivenElementBoundaryNumbers_quadrilateral(Mesh& mesh)
  {
    using namespace std;
    mesh.nNodes_elementBoundary = 2;
    mesh.nElementBoundaries_element = 4;
    assert(mesh.elementBoundariesArray);
    double start,stop;
    //cout<<"Constructing element boundary map"<<endl;
    map<NodeTuple<2>,
      ElementNeighbors> elementBoundaryElements;
    map<NodeTuple<2>,
      int> elementBoundaryIds;
    start=CurrentTime();
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<4;ebN++)
        {
          register int ebN_global = mesh.elementBoundariesArray[eN*4+ebN];
          register int nodes[2];
          nodes[0] = mesh.elementNodesArray[eN*4+ebN];
          nodes[1] = mesh.elementNodesArray[eN*4+(ebN+1)%4];
          NodeTuple<2> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
              assert(elementBoundaryIds[ebt] == ebN_global);
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
              elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
    
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Populating arrays"<<endl;
    start = CurrentTime();
    for(map<NodeTuple<2>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++)
      {
        int ebN = elementBoundaryIds[eb->first];
        mesh.elementBoundaryNodesArray[ebN*2 + 0] = eb->first.nodes[0];
        mesh.elementBoundaryNodesArray[ebN*2 + 1] = eb->first.nodes[1];
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left;
          }
        else
          exteriorElementBoundaries.insert(ebN);          
        assert(mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] == ebN);
        if (eb->second.right != -1)
          {
            assert(mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] == ebN);
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    set<NodeTuple<2> > edges;
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        int nodes[2];
        for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
          {
            nodes[0] = mesh.elementNodesArray[eN*4+nN_L];
            nodes[1] = mesh.elementNodesArray[eN*4+(nN_L+1)%4];
            edges.insert(NodeTuple<2>(nodes));
          }
      }
    mesh.nEdges_global = edges.size();
    mesh.edgeNodesArray = new int[mesh.nEdges_global*2];
    int edgeN=0;
    for (set<NodeTuple<2> >::iterator edgeTuple_p=edges.begin();edgeTuple_p != edges.end();edgeTuple_p++,edgeN++)
      {
        mesh.edgeNodesArray[edgeN*2+0] = edgeTuple_p->nodes[0];
        mesh.edgeNodesArray[edgeN*2+1] = edgeTuple_p->nodes[1];
      }
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1]+nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    stop = CurrentTime();
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    //repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    //set actual values elsewhere
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
  int constructElementBoundaryElementsArrayWithGivenElementBoundaryNumbers_tetrahedron(Mesh& mesh)
  {
    mesh.nNodes_elementBoundary = 3;
    mesh.nElementBoundaries_element = 4;
    assert(mesh.elementBoundariesArray);
    using namespace std;
    double start,stop;
    map<NodeTuple<3>,
      ElementNeighbors> elementBoundaryElements;
    map<NodeTuple<3>,
      int> elementBoundaryIds;
    start=CurrentTime();
    //cout<<"Extracting boundary elements"<<endl;
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<mesh.nElementBoundaries_element;ebN++)
        {
          register int ebN_global = mesh.elementBoundariesArray[eN*mesh.nElementBoundaries_element+ebN];
          register int nodes[3];
          nodes[0] = mesh.elementNodesArray[eN*4+((ebN+1)%4)];
          nodes[1] = mesh.elementNodesArray[eN*4+((ebN+2)%4)];
          nodes[2] = mesh.elementNodesArray[eN*4+((ebN+3)%4)];
          NodeTuple<3> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
              assert(elementBoundaryIds[ebt] == ebN_global);
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
              elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
 
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Generating elementBoundaryElementsArray and elementBoundaryNodesArray"<<endl;
    start = CurrentTime();
    for(map<NodeTuple<3>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++)
      {
        int ebN = elementBoundaryIds[eb->first];
        mesh.elementBoundaryNodesArray[ebN*3 + 0] = eb->first.nodes[0];
        mesh.elementBoundaryNodesArray[ebN*3 + 1] = eb->first.nodes[1];
        mesh.elementBoundaryNodesArray[ebN*3 + 2] = eb->first.nodes[2];
 
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left; 
          }
        else
          exteriorElementBoundaries.insert(ebN);          
        assert(mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] == ebN);
        if (eb->second.right != -1)
          {
            assert(mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] == ebN);
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    set<NodeTuple<2> > edges;
    //std::cout<<"extracting edges"<<std::endl;
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        int nodes[2];
        for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
          for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
            {
              nodes[0] = mesh.elementNodesArray[eN*4+nN_L];
              nodes[1] = mesh.elementNodesArray[eN*4+nN_R];
              edges.insert(NodeTuple<2>(nodes));
            }
      }
    mesh.nEdges_global = edges.size();
    mesh.edgeNodesArray = new int[mesh.nEdges_global*2];
    set<NodeTuple<2> >::iterator edge_p=edges.begin();
    for (int edgeN=0;edgeN<int(edges.size());edgeN++)
      {
        mesh.edgeNodesArray[edgeN*2+0] = edge_p->nodes[0];
        mesh.edgeNodesArray[edgeN*2+1] = edge_p->nodes[1];
        edge_p++;
      }
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1] + nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    stop = CurrentTime();
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    //mwf repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //mwf end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
 
 
  int constructElementBoundaryElementsArrayWithGivenElementBoundaryAndEdgeNumbers_edge(Mesh& mesh)
  {
    mesh.nNodes_elementBoundary = 1;
    mesh.nElementBoundaries_element = 2;
    assert(mesh.elementBoundariesArray);
    using namespace std;
    double start,stop;
    //cout<<"Constructing element boundary map"<<endl;
    map<NodeTuple<1>,
      ElementNeighbors> elementBoundaryElements;
    map<NodeTuple<1>,
      int> elementBoundaryIds;
    start=CurrentTime();
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<2;ebN++)
        {
          register int ebN_global = mesh.elementBoundariesArray[eN*2+ebN];
          register int nodes[1];
          nodes[0] = mesh.elementNodesArray[eN*2+(ebN+1)%2];
          NodeTuple<1> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
              assert(elementBoundaryIds[ebt] == ebN_global);
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
              elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
 
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Populating Arrays"<<endl;
    start = CurrentTime();
    for(map<NodeTuple<1>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++)
      {
        int ebN = elementBoundaryIds[eb->first];
        mesh.elementBoundaryNodesArray[ebN] = eb->first.nodes[0];
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left;
          }
        else
          exteriorElementBoundaries.insert(ebN);
        assert(mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] == ebN);
        if (eb->second.right != -1)
          {
            assert(mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] == ebN);
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
 
    assert(mesh.edgeNodesArray);
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        nodeStar[mesh.edgeNodesArray[eN*2+0]].insert(mesh.edgeNodesArray[eN*2+1]);
        nodeStar[mesh.edgeNodesArray[eN*2+1]].insert(mesh.edgeNodesArray[eN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1]+nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    stop = CurrentTime();
    //repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
 
 
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
  int constructElementBoundaryElementsArrayWithGivenElementBoundaryAndEdgeNumbers_triangle(Mesh& mesh)
  {
    using namespace std;
    mesh.nNodes_elementBoundary = 2;
    mesh.nElementBoundaries_element = 3;
    assert(mesh.elementBoundariesArray);
    double start,stop;
    //cout<<"Constructing element boundary map"<<endl;
    map<NodeTuple<2>,
      ElementNeighbors> elementBoundaryElements;
    map<NodeTuple<2>,
      int> elementBoundaryIds;
    start=CurrentTime();
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<3;ebN++)
        {
          register int ebN_global = mesh.elementBoundariesArray[eN*3+ebN];
          register int nodes[2];
          nodes[0] = mesh.elementNodesArray[eN*3+(ebN+1)%3];
          nodes[1] = mesh.elementNodesArray[eN*3+(ebN+2)%3];
          NodeTuple<2> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
              assert(elementBoundaryIds[ebt] == ebN_global);
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
              elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
    
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Populating arrays"<<endl;
    start = CurrentTime();
    for(map<NodeTuple<2>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++)
      {
        int ebN = elementBoundaryIds[eb->first];
        mesh.elementBoundaryNodesArray[ebN*2 + 0] = eb->first.nodes[0];
        mesh.elementBoundaryNodesArray[ebN*2 + 1] = eb->first.nodes[1];
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left;
          }
        else
          exteriorElementBoundaries.insert(ebN);          
        assert(mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] == ebN);
        if (eb->second.right != -1)
          {
            assert(mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] == ebN);
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    assert(mesh.edgeNodesArray);
    //    set<NodeTuple<2> > edges;
    //     for (int eN=0;eN<mesh.nElements_global;eN++)
    //       {
    //         int nodes[2];
    //         for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
    //           for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
    //             {
    //               nodes[0] = mesh.elementNodesArray[eN*3+nN_L];
    //               nodes[1] = mesh.elementNodesArray[eN*3+nN_R];
    //               edges.insert(NodeTuple<2>(nodes));
    //             }
    //       }
    //     mesh.nEdges_global = edges.size();
    //     mesh.edgeNodesArray = new int[mesh.nEdges_global*2];
    //     int edgeN=0;
    //     for (set<NodeTuple<2> >::iterator edgeTuple_p=edges.begin();edgeTuple_p != edges.end();edgeTuple_p++,edgeN++)
    //       {
    //         mesh.edgeNodesArray[edgeN*2+0] = edgeTuple_p->nodes[0];
    //         mesh.edgeNodesArray[edgeN*2+1] = edgeTuple_p->nodes[1];
    //       }
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1]+nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    stop = CurrentTime();
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    //repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    //set actual values elsewhere
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
  int constructElementBoundaryElementsArrayWithGivenElementBoundaryAndEdgeNumbers_quadrilateral(Mesh& mesh)
  {
    using namespace std;
    mesh.nNodes_elementBoundary = 2;
    mesh.nElementBoundaries_element = 4;
    assert(mesh.elementBoundariesArray);
    double start,stop;
    //cout<<"Constructing element boundary map"<<endl;
    map<NodeTuple<2>,
      ElementNeighbors> elementBoundaryElements;
    map<NodeTuple<2>,
      int> elementBoundaryIds;
    start=CurrentTime();
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<4;ebN++)
        {
          register int ebN_global = mesh.elementBoundariesArray[eN*4+ebN];
          register int nodes[2];
          nodes[0] = mesh.elementNodesArray[eN*4+ebN];
          nodes[1] = mesh.elementNodesArray[eN*4+(ebN+1)%4];
          NodeTuple<2> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
              assert(elementBoundaryIds[ebt] == ebN_global);
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
              elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
    
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Populating arrays"<<endl;
    start = CurrentTime();
    for(map<NodeTuple<2>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++)
      {
        int ebN = elementBoundaryIds[eb->first];
        mesh.elementBoundaryNodesArray[ebN*2 + 0] = eb->first.nodes[0];
        mesh.elementBoundaryNodesArray[ebN*2 + 1] = eb->first.nodes[1];
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left;
          }
        else
          exteriorElementBoundaries.insert(ebN);          
        assert(mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] == ebN);
        if (eb->second.right != -1)
          {
            assert(mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] == ebN);
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    assert(mesh.edgeNodesArray);
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1]+nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    stop = CurrentTime();
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    //repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    //set actual values elsewhere
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
  int constructElementBoundaryElementsArrayWithGivenElementBoundaryAndEdgeNumbers_tetrahedron(Mesh& mesh)
  {
    mesh.nNodes_elementBoundary = 3;
    mesh.nElementBoundaries_element = 4;
    assert(mesh.elementBoundariesArray);
    using namespace std;
    double start,stop;
    map<NodeTuple<3>,
      ElementNeighbors> elementBoundaryElements;
    map<NodeTuple<3>,
      int> elementBoundaryIds;
    start=CurrentTime();
    //cout<<"Extracting boundary elements"<<endl;
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<mesh.nElementBoundaries_element;ebN++)
        {
          register int ebN_global = mesh.elementBoundariesArray[eN*mesh.nElementBoundaries_element+ebN];
          register int nodes[3];
          nodes[0] = mesh.elementNodesArray[eN*4+((ebN+1)%4)];
          nodes[1] = mesh.elementNodesArray[eN*4+((ebN+2)%4)];
          nodes[2] = mesh.elementNodesArray[eN*4+((ebN+3)%4)];
 
          NodeTuple<3> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
              assert(elementBoundaryIds[ebt] == ebN_global);
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
              elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
 
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Generating elementBoundaryElementsArray and elementBoundaryNodesArray"<<endl;
    start = CurrentTime();
    for(map<NodeTuple<3>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++)
      {
        int ebN = elementBoundaryIds[eb->first];
        mesh.elementBoundaryNodesArray[ebN*3 + 0] = eb->first.nodes[0];
        mesh.elementBoundaryNodesArray[ebN*3 + 1] = eb->first.nodes[1];
        mesh.elementBoundaryNodesArray[ebN*3 + 2] = eb->first.nodes[2];
 
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left; 
          }
        else
          exteriorElementBoundaries.insert(ebN);          
        assert(mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] == ebN);
        if (eb->second.right != -1)
          {
            assert(mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] == ebN);
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    assert(mesh.edgeNodesArray);
    //     set<NodeTuple<2> > edges;
    //     //std::cout<<"extracting edges"<<std::endl;
    //     for (int eN=0;eN<mesh.nElements_global;eN++)
    //       {
    //         int nodes[2];
    //         for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
    //           for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
    //             {
    //               nodes[0] = mesh.elementNodesArray[eN*4+nN_L];
    //               nodes[1] = mesh.elementNodesArray[eN*4+nN_R];
    //               edges.insert(NodeTuple<2>(nodes));
    //             }
    //       }
    //     mesh.nEdges_global = edges.size();
    //     mesh.edgeNodesArray = new int[mesh.nEdges_global*2];
    //     set<NodeTuple<2> >::iterator edge_p=edges.begin();
    //     for (int edgeN=0;edgeN<int(edges.size());edgeN++)
    //       {
    //         mesh.edgeNodesArray[edgeN*2+0] = edge_p->nodes[0];
    //         mesh.edgeNodesArray[edgeN*2+1] = edge_p->nodes[1];
    //         edge_p++;
    //       }
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1] + nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    //mwf repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //mwf end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
  int constructElementBoundaryElementsArrayWithGivenElementBoundaryAndEdgeNumbers_hexahedron(Mesh& mesh)
  {
    mesh.nNodes_elementBoundary = 4;
    mesh.nElementBoundaries_element = 6;
    assert(mesh.elementBoundariesArray);
    using namespace std;
    double start,stop;
    map<NodeTuple<4>,
      ElementNeighbors> elementBoundaryElements;
    map<NodeTuple<4>,
      int> elementBoundaryIds;
    start=CurrentTime();
    
    int lface[6][4] = {{0,1,2,3},
                       {0,1,5,4},
                       {1,2,6,5},
                       {2,3,7,6},
                       {3,0,4,7},
                       {4,5,6,7}};    
    
    //cout<<"Extracting boundary elements"<<endl;
    for(int eN=0;eN<mesh.nElements_global;eN++)
      for(int ebN=0;ebN<mesh.nElementBoundaries_element;ebN++)
        {
          register int ebN_global = mesh.elementBoundariesArray[eN*mesh.nElementBoundaries_element+ebN];
          register int nodes[4];
          nodes[0] = mesh.elementNodesArray[eN*8+lface[ebN][0]];
          nodes[1] = mesh.elementNodesArray[eN*8+lface[ebN][1]];
          nodes[2] = mesh.elementNodesArray[eN*8+lface[ebN][2]];      
          nodes[3] = mesh.elementNodesArray[eN*8+lface[ebN][3]];
          NodeTuple<4> ebt(nodes);
          if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
            {
              elementBoundaryElements[ebt].right=eN;
              elementBoundaryElements[ebt].right_ebN_element=ebN;
            }
          else
            {
              elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,ElementNeighbors(eN,ebN)));
              elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
            }
        }
    stop = CurrentTime();
    //cout<<"Elapsed time for building element boundary elements map= "<<(stop-start)<<"s"<<endl;
    mesh.nElementBoundaries_global = elementBoundaryElements.size();
    //cout<<"nElementBoundaries_global = "<<mesh.nElementBoundaries_global<<endl;
 
    //cout<<"Allocating Arrays"<<endl;
    start = CurrentTime();
    set<int> interiorElementBoundaries,exteriorElementBoundaries;
    mesh.elementBoundaryNodesArray =  new int[mesh.nElementBoundaries_global*mesh.nNodes_elementBoundary];
    mesh.elementBoundaryElementsArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementBoundaryLocalElementBoundariesArray = new int[mesh.nElementBoundaries_global*2];
    mesh.elementNeighborsArray = new int[mesh.nElements_global*mesh.nElementBoundaries_element];
    stop = CurrentTime();
    //cout<<"Elapsed time for allocating arrays = "<<(stop-start)<<"s"<<endl;
 
    //cout<<"Generating elementBoundaryElementsArray and elementBoundaryNodesArray"<<endl;
    start = CurrentTime();
    for(map<NodeTuple<4>,ElementNeighbors>::iterator eb=elementBoundaryElements.begin();
        eb != elementBoundaryElements.end();
        eb++)
      {
        int ebN = elementBoundaryIds[eb->first];
        mesh.elementBoundaryNodesArray[ebN*4 + 0] = eb->first.nodes[0];
        mesh.elementBoundaryNodesArray[ebN*4 + 1] = eb->first.nodes[1];
        mesh.elementBoundaryNodesArray[ebN*4 + 2] = eb->first.nodes[2];
        mesh.elementBoundaryNodesArray[ebN*4 + 3] = eb->first.nodes[3];
 
        mesh.elementBoundaryElementsArray[ebN*2 + 0] = eb->second.left;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 0] = eb->second.left_ebN_element;
        mesh.elementBoundaryElementsArray[ebN*2 + 1] = eb->second.right;
        mesh.elementBoundaryLocalElementBoundariesArray[ebN*2 + 1] = eb->second.right_ebN_element;
        mesh.elementNeighborsArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] = eb->second.right; 
        if(eb->second.right != -1)
          {
            interiorElementBoundaries.insert(ebN);
            mesh.elementNeighborsArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] = eb->second.left; 
          }
        else
          exteriorElementBoundaries.insert(ebN);          
        assert(mesh.elementBoundariesArray[eb->second.left*mesh.nElementBoundaries_element + eb->second.left_ebN_element] == ebN);
        if (eb->second.right != -1)
          {
            assert(mesh.elementBoundariesArray[eb->second.right*mesh.nElementBoundaries_element + eb->second.right_ebN_element] == ebN);
          }
      }
    mesh.nInteriorElementBoundaries_global = interiorElementBoundaries.size();
    mesh.interiorElementBoundariesArray = new int[mesh.nInteriorElementBoundaries_global];
    mesh.nExteriorElementBoundaries_global = exteriorElementBoundaries.size();
    mesh.exteriorElementBoundariesArray = new int[mesh.nExteriorElementBoundaries_global];
    int ebNI=0,ebNE=0;
    for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
      mesh.interiorElementBoundariesArray[ebNI] = *ebN;
    for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
      mesh.exteriorElementBoundariesArray[ebNE] = *ebN;
    assert(mesh.edgeNodesArray);
    //     set<NodeTuple<2> > edges;
    //     //std::cout<<"extracting edges"<<std::endl;
    //     for (int eN=0;eN<mesh.nElements_global;eN++)
    //       {
    //         int nodes[2];
    //         for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
    //           for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
    //             {
    //               nodes[0] = mesh.elementNodesArray[eN*4+nN_L];
    //               nodes[1] = mesh.elementNodesArray[eN*4+nN_R];
    //               edges.insert(NodeTuple<2>(nodes));
    //             }
    //       }
    //     mesh.nEdges_global = edges.size();
    //     mesh.edgeNodesArray = new int[mesh.nEdges_global*2];
    //     set<NodeTuple<2> >::iterator edge_p=edges.begin();
    //     for (int edgeN=0;edgeN<int(edges.size());edgeN++)
    //       {
    //         mesh.edgeNodesArray[edgeN*2+0] = edge_p->nodes[0];
    //         mesh.edgeNodesArray[edgeN*2+1] = edge_p->nodes[1];
    //         edge_p++;
    //       }
 
    vector<set<int> > nodeStar(mesh.nNodes_global);
    for (int edgeN=0;edgeN<mesh.nEdges_global;edgeN++)
      {
        nodeStar[mesh.edgeNodesArray[edgeN*2+0]].insert(mesh.edgeNodesArray[edgeN*2+1]);
        nodeStar[mesh.edgeNodesArray[edgeN*2+1]].insert(mesh.edgeNodesArray[edgeN*2+0]);
      }
    mesh.nodeStarOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeStarOffsets[0] = 0;
    for (int nN=1;nN<mesh.nNodes_global+1;nN++)
      mesh.nodeStarOffsets[nN] = mesh.nodeStarOffsets[nN-1] + nodeStar[nN-1].size();
    mesh.nodeStarArray = new int[mesh.nodeStarOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0;nN<mesh.nNodes_global;nN++)
      for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
        mesh.nodeStarArray[offset] = *nN_star;
    stop = CurrentTime();
    mesh.max_nNodeNeighbors_node=0;
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      mesh.max_nNodeNeighbors_node=max(mesh.max_nNodeNeighbors_node,mesh.nodeStarOffsets[nN+1]-mesh.nodeStarOffsets[nN]);
    //mwf repeat for node-->elements arrays
    vector<set<int> > nodeElementsStar(mesh.nNodes_global);
    for (int eN = 0; eN < mesh.nElements_global; eN++)
      {
        for (int nN = 0; nN < mesh.nNodes_element; nN++)
          nodeElementsStar[mesh.elementNodesArray[eN*mesh.nNodes_element+nN]].insert(eN);
      }
    mesh.nodeElementOffsets = new int[mesh.nNodes_global+1];
    mesh.nodeElementOffsets[0] = 0;
    for (int nN = 0; nN < mesh.nNodes_global; nN++)
      mesh.nodeElementOffsets[nN+1] = mesh.nodeElementOffsets[nN]+nodeElementsStar[nN].size();
    mesh.nodeElementsArray  = new int[mesh.nodeElementOffsets[mesh.nNodes_global]];
    for (int nN=0,offset=0; nN < mesh.nNodes_global; nN++)
      {
        for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
             eN_star++,offset++)
          {
            mesh.nodeElementsArray[offset] = *eN_star;
          }
      }
    //mwf end node-->elements construction
    mesh.elementBoundaryMaterialTypes = new int[mesh.nElementBoundaries_global];
    //if nodeMaterial is DEFAULT, go ahead and set to interior or exterior
    //depending on which boundary node belongs to. 
    //If node on at least one exterior boundary then it's exterior
    for (int ebNE = 0; ebNE < mesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = mesh.exteriorElementBoundariesArray[ebNE];
        mesh.elementBoundaryMaterialTypes[ebN] = EXTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = EXTERIOR_NODE_MATERIAL;
          }
      }
    for (int ebNI = 0; ebNI < mesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = mesh.interiorElementBoundariesArray[ebNI];
        mesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        for (int nN_local = 0; nN_local < mesh.nNodes_elementBoundary; nN_local++)
          {
            int nN = mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_local];
            if (mesh.nodeMaterialTypes[nN] == DEFAULT_NODE_MATERIAL)
              mesh.nodeMaterialTypes[nN] = INTERIOR_NODE_MATERIAL;
          }
      }
    //cout<<"Elapsed time for populating arrays = "<<(stop-start)<<"s"<<endl;
    return 0;
  }
 
 
  inline double edgeLength(int nL,int nR, const double* nodeArray)
  {
    register double dx,dy,dz;
    dx = nodeArray[nL*3+0] - nodeArray[nR*3+0];
    dy = nodeArray[nL*3+1] - nodeArray[nR*3+1];
    dz = nodeArray[nL*3+2] - nodeArray[nR*3+2];
    return sqrt(dx*dx+dy*dy+dz*dz);
  }
 
  inline double triangleArea(int n0, int n1, int n2, const double* nodeArray)
  {
    register double va[3],vb[3];
    va[0] = nodeArray[n1*3+0] - nodeArray[n0*3+0];
    va[1] = nodeArray[n1*3+1] - nodeArray[n0*3+1];
    va[2] = nodeArray[n1*3+2] - nodeArray[n0*3+2];
    vb[0] = nodeArray[n2*3+0] - nodeArray[n0*3+0];
    vb[1] = nodeArray[n2*3+1] - nodeArray[n0*3+1];
    vb[2] = nodeArray[n2*3+2] - nodeArray[n0*3+2];
    return 0.5*fabs(va[1]*vb[2] - vb[1]*va[2] - (va[0]*vb[2] - vb[0]*va[2]) + (va[0]*vb[1] - va[1]*vb[0]));
  }
 
  inline double tetrahedronVolume(int n0, int n1, int n2, int n3, const double* nodeArray)
  {
    register double t[3][3];
    t[0][0] = nodeArray[n1*3+0] - nodeArray[n0*3+0];
    t[0][1] = nodeArray[n1*3+1] - nodeArray[n0*3+1];
    t[0][2] = nodeArray[n1*3+2] - nodeArray[n0*3+2];
 
    t[1][0] = nodeArray[n2*3+0] - nodeArray[n0*3+0];
    t[1][1] = nodeArray[n2*3+1] - nodeArray[n0*3+1];
    t[1][2] = nodeArray[n2*3+2] - nodeArray[n0*3+2];
 
    t[2][0] = nodeArray[n3*3+0] - nodeArray[n0*3+0];
    t[2][1] = nodeArray[n3*3+1] - nodeArray[n0*3+1];
    t[2][2] = nodeArray[n3*3+2] - nodeArray[n0*3+2];
    return fabs(t[0][0]*(t[1][1]*t[2][2] - t[1][2]*t[2][1]) - \
                t[0][1]*(t[1][0]*t[2][2] - t[1][2]*t[2][0]) + \
                t[0][2]*(t[1][0]*t[2][1] - t[1][1]*t[2][0]))/6.0;
  }
 
  int allocateGeometricInfo_tetrahedron(Mesh& mesh)
  {
    mesh.elementDiametersArray = new double[mesh.nElements_global];
    mesh.elementInnerDiametersArray = new double[mesh.nElements_global];
    mesh.elementBoundaryDiametersArray = new double[mesh.nElementBoundaries_global];
    mesh.elementBarycentersArray  = new double[mesh.nElements_global*3];
    mesh.elementBoundaryBarycentersArray = new double[mesh.nElementBoundaries_global*3];
    mesh.nodeDiametersArray = new double[mesh.nNodes_global];
    mesh.nodeSupportArray = new double[mesh.nNodes_global];
    return 0;
  }
 
  int computeGeometricInfo_tetrahedron(Mesh& mesh)
  {
    memset(mesh.elementDiametersArray,0,mesh.nElements_global*sizeof(double));
    memset(mesh.elementInnerDiametersArray,0,mesh.nElements_global*sizeof(double));
    memset(mesh.elementBoundaryDiametersArray,0,mesh.nElementBoundaries_global*sizeof(double));
    memset(mesh.elementBarycentersArray,0,mesh.nElements_global*3*sizeof(double));
    memset(mesh.elementBoundaryBarycentersArray,0,mesh.nElementBoundaries_global*3*sizeof(double));
    memset(mesh.nodeDiametersArray,0,mesh.nNodes_global*sizeof(double));
    memset(mesh.nodeSupportArray,0,mesh.nNodes_global*sizeof(double));
    mesh.hMin = edgeLength(mesh.elementNodesArray[0],
                           mesh.elementNodesArray[1],
                           mesh.nodeArray);
    const double nNperElemInv = 1.0/double(mesh.nNodes_element);
    const double nNperElemBInv= 1.0/double(mesh.nNodes_elementBoundary);
    for (int ebN=0;ebN<mesh.nElementBoundaries_global;ebN++)
      {
        mesh.elementBoundaryBarycentersArray[ebN*3 + 0] = 0.0;
        mesh.elementBoundaryBarycentersArray[ebN*3 + 1] = 0.0;
        mesh.elementBoundaryBarycentersArray[ebN*3 + 2] = 0.0;
 
        for (int nN_L=0;nN_L<mesh.nNodes_elementBoundary;nN_L++)
          {
            mesh.elementBoundaryBarycentersArray[ebN*3 + 0] += 
              nNperElemBInv*mesh.nodeArray[mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_L]*3 + 0];
            mesh.elementBoundaryBarycentersArray[ebN*3 + 1] += 
              nNperElemBInv*mesh.nodeArray[mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_L]*3 + 1];
            mesh.elementBoundaryBarycentersArray[ebN*3 + 2] += 
              nNperElemBInv*mesh.nodeArray[mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN_L]*3 + 2];
 
            for (int nN_R=nN_L+1;nN_R<mesh.nNodes_elementBoundary;nN_R++) {
              mesh.elementBoundaryDiametersArray[ebN] = fmax(mesh.elementBoundaryDiametersArray[ebN],
                                                             edgeLength(mesh.elementBoundaryNodesArray[ebN*3+nN_L],
                                                                        mesh.elementBoundaryNodesArray[ebN*3+nN_R],
                                                                        mesh.nodeArray));
            }
          }
      } 
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        double volume, surfaceArea=0.0, hMin=0.0,hMax=0.0;
        volume = tetrahedronVolume(mesh.elementNodesArray[eN*4+0],
                                   mesh.elementNodesArray[eN*4+1],
                                   mesh.elementNodesArray[eN*4+2],
                                   mesh.elementNodesArray[eN*4+3],
                                   mesh.nodeArray);
        //loop over faces to get surface error
        for (int ebN=0;ebN<4;ebN++)
          {
            surfaceArea += triangleArea(mesh.elementNodesArray[eN*4+(ebN+1)%4],
                                        mesh.elementNodesArray[eN*4+(ebN+2)%4],
                                        mesh.elementNodesArray[eN*4+(ebN+3)%4],
                                        mesh.nodeArray);
          }
        hMin = 6.0*volume/surfaceArea;
        //hMax
        for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
          for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
            hMax = fmax(hMax,
                        edgeLength(mesh.elementNodesArray[eN*4+nN_L],
                                   mesh.elementNodesArray[eN*4+nN_R],
                                   mesh.nodeArray));
        mesh.elementInnerDiametersArray[eN] = hMin;
        mesh.elementDiametersArray[eN] = hMax;
        mesh.volume += volume;
        mesh.sigmaMax = fmax(hMax/hMin,mesh.sigmaMax);
        mesh.h = fmax(hMax,mesh.h);
        mesh.hMin = fmin(hMin,mesh.hMin);
        
        mesh.elementBarycentersArray[eN*3 + 0] = 0.0; 
        mesh.elementBarycentersArray[eN*3 + 1] = 0.0; 
        mesh.elementBarycentersArray[eN*3 + 2] = 0.0; 
        for (int nN=0;nN<mesh.nNodes_element;nN++)
          {
            mesh.elementBarycentersArray[eN*3 + 0] += 
              nNperElemInv*mesh.nodeArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]*3 + 0];
            mesh.elementBarycentersArray[eN*3 + 1] += 
              nNperElemInv*mesh.nodeArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]*3 + 1];
            mesh.elementBarycentersArray[eN*3 + 2] += 
              nNperElemInv*mesh.nodeArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]*3 + 2];
            mesh.nodeDiametersArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]] += hMax*volume;
            mesh.nodeSupportArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]] += volume;
          }
      }
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      {
        mesh.nodeDiametersArray[nN] /= mesh.nodeSupportArray[nN];
      }
//    printf("volume = %12.5e \n",mesh.volume);
//    printf("h = %12.5e \n",mesh.h);
//    printf("hMin = %12.5e \n",mesh.hMin);
//    printf("sigmaMax = %12.5e \n",mesh.sigmaMax);
    return 0;
  }
  inline double hexahedronVolume(int n0, int n1, int n2, int n3, int n4, int n5, int n6, int n7, const double* nodeArray)
  {
    register double t[3];
    t[0] = nodeArray[n0*3+0] - nodeArray[n1*3+0];
    t[1] = nodeArray[n0*3+1] - nodeArray[n3*3+1];
    t[2] = nodeArray[n0*3+2] - nodeArray[n4*3+2];
    
    return fabs(t[0]*t[1]*t[2]);
  }
  int allocateGeometricInfo_hexahedron(Mesh& mesh)
  {
    mesh.elementDiametersArray = new double[mesh.nElements_global];
    mesh.elementInnerDiametersArray = new double[mesh.nElements_global];
    mesh.elementBoundaryDiametersArray = new double[mesh.nElementBoundaries_global];
    mesh.elementBarycentersArray  = new double[mesh.nElements_global*3];
    mesh.elementBoundaryBarycentersArray = new double[mesh.nElementBoundaries_global*3];
    mesh.nodeDiametersArray = new double[mesh.nNodes_global];
    mesh.nodeSupportArray = new double[mesh.nNodes_global];
    return 0;
  }
 
  int computeGeometricInfo_hexahedron(Mesh& mesh)
  {
    memset(mesh.elementDiametersArray,0,mesh.nElements_global*sizeof(double));
    memset(mesh.elementInnerDiametersArray,0,mesh.nElements_global*sizeof(double));
    memset(mesh.elementBoundaryDiametersArray,0,mesh.nElementBoundaries_global*sizeof(double));
    memset(mesh.nodeDiametersArray,0,mesh.nNodes_global*sizeof(double));
    memset(mesh.nodeSupportArray,0,mesh.nNodes_global*sizeof(double));
 
    mesh.hMin = edgeLength(mesh.elementNodesArray[0],
                           mesh.elementNodesArray[1],
                           mesh.nodeArray);
    //const double nNperElemInv = 1.0/double(mesh.nNodes_element);
    //const double nNperElemBInv= 1.0/double(mesh.nNodes_elementBoundary);
    for (int ebN=0;ebN<mesh.nElementBoundaries_global;ebN++)
      {
 
 
        for (int nN_L=0;nN_L<mesh.nNodes_elementBoundary;nN_L++)
          {
            for (int nN_R=nN_L+1;nN_R<mesh.nNodes_elementBoundary;nN_R++)
              mesh.elementBoundaryDiametersArray[ebN] = fmax(mesh.elementBoundaryDiametersArray[ebN],
                                                             edgeLength(mesh.elementBoundaryNodesArray[ebN*4+nN_L],
                                                                                mesh.elementBoundaryNodesArray[ebN*4+nN_R],
                                                                                mesh.nodeArray));
            }
                  
      } 
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        double volume, hMin=0.0,hMax=0.0;
        volume = hexahedronVolume(mesh.elementNodesArray[eN*8+0],
                                  mesh.elementNodesArray[eN*8+1],
                                  mesh.elementNodesArray[eN*8+2],
                                  mesh.elementNodesArray[eN*8+3],
                                  mesh.elementNodesArray[eN*8+4],
                                  mesh.elementNodesArray[eN*8+5],
                                  mesh.elementNodesArray[eN*8+6],
                                  mesh.elementNodesArray[eN*8+7],
                                  mesh.nodeArray);
        
        //hMax
        hMax = 0.0;
        hMin = 9e99;
        for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
          for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
          {
            hMax = fmax(hMax,
                        edgeLength(mesh.elementNodesArray[eN*8+nN_L],
                                   mesh.elementNodesArray[eN*8+nN_R],
                                   mesh.nodeArray));
            hMin = fmin(hMin,
                        edgeLength(mesh.elementNodesArray[eN*8+nN_L],
                                   mesh.elementNodesArray[eN*8+nN_R],
                                   mesh.nodeArray));
        }
        mesh.elementInnerDiametersArray[eN] = hMin;
        mesh.elementDiametersArray[eN] = hMax;
        mesh.volume += volume;
        mesh.sigmaMax = fmax(hMax/hMin,mesh.sigmaMax);
        mesh.h = fmax(hMax,mesh.h);
        mesh.hMin = fmin(hMin,mesh.hMin);
 
        for (int nN=0;nN<mesh.nNodes_element;nN++)
          {
            mesh.nodeDiametersArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]] += hMax*volume;
            mesh.nodeSupportArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]] += volume;
          }
      }
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      {
        mesh.nodeDiametersArray[nN] /= mesh.nodeSupportArray[nN];
      }
 
    //printf("volume = %12.5e \n",mesh.volume);
    //printf("h = %12.5e \n",mesh.h);
    //printf("hMin = %12.5e \n",mesh.hMin);
    //printf("sigmaMax = %12.5e \n",mesh.sigmaMax);
    return 0;
  }
 
  int allocateGeometricInfo_NURBS(Mesh& mesh)
  {
    allocateGeometricInfo_hexahedron(mesh);
    return 0;
  }
   
  int computeGeometricInfo_NURBS(Mesh& mesh)
  {
    computeGeometricInfo_hexahedron(mesh);
    return 0;
  }
 
 
  int allocateGeometricInfo_triangle(Mesh& mesh)
  {
    mesh.elementDiametersArray = new double[mesh.nElements_global];
    mesh.elementInnerDiametersArray = new double[mesh.nElements_global];
    mesh.elementBoundaryDiametersArray = new double[mesh.nElementBoundaries_global];
    mesh.elementBarycentersArray  = new double[mesh.nElements_global*3];
    mesh.elementBoundaryBarycentersArray = new double[mesh.nElementBoundaries_global*3];
    mesh.nodeDiametersArray = new double[mesh.nNodes_global];
    mesh.nodeSupportArray = new double[mesh.nNodes_global];
    return 0;
  }
 
  int computeGeometricInfo_triangle(Mesh& mesh)
  {
    memset(mesh.elementDiametersArray,0,mesh.nElements_global*sizeof(double));
    memset(mesh.elementInnerDiametersArray,0,mesh.nElements_global*sizeof(double));
    memset(mesh.elementBoundaryDiametersArray,0,mesh.nElementBoundaries_global*sizeof(double));
    memset(mesh.elementBarycentersArray,0,mesh.nElements_global*3*sizeof(double));
    memset(mesh.elementBoundaryBarycentersArray,0,mesh.nElementBoundaries_global*3*sizeof(double));
    memset(mesh.nodeDiametersArray,0,mesh.nNodes_global*sizeof(double));
    memset(mesh.nodeSupportArray,0,mesh.nNodes_global*sizeof(double));
    mesh.hMin = edgeLength(mesh.elementNodesArray[0],
                           mesh.elementNodesArray[1],
                           mesh.nodeArray);
    for (int ebN=0;ebN<mesh.nElementBoundaries_global;ebN++)
      mesh.elementBoundaryDiametersArray[ebN] = edgeLength(mesh.elementBoundaryNodesArray[ebN*2+0],
                                                           mesh.elementBoundaryNodesArray[ebN*2+1],
                                                           mesh.nodeArray);
    const double nNperElemInv = 1.0/double(mesh.nNodes_element);
    const double nNperElemBInv= 1.0/double(mesh.nNodes_elementBoundary);
 
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        double area=0.0, perimeter=0.0,h=0.0, hMin=0.0,hMax=0.0;
        area = triangleArea(mesh.elementNodesArray[eN*3 + 0],
                            mesh.elementNodesArray[eN*3 + 1],
                            mesh.elementNodesArray[eN*3 + 2],
                            mesh.nodeArray);
        //hMax
        for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
          for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
            {
              h = edgeLength(mesh.elementNodesArray[eN*3+nN_L],
                             mesh.elementNodesArray[eN*3+nN_R],
                             mesh.nodeArray);
              perimeter += h;
              hMax = fmax(hMax,h);
            } 
        hMin = 4.0*area/perimeter; 
        mesh.elementInnerDiametersArray[eN] = hMin;
        mesh.elementDiametersArray[eN] = hMax;
        mesh.volume += area;
        mesh.sigmaMax = fmax(hMax/hMin,mesh.sigmaMax);
        mesh.h = fmax(hMax,mesh.h);
        mesh.hMin = fmin(hMin,mesh.hMin);
 
        mesh.elementBarycentersArray[eN*3 + 0] = 0.0; 
        mesh.elementBarycentersArray[eN*3 + 1] = 0.0; 
        mesh.elementBarycentersArray[eN*3 + 2] = 0.0; 
        for (int nN=0;nN<mesh.nNodes_element;nN++)
          {
            mesh.elementBarycentersArray[eN*3 + 0] += 
              nNperElemInv*mesh.nodeArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]*3 + 0];
            mesh.elementBarycentersArray[eN*3 + 1] += 
              nNperElemInv*mesh.nodeArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]*3 + 1];
            mesh.elementBarycentersArray[eN*3 + 2] += 
              nNperElemInv*mesh.nodeArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]*3 + 2];
            mesh.nodeDiametersArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]] += hMax*area;
            mesh.nodeSupportArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]] += area;
          }
 
      }
 
    for (int ebN=0;ebN<mesh.nElementBoundaries_global;ebN++)
      {
        mesh.elementBoundaryBarycentersArray[ebN*3 + 0] = 0.0;
        mesh.elementBoundaryBarycentersArray[ebN*3 + 1] = 0.0;
        mesh.elementBoundaryBarycentersArray[ebN*3 + 2] = 0.0;
 
        for (int nN=0;nN<mesh.nNodes_elementBoundary;nN++)
          {
            mesh.elementBoundaryBarycentersArray[ebN*3 + 0] += 
              nNperElemBInv*mesh.nodeArray[mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN]*3 + 0];
            mesh.elementBoundaryBarycentersArray[ebN*3 + 1] += 
              nNperElemBInv*mesh.nodeArray[mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN]*3 + 1];
            mesh.elementBoundaryBarycentersArray[ebN*3 + 2] += 
              nNperElemBInv*mesh.nodeArray[mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN]*3 + 2];
          }
      }
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      {
        mesh.nodeDiametersArray[nN] /= mesh.nodeSupportArray[nN];
      }
    //printf("volume = %12.5e \n",mesh.volume);
    //printf("h = %12.5e \n",mesh.h);
    //printf("hMin = %12.5e \n",mesh.hMin);
    //printf("sigmaMax = %12.5e \n",mesh.sigmaMax);
    return 0;
  }
  
  int allocateGeometricInfo_quadrilateral(Mesh& mesh)
  {
    mesh.elementDiametersArray = new double[mesh.nElements_global];
    mesh.elementInnerDiametersArray = new double[mesh.nElements_global];
    mesh.elementBoundaryDiametersArray = new double[mesh.nElementBoundaries_global];
    mesh.nodeDiametersArray = new double[mesh.nNodes_global];
    mesh.nodeSupportArray = new double[mesh.nNodes_global];
    return 0;
  }
 
  int computeGeometricInfo_quadrilateral(Mesh& mesh)
  {
    memset(mesh.elementDiametersArray,0,mesh.nElements_global*sizeof(double));
    memset(mesh.elementInnerDiametersArray,0,mesh.nElements_global*sizeof(double));
    memset(mesh.elementBoundaryDiametersArray,0,mesh.nElementBoundaries_global*sizeof(double));
    memset(mesh.nodeDiametersArray,0,mesh.nNodes_global*sizeof(double));
    memset(mesh.nodeSupportArray,0,mesh.nNodes_global*sizeof(double));
    mesh.hMin = edgeLength(mesh.elementNodesArray[0],
                           mesh.elementNodesArray[1],
                           mesh.nodeArray);
    for (int ebN=0;ebN<mesh.nElementBoundaries_global;ebN++)
      mesh.elementBoundaryDiametersArray[ebN] = edgeLength(mesh.elementBoundaryNodesArray[ebN*2+0],
                                                           mesh.elementBoundaryNodesArray[ebN*2+1],
                                                           mesh.nodeArray);
    const double nNperElemInv = 1.0/double(mesh.nNodes_element);
    const double nNperElemBInv= 1.0/double(mesh.nNodes_elementBoundary);
 
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        double area=0.0, perimeter=0.0,h=0.0, hMin=1.0e6,hMax=0.0;
        area = triangleArea(mesh.elementNodesArray[eN*4 + 0],
                            mesh.elementNodesArray[eN*4 + 1],
                            mesh.elementNodesArray[eN*4 + 2],
                            mesh.nodeArray) +
          triangleArea(mesh.elementNodesArray[eN*4 + 0],
                       mesh.elementNodesArray[eN*4 + 2],
                       mesh.elementNodesArray[eN*4 + 3],
                       mesh.nodeArray);
          
        //hMax &&  hMin
        for (int nN_L=0;nN_L<mesh.nNodes_element;nN_L++)
          for (int nN_R=nN_L+1;nN_R<mesh.nNodes_element;nN_R++)
            {
              h = edgeLength(mesh.elementNodesArray[eN*4+nN_L],
                             mesh.elementNodesArray[eN*4+nN_R],
                             mesh.nodeArray);
              perimeter += h;
              hMax = fmax(hMax,h);
              hMin = fmin(hMin,h);
            } 
        mesh.elementInnerDiametersArray[eN] = hMin;
        mesh.elementDiametersArray[eN] = hMax;
        mesh.volume += area;
        mesh.sigmaMax = fmax(hMax/hMin,mesh.sigmaMax);
        mesh.h = fmax(hMax,mesh.h);
        mesh.hMin = fmin(hMin,mesh.hMin);
        for (int nN=0;nN<mesh.nNodes_element;nN++)
          {
            mesh.nodeDiametersArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]] += hMax*area;
            mesh.nodeSupportArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]] += area;
          }
      }
 
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      {
        mesh.nodeDiametersArray[nN] /= mesh.nodeSupportArray[nN];
      }
    return 0;
  }
  
  int allocateGeometricInfo_edge(Mesh& mesh)
  {
    mesh.elementDiametersArray = new double[mesh.nElements_global];
    mesh.elementInnerDiametersArray = new double[mesh.nElements_global];
    mesh.elementBoundaryDiametersArray = new double[mesh.nElementBoundaries_global];
    mesh.elementBarycentersArray  = new double[mesh.nElements_global*3];
    mesh.elementBoundaryBarycentersArray = new double[mesh.nElementBoundaries_global*3];
    mesh.nodeDiametersArray = new double[mesh.nNodes_global];
    mesh.nodeSupportArray = new double[mesh.nNodes_global];
    return 0;
  }
 
  int computeGeometricInfo_edge(Mesh& mesh)
  {
    memset(mesh.elementDiametersArray,0,mesh.nElements_global*sizeof(double));
    memset(mesh.elementInnerDiametersArray,0,mesh.nElements_global*sizeof(double));
    memset(mesh.elementBoundaryDiametersArray,0,mesh.nElementBoundaries_global*sizeof(double));
    memset(mesh.elementBarycentersArray,0,mesh.nElements_global*3*sizeof(double));
    memset(mesh.elementBoundaryBarycentersArray,0,mesh.nElementBoundaries_global*3*sizeof(double));
    memset(mesh.nodeDiametersArray,0,mesh.nNodes_global*sizeof(double));
    memset(mesh.nodeSupportArray,0,mesh.nNodes_global*sizeof(double));
    mesh.hMin = edgeLength(mesh.elementNodesArray[0],
                           mesh.elementNodesArray[1],
                           mesh.nodeArray);
    const double nNperElemInv = 1.0/double(mesh.nNodes_element);
    const double nNperElemBInv= 1.0/double(mesh.nNodes_elementBoundary);
    for (int eN=0;eN<mesh.nElements_global;eN++)
      {
        mesh.elementDiametersArray[eN] = edgeLength(mesh.elementNodesArray[eN*2+0],
                                                    mesh.elementNodesArray[eN*2+1],
                                                    mesh.nodeArray);
        mesh.elementInnerDiametersArray[eN] = mesh.elementDiametersArray[eN];
        mesh.h = fmax(mesh.h,mesh.elementDiametersArray[eN]);
        mesh.hMin = fmin(mesh.hMin,mesh.elementDiametersArray[eN]);
        mesh.volume += mesh.elementDiametersArray[eN];
 
        mesh.elementBarycentersArray[eN*3 + 0] = 0.0; 
        mesh.elementBarycentersArray[eN*3 + 1] = 0.0; 
        mesh.elementBarycentersArray[eN*3 + 2] = 0.0; 
        for (int nN=0;nN<mesh.nNodes_element;nN++)
          {
            mesh.elementBarycentersArray[eN*3 + 0] += 
              nNperElemInv*mesh.nodeArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]*3 + 0];
            mesh.elementBarycentersArray[eN*3 + 1] += 
              nNperElemInv*mesh.nodeArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]*3 + 1];
            mesh.elementBarycentersArray[eN*3 + 2] += 
              nNperElemInv*mesh.nodeArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]*3 + 2];
            mesh.nodeDiametersArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]] += mesh.elementDiametersArray[eN]*mesh.elementDiametersArray[eN];
            mesh.nodeSupportArray[mesh.elementNodesArray[eN*mesh.nNodes_element + nN]] += mesh.elementDiametersArray[eN];
          }
      }
    for (int ebN=0;ebN<mesh.nElementBoundaries_global;ebN++)
      {
        mesh.elementBoundaryDiametersArray[ebN] = 1.0;
 
        mesh.elementBoundaryBarycentersArray[ebN*3 + 0] = 0.0;
        mesh.elementBoundaryBarycentersArray[ebN*3 + 1] = 0.0;
        mesh.elementBoundaryBarycentersArray[ebN*3 + 2] = 0.0;
 
        for (int nN=0;nN<mesh.nNodes_elementBoundary;nN++)
          {
            mesh.elementBoundaryBarycentersArray[ebN*3 + 0] += 
              nNperElemBInv*mesh.nodeArray[mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN]*3 + 0];
            mesh.elementBoundaryBarycentersArray[ebN*3 + 1] += 
              nNperElemBInv*mesh.nodeArray[mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN]*3 + 1];
            mesh.elementBoundaryBarycentersArray[ebN*3 + 2] += 
              nNperElemBInv*mesh.nodeArray[mesh.elementBoundaryNodesArray[ebN*mesh.nNodes_elementBoundary+nN]*3 + 2];
          }
      }    
    for (int nN=0;nN<mesh.nNodes_global;nN++)
      {
        mesh.nodeDiametersArray[nN] /= mesh.nodeSupportArray[nN];
      }
    mesh.sigmaMax = 1.0;
    //printf("volume = %12.5e \n",mesh.volume);
    //printf("h = %12.5e \n",mesh.h);
    //printf("hMin = %12.5e \n",mesh.hMin);
    //printf("sigmaMax = %12.5e \n",mesh.sigmaMax);
    return 0;
  }
  //allocate node and element-node connectivity tables so they can be filled in externally
  int allocateNodeAndElementNodeDataStructures(Mesh& mesh, int nElements_global, int nNodes_global, int nNodes_element)
  {
    assert(!mesh.nodeArray);
    assert(!mesh.elementNodesArray);
    assert(!mesh.elementMaterialTypes);
    assert(!mesh.nodeMaterialTypes);
 
    mesh.nElements_global = nElements_global;
    mesh.nNodes_global    = nNodes_global;
    mesh.nNodes_element   = nNodes_element;
 
    mesh.nodeArray         = new double[mesh.nNodes_global*3];
    mesh.nodeMaterialTypes = new int[mesh.nNodes_global];
    mesh.elementNodesArray = new int[mesh.nElements_global*mesh.nNodes_element];
    mesh.elementMaterialTypes = new int[mesh.nElements_global];
 
    return 0;
  }
 
  //mwftodo get global refinement to preserve element boundary type   
  int globallyRefineEdgeMesh(const int& nLevels, Mesh& mesh, MultilevelMesh& multilevelMesh, bool averageNewNodeFlags)
  {
    using namespace  std;
    multilevelMesh.nLevels = nLevels;
    multilevelMesh.meshArray = new Mesh[nLevels];
    for(int i=0;i<nLevels;i++)
      initializeMesh(multilevelMesh.meshArray[i]);
    multilevelMesh.elementChildrenArray = new int*[nLevels];
    multilevelMesh.elementChildrenOffsets = new int*[nLevels];
    multilevelMesh.elementParentsArray = new int*[nLevels];
    multilevelMesh.meshArray[0] = mesh; //shallow copy
    for(int i=1;i<nLevels;i++)
      {
        //cout<<"Refinement Level "<<i<<endl;
        set<Node> newNodeSet;
        set<Node>::iterator nodeItr;
        pair<set<Node>::iterator,bool> ret;
        multilevelMesh.meshArray[i].nNodes_element=2;
        //2 children per parent
        multilevelMesh.meshArray[i].nElements_global  = 2*multilevelMesh.meshArray[i-1].nElements_global;
        multilevelMesh.meshArray[i].elementNodesArray = new int[multilevelMesh.meshArray[i].nElements_global*2];
        multilevelMesh.elementChildrenArray[i-1]      = new int[multilevelMesh.meshArray[i-1].nElements_global*2];
        multilevelMesh.elementChildrenOffsets[i-1]  = new int[multilevelMesh.meshArray[i-1].nElements_global+1];
        multilevelMesh.elementParentsArray[i]         = new int[multilevelMesh.meshArray[i].nElements_global];
        multilevelMesh.meshArray[i].elementMaterialTypes = new int[multilevelMesh.meshArray[i].nElements_global];
        int nN_new = multilevelMesh.meshArray[i-1].nNodes_global;
        multilevelMesh.elementChildrenOffsets[i-1][0] = 0;
        for(int eN_parent=0,eN=0;eN_parent<multilevelMesh.meshArray[i-1].nElements_global;eN_parent++)
          {
            multilevelMesh.elementChildrenOffsets[i-1][eN_parent+1] = multilevelMesh.elementChildrenOffsets[i-1][eN_parent]+2;
            multilevelMesh.elementChildrenArray[i-1][2*eN_parent + 0 ] = eN + 0;
            multilevelMesh.elementChildrenArray[i-1][2*eN_parent + 1 ] = eN + 1;
            multilevelMesh.elementParentsArray[i][eN + 0] = eN_parent;
            multilevelMesh.elementParentsArray[i][eN + 1] = eN_parent;
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 0] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 1] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            Node midpoints[1];
            midpoint(multilevelMesh.meshArray[i-1].nodeArray + multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*2 + 0]*3,
                     multilevelMesh.meshArray[i-1].nodeArray + multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*2 + 1]*3,
                     midpoints[0]);
            midpoints[0].nN = nN_new;
            if (multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*2 + 0]] ==
                multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*2 + 1]])
              midpoints[0].flag = multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*2 + 0]];
            else if (averageNewNodeFlags)
              midpoints[0].flag = 0.5*(multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*2 + 0]] +
                                       multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*2 + 1]]);
            else
              midpoints[0].flag = DEFAULT_NODE_MATERIAL;
 
            newNodeSet.insert(midpoints[0]);
            nN_new++;
            //the two new edges
            eN = newEdge(eN,multilevelMesh.meshArray[i].elementNodesArray,
                         multilevelMesh.meshArray[i-1].elementNodesArray[2*eN_parent+0],
                         midpoints[0].nN);
            eN = newEdge(eN,multilevelMesh.meshArray[i].elementNodesArray,
                         midpoints[0].nN,
                         multilevelMesh.meshArray[i-1].elementNodesArray[2*eN_parent+1]);
          }
        multilevelMesh.elementChildrenOffsets[i-1][multilevelMesh.meshArray[i-1].nElements_global] = multilevelMesh.elementChildrenOffsets[i-1][multilevelMesh.meshArray[i-1].nElements_global-1]+2;
        assert(unsigned(nN_new) == (newNodeSet.size()+multilevelMesh.meshArray[i-1].nNodes_global));
        multilevelMesh.meshArray[i].nNodes_global = nN_new;
        multilevelMesh.meshArray[i].nodeArray = new double[multilevelMesh.meshArray[i].nNodes_global*3];
        for(int nN=0;nN<multilevelMesh.meshArray[i-1].nNodes_global;nN++)
          {
            multilevelMesh.meshArray[i].nodeArray[nN*3+0] = multilevelMesh.meshArray[i-1].nodeArray[nN*3+0];
            multilevelMesh.meshArray[i].nodeArray[nN*3+1] = multilevelMesh.meshArray[i-1].nodeArray[nN*3+1];
            multilevelMesh.meshArray[i].nodeArray[nN*3+2] = multilevelMesh.meshArray[i-1].nodeArray[nN*3+2];
          }
        for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
          {
            multilevelMesh.meshArray[i].nodeArray[nodeItr->nN*3+0] = nodeItr->x;
            multilevelMesh.meshArray[i].nodeArray[nodeItr->nN*3+1] = nodeItr->y;
            multilevelMesh.meshArray[i].nodeArray[nodeItr->nN*3+2] = nodeItr->z;
          }
        multilevelMesh.meshArray[i].nodeMaterialTypes = new int[multilevelMesh.meshArray[i].nNodes_global];
        std::copy(multilevelMesh.meshArray[i-1].nodeMaterialTypes,
                  multilevelMesh.meshArray[i-1].nodeMaterialTypes+multilevelMesh.meshArray[i-1].nNodes_global,
                  multilevelMesh.meshArray[i].nodeMaterialTypes);
        //new nodes get default material type, should be set on interior and 
        //boundary in constructElementBoundaryElementsArray_*
        if (multilevelMesh.meshArray[i].nNodes_global > multilevelMesh.meshArray[i-1].nNodes_global)
          memset(multilevelMesh.meshArray[i].nodeMaterialTypes+multilevelMesh.meshArray[i-1].nNodes_global,
                 DEFAULT_NODE_MATERIAL,
                 (multilevelMesh.meshArray[i].nNodes_global-multilevelMesh.meshArray[i-1].nNodes_global)*sizeof(int));
 
        for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
          {
            multilevelMesh.meshArray[i].nodeMaterialTypes[nodeItr->nN] = nodeItr->flag;
          }
      }
    return 0;
  }
 
  int globallyRefineTriangularMesh(const int& nLevels, Mesh& mesh, MultilevelMesh& multilevelMesh, bool averageNewNodeFlags)
  {
    using namespace  std;
    multilevelMesh.nLevels = nLevels;
    multilevelMesh.meshArray = new Mesh[nLevels];
    for(int i=0;i<nLevels;i++)
      initializeMesh(multilevelMesh.meshArray[i]);
    multilevelMesh.elementChildrenArray = new int*[nLevels];
    multilevelMesh.elementChildrenOffsets = new int*[nLevels];
    multilevelMesh.elementParentsArray = new int*[nLevels];
    multilevelMesh.meshArray[0] = mesh; //shallow copy
    for(int i=1;i<nLevels;i++)
      {
        //cout<<"Refinement Level "<<i<<endl;
        set<Node> newNodeSet;
        set<Node>::iterator nodeItr;
        pair<set<Node>::iterator,bool> ret;
        multilevelMesh.meshArray[i].nNodes_element=3;
        //4 children per parent
        multilevelMesh.meshArray[i].nElements_global  = 4*multilevelMesh.meshArray[i-1].nElements_global;
        multilevelMesh.meshArray[i].elementNodesArray = new int[multilevelMesh.meshArray[i].nElements_global*3];
        multilevelMesh.elementChildrenArray[i-1]      = new int[multilevelMesh.meshArray[i-1].nElements_global*4];
        multilevelMesh.elementChildrenOffsets[i-1]  = new int[multilevelMesh.meshArray[i-1].nElements_global+1];
        multilevelMesh.elementParentsArray[i]         = new int[multilevelMesh.meshArray[i].nElements_global];
        multilevelMesh.meshArray[i].elementMaterialTypes = new int[multilevelMesh.meshArray[i].nElements_global];
        multilevelMesh.elementChildrenOffsets[i-1][0] = 0;
        int nN_new = multilevelMesh.meshArray[i-1].nNodes_global;
        for(int eN_parent=0,eN=0;eN_parent<multilevelMesh.meshArray[i-1].nElements_global;eN_parent++)
          {
            multilevelMesh.elementChildrenOffsets[i-1][eN_parent+1] = multilevelMesh.elementChildrenOffsets[i-1][eN_parent] + 4;
            multilevelMesh.elementChildrenArray[i-1][4*eN_parent + 0 ] = eN + 0;
            multilevelMesh.elementChildrenArray[i-1][4*eN_parent + 1 ] = eN + 1;
            multilevelMesh.elementChildrenArray[i-1][4*eN_parent + 2 ] = eN + 2;
            multilevelMesh.elementChildrenArray[i-1][4*eN_parent + 3 ] = eN + 3;
            multilevelMesh.elementParentsArray[i][eN + 0] = eN_parent;
            multilevelMesh.elementParentsArray[i][eN + 1] = eN_parent;
            multilevelMesh.elementParentsArray[i][eN + 2] = eN_parent;
            multilevelMesh.elementParentsArray[i][eN + 3] = eN_parent;
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 0] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 1] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 2] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 3] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            Node midpoints[3];
            for(int nN_element_0=0,nN_midpoint=0;nN_element_0<mesh.nNodes_element;nN_element_0++)
              for(int nN_element_1=nN_element_0+1;nN_element_1<mesh.nNodes_element;nN_element_1++,nN_midpoint++)
                {
                  midpoint(multilevelMesh.meshArray[i-1].nodeArray + multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*3 + nN_element_0]*3,
                           multilevelMesh.meshArray[i-1].nodeArray + multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*3 + nN_element_1]*3,
                           midpoints[nN_midpoint]);
                  nodeItr = newNodeSet.find(midpoints[nN_midpoint]);
                  if(nodeItr == newNodeSet.end())
                    {
                      midpoints[nN_midpoint].nN = nN_new;
                      if (multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*3 + nN_element_0]] ==
                          multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*3 + nN_element_1]])
                        midpoints[nN_midpoint].flag = multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*3 + nN_element_0]];
                      else if (averageNewNodeFlags)
                        midpoints[nN_midpoint].flag = 0.5*(multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*3 + nN_element_0]] +
                                                           multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*3 + nN_element_1]]);
                      else
                        midpoints[nN_midpoint].flag = DEFAULT_NODE_MATERIAL;
 
                      newNodeSet.insert(midpoints[nN_midpoint]);
                      nN_new++;
                    }
                  else
                    midpoints[nN_midpoint].nN = nodeItr->nN;
                }
            //the triangles formed by chopping the points off the parent
            eN = newTriangle(eN,multilevelMesh.meshArray[i].elementNodesArray,
                             multilevelMesh.meshArray[i-1].elementNodesArray[3*eN_parent+0],
                             midpoints[0].nN,
                             midpoints[1].nN);
            eN = newTriangle(eN,multilevelMesh.meshArray[i].elementNodesArray,
                             multilevelMesh.meshArray[i-1].elementNodesArray[3*eN_parent+1],
                             midpoints[0].nN,
                             midpoints[2].nN);
            eN = newTriangle(eN,multilevelMesh.meshArray[i].elementNodesArray,
                             multilevelMesh.meshArray[i-1].elementNodesArray[3*eN_parent+2],
                             midpoints[1].nN,
                             midpoints[2].nN);
            eN = newTriangle(eN,multilevelMesh.meshArray[i].elementNodesArray,
                             midpoints[0].nN,
                             midpoints[1].nN,
                             midpoints[2].nN);
          }
        multilevelMesh.elementChildrenOffsets[i-1][multilevelMesh.meshArray[i-1].nElements_global] = multilevelMesh.elementChildrenOffsets[i-1][multilevelMesh.meshArray[i-1].nElements_global-1]+4;
        assert(unsigned(nN_new) == (newNodeSet.size()+multilevelMesh.meshArray[i-1].nNodes_global));
        multilevelMesh.meshArray[i].nNodes_global = nN_new;
        multilevelMesh.meshArray[i].nodeArray = new double[multilevelMesh.meshArray[i].nNodes_global*3];
        for(int nN=0;nN<multilevelMesh.meshArray[i-1].nNodes_global;nN++)
          {
            multilevelMesh.meshArray[i].nodeArray[nN*3+0] = multilevelMesh.meshArray[i-1].nodeArray[nN*3+0];
            multilevelMesh.meshArray[i].nodeArray[nN*3+1] = multilevelMesh.meshArray[i-1].nodeArray[nN*3+1];
            multilevelMesh.meshArray[i].nodeArray[nN*3+2] = multilevelMesh.meshArray[i-1].nodeArray[nN*3+2];
          }
        for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
          {
            multilevelMesh.meshArray[i].nodeArray[nodeItr->nN*3+0] = nodeItr->x;
            multilevelMesh.meshArray[i].nodeArray[nodeItr->nN*3+1] = nodeItr->y;
            multilevelMesh.meshArray[i].nodeArray[nodeItr->nN*3+2] = nodeItr->z;
          }
        multilevelMesh.meshArray[i].nodeMaterialTypes = new int[multilevelMesh.meshArray[i].nNodes_global];
        std::copy(multilevelMesh.meshArray[i-1].nodeMaterialTypes,
                  multilevelMesh.meshArray[i-1].nodeMaterialTypes+multilevelMesh.meshArray[i-1].nNodes_global,
                  multilevelMesh.meshArray[i].nodeMaterialTypes);
        //new nodes get default material type, should be set on interior and 
        //boundary in constructElementBoundaryElementsArray_*
        if (multilevelMesh.meshArray[i].nNodes_global > multilevelMesh.meshArray[i-1].nNodes_global)
          memset(multilevelMesh.meshArray[i].nodeMaterialTypes+multilevelMesh.meshArray[i-1].nNodes_global,
                 DEFAULT_NODE_MATERIAL,
                 (multilevelMesh.meshArray[i].nNodes_global-multilevelMesh.meshArray[i-1].nNodes_global)*sizeof(int));
        
        for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
          {
            //mwf hack 
            //std::cout<<"refined triangular mesh: nN parent= "<<multilevelMesh.meshArray[i-1].nNodes_global<<" nN child= "<<multilevelMesh.meshArray[i].nNodes_global
            //         <<" nodeItr->nN,flag,x= "<<nodeItr->nN<<" , "<<nodeItr->flag <<" , " << nodeItr->x <<" , "<< nodeItr->y << " , "<< nodeItr->z << std::endl;
            multilevelMesh.meshArray[i].nodeMaterialTypes[nodeItr->nN] = nodeItr->flag;
          }
        
      }
    return 0;
  }
  
  int globallyRefineQuadrilateralMesh(const int& nLevels, Mesh& mesh, MultilevelMesh& multilevelMesh, bool averageNewNodeFlags)
  {
    using namespace  std;
    multilevelMesh.nLevels = nLevels;
    multilevelMesh.meshArray = new Mesh[nLevels];
    for(int i=0;i<nLevels;i++)
      initializeMesh(multilevelMesh.meshArray[i]);
    multilevelMesh.elementChildrenArray = new int*[nLevels];
    multilevelMesh.elementChildrenOffsets = new int*[nLevels];
    multilevelMesh.elementParentsArray = new int*[nLevels];
    multilevelMesh.meshArray[0] = mesh; //shallow copy
    for(int i=1;i<nLevels;i++)
      {
        std::cout<<"Quad refinement not imlemented "<<i<<endl;
        assert(false);
      }
    return 0;
  }
 
  int globallyRefineTetrahedralMesh(const int& nLevels, Mesh& mesh, MultilevelMesh& multilevelMesh, bool averageNewNodeFlags)
  {
    using namespace std;
    multilevelMesh.nLevels = nLevels;
    multilevelMesh.meshArray = new Mesh[nLevels];
    for(int i=0;i<nLevels;i++)
      initializeMesh(multilevelMesh.meshArray[i]);
    multilevelMesh.elementChildrenArray = new int*[nLevels];
    multilevelMesh.elementChildrenOffsets = new int*[nLevels];
    multilevelMesh.elementParentsArray = new int*[nLevels];
    multilevelMesh.meshArray[0] = mesh; //shallow copy
    for(int i=1;i<nLevels;i++)
      {
        //cout<<"Refinement Level "<<i<<endl;
        set<Node> newNodeSet;
        set<Node>::iterator nodeItr;
        pair<set<Node>::iterator,bool> ret;
        multilevelMesh.meshArray[i].nNodes_element=4;
        //8 children per parent
        multilevelMesh.meshArray[i].nElements_global  = 8*multilevelMesh.meshArray[i-1].nElements_global;
        multilevelMesh.meshArray[i].elementNodesArray = new int[multilevelMesh.meshArray[i].nElements_global*4];
        multilevelMesh.elementChildrenArray[i-1]      = new int[multilevelMesh.meshArray[i-1].nElements_global*8];
        multilevelMesh.elementChildrenOffsets[i-1]  = new int[multilevelMesh.meshArray[i-1].nElements_global+1];
        multilevelMesh.elementParentsArray[i]         = new int[multilevelMesh.meshArray[i].nElements_global];
        multilevelMesh.meshArray[i].elementMaterialTypes = new int[multilevelMesh.meshArray[i].nElements_global];
        multilevelMesh.elementChildrenOffsets[i-1][0] = 0;
        int nN_new = multilevelMesh.meshArray[i-1].nNodes_global;
        for(int eN_parent=0,eN=0;eN_parent<multilevelMesh.meshArray[i-1].nElements_global;eN_parent++)
          {
            multilevelMesh.elementChildrenOffsets[i-1][eN_parent+1] = multilevelMesh.elementChildrenOffsets[i-1][eN_parent]+8;
            multilevelMesh.elementChildrenArray[i-1][8*eN_parent + 0 ] = eN + 0;
            multilevelMesh.elementChildrenArray[i-1][8*eN_parent + 1 ] = eN + 1;
            multilevelMesh.elementChildrenArray[i-1][8*eN_parent + 2 ] = eN + 2;
            multilevelMesh.elementChildrenArray[i-1][8*eN_parent + 3 ] = eN + 3;
            multilevelMesh.elementChildrenArray[i-1][8*eN_parent + 4 ] = eN + 4;
            multilevelMesh.elementChildrenArray[i-1][8*eN_parent + 5 ] = eN + 5;
            multilevelMesh.elementChildrenArray[i-1][8*eN_parent + 6 ] = eN + 6;
            multilevelMesh.elementChildrenArray[i-1][8*eN_parent + 7 ] = eN + 7; 
            multilevelMesh.elementParentsArray[i][eN + 0] = eN_parent;
            multilevelMesh.elementParentsArray[i][eN + 1] = eN_parent;
            multilevelMesh.elementParentsArray[i][eN + 2] = eN_parent;
            multilevelMesh.elementParentsArray[i][eN + 3] = eN_parent;
            multilevelMesh.elementParentsArray[i][eN + 4] = eN_parent;
            multilevelMesh.elementParentsArray[i][eN + 5] = eN_parent;
            multilevelMesh.elementParentsArray[i][eN + 6] = eN_parent;
            multilevelMesh.elementParentsArray[i][eN + 7] = eN_parent; 
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 0] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 1] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 2] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 3] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 4] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 5] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 6] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent];
            multilevelMesh.meshArray[i].elementMaterialTypes[eN + 7] = multilevelMesh.meshArray[i-1].elementMaterialTypes[eN_parent]; 
            
            Node midpoints[6];
            double mind;
            int mindN;
            for(int nN_element_0=0,nN_midpoint=0;nN_element_0<mesh.nNodes_element;nN_element_0++)
              for(int nN_element_1=nN_element_0+1;nN_element_1<mesh.nNodes_element;nN_element_1++,nN_midpoint++)
                {
                  midpoint(multilevelMesh.meshArray[i-1].nodeArray + multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*4 + nN_element_0]*3,
                           multilevelMesh.meshArray[i-1].nodeArray + multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*4 + nN_element_1]*3,
                           midpoints[nN_midpoint]);
                  nodeItr = newNodeSet.find(midpoints[nN_midpoint]);
                  if(nodeItr == newNodeSet.end())
                    {
                      midpoints[nN_midpoint].nN = nN_new;
                      if (multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*4 + nN_element_0]] ==
                          multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*4 + nN_element_1]])
                        midpoints[nN_midpoint].flag = multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*4 + nN_element_0]];
                      else if (averageNewNodeFlags)
                        midpoints[nN_midpoint].flag = 0.5*(multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*4 + nN_element_0]] +
                                                           multilevelMesh.meshArray[i-1].nodeMaterialTypes[multilevelMesh.meshArray[i-1].elementNodesArray[eN_parent*4 + nN_element_1]]);
                      else
                        midpoints[nN_midpoint].flag = DEFAULT_NODE_MATERIAL;
 
                      newNodeSet.insert(midpoints[nN_midpoint]);
                      nN_new++;
                    }
                  else
                    midpoints[nN_midpoint].nN = nodeItr->nN;
                }
            //the tets formed by chopping the points of the parent
            eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                multilevelMesh.meshArray[i-1].elementNodesArray[4*eN_parent+0],
                                midpoints[0].nN,
                                midpoints[1].nN,
                                midpoints[2].nN);
            eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                multilevelMesh.meshArray[i-1].elementNodesArray[4*eN_parent+1],
                                midpoints[0].nN,
                                midpoints[3].nN,
                                midpoints[4].nN);
            eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                multilevelMesh.meshArray[i-1].elementNodesArray[4*eN_parent+2],
                                midpoints[1].nN,
                                midpoints[3].nN,
                                midpoints[5].nN);
            eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                multilevelMesh.meshArray[i-1].elementNodesArray[4*eN_parent+3],
                                midpoints[2].nN,
                                midpoints[4].nN,
                                midpoints[5].nN);
            mind=edgeLength(midpoints[0],midpoints[5]);
            mindN=0;
            double dd;
            for(int dN=1;dN<3;dN++)
              {
                dd = edgeLength(midpoints[dN],midpoints[5-dN]);
                if (dd < mind)
                  {
                    mind = dd;
                    mindN = dN;
                  }
                else if (dd == mind)
                  {
                    if(midpoints[dN] < midpoints[mindN])
                      mindN = dN;
                    else if (!(midpoints[mindN] < midpoints[dN]) && (midpoints[5-dN] < midpoints[5-mindN]))
                      mindN = dN;
                  }
              }
            if(mindN == 0)
              {
                eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                    midpoints[0].nN, 
                                    midpoints[5].nN, 
                                    midpoints[2].nN, 
                                    midpoints[4].nN);
                eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                    midpoints[0].nN, 
                                    midpoints[5].nN, 
                                    midpoints[2].nN, 
                                    midpoints[1].nN);
                eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                    midpoints[0].nN, 
                                    midpoints[5].nN, 
                                    midpoints[1].nN, 
                                    midpoints[3].nN);
                eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                    midpoints[0].nN, 
                                    midpoints[5].nN, 
                                    midpoints[3].nN, 
                                    midpoints[4].nN);
              }
            else if (mindN == 1)
              {
                eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                    midpoints[1].nN, 
                                    midpoints[4].nN, 
                                    midpoints[2].nN, 
                                    midpoints[5].nN);
                eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                    midpoints[1].nN, 
                                    midpoints[4].nN, 
                                    midpoints[5].nN, 
                                    midpoints[3].nN);
                eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                    midpoints[1].nN, 
                                    midpoints[4].nN, 
                                    midpoints[3].nN, 
                                    midpoints[0].nN);
                eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                    midpoints[1].nN, 
                                    midpoints[4].nN, 
                                    midpoints[0].nN, 
                                    midpoints[2].nN);
              }
            else
              {
                eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                    midpoints[2].nN, 
                                    midpoints[3].nN, 
                                    midpoints[0].nN, 
                                    midpoints[4].nN);
                eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                    midpoints[2].nN, 
                                    midpoints[3].nN, 
                                    midpoints[4].nN, 
                                    midpoints[5].nN);
                eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                    midpoints[2].nN, 
                                    midpoints[3].nN, 
                                    midpoints[5].nN, 
                                    midpoints[1].nN);
                eN = newTetrahedron(eN,multilevelMesh.meshArray[i].elementNodesArray,
                                    midpoints[2].nN, 
                                    midpoints[3].nN, 
                                    midpoints[1].nN, 
                                    midpoints[0].nN);
              }
          }
        multilevelMesh.elementChildrenOffsets[i-1][multilevelMesh.meshArray[i-1].nElements_global] = multilevelMesh.elementChildrenOffsets[i-1][multilevelMesh.meshArray[i-1].nElements_global-1]+8;
        assert(unsigned(nN_new) == (newNodeSet.size()+multilevelMesh.meshArray[i-1].nNodes_global));
        multilevelMesh.meshArray[i].nNodes_global = nN_new;
        multilevelMesh.meshArray[i].nodeArray = new double[multilevelMesh.meshArray[i].nNodes_global*3];
        for(int nN=0;nN<multilevelMesh.meshArray[i-1].nNodes_global;nN++)
          {
            multilevelMesh.meshArray[i].nodeArray[nN*3+0] = multilevelMesh.meshArray[i-1].nodeArray[nN*3+0];
            multilevelMesh.meshArray[i].nodeArray[nN*3+1] = multilevelMesh.meshArray[i-1].nodeArray[nN*3+1];
            multilevelMesh.meshArray[i].nodeArray[nN*3+2] = multilevelMesh.meshArray[i-1].nodeArray[nN*3+2];
          }
        for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
          {
            multilevelMesh.meshArray[i].nodeArray[nodeItr->nN*3+0] = nodeItr->x;
            multilevelMesh.meshArray[i].nodeArray[nodeItr->nN*3+1] = nodeItr->y;
            multilevelMesh.meshArray[i].nodeArray[nodeItr->nN*3+2] = nodeItr->z;
          }
        //         //cout<<"re-ordeing nodes llllllllllllllllllll"<<endl;
        //         for (int eN=0;eN<multilevelMesh.meshArray[i].nElements_global;eN++)
        //           {
        //             register int n0,n1,n2,n3;
        //             register double t[3][3],det;
            
        //             n0 = multilevelMesh.meshArray[i].elementNodesArray[eN*4+0];
        //             n1 = multilevelMesh.meshArray[i].elementNodesArray[eN*4+1];
        //             n2 = multilevelMesh.meshArray[i].elementNodesArray[eN*4+2];
        //             n3 = multilevelMesh.meshArray[i].elementNodesArray[eN*4+3];
            
        //             t[0][0] = multilevelMesh.meshArray[i].nodeArray[n1*3+0] - multilevelMesh.meshArray[i].nodeArray[n0*3+0];
        //             t[0][1] = multilevelMesh.meshArray[i].nodeArray[n1*3+1] - multilevelMesh.meshArray[i].nodeArray[n0*3+1];
        //             t[0][2] = multilevelMesh.meshArray[i].nodeArray[n1*3+2] - multilevelMesh.meshArray[i].nodeArray[n0*3+2];
            
        //             t[1][0] = multilevelMesh.meshArray[i].nodeArray[n2*3+0] - multilevelMesh.meshArray[i].nodeArray[n0*3+0];
        //             t[1][1] = multilevelMesh.meshArray[i].nodeArray[n2*3+1] - multilevelMesh.meshArray[i].nodeArray[n0*3+1];
        //             t[1][2] = multilevelMesh.meshArray[i].nodeArray[n2*3+2] - multilevelMesh.meshArray[i].nodeArray[n0*3+2];
            
        //             t[2][0] = multilevelMesh.meshArray[i].nodeArray[n3*3+0] - multilevelMesh.meshArray[i].nodeArray[n0*3+0];
        //             t[2][1] = multilevelMesh.meshArray[i].nodeArray[n3*3+1] - multilevelMesh.meshArray[i].nodeArray[n0*3+1];
        //             t[2][2] = multilevelMesh.meshArray[i].nodeArray[n3*3+2] - multilevelMesh.meshArray[i].nodeArray[n0*3+2];
            
        //             det = t[0][0]*(t[1][1]*t[2][2] - t[1][2]*t[2][1]) -   
        //               t[0][1]*(t[1][0]*t[2][2] - t[1][2]*t[2][0]) +       
        //               t[0][2]*(t[1][0]*t[2][1] - t[1][1]*t[2][0]);
            
        //             if(det < 0.0)
        //               {
        //                 multilevelMesh.meshArray[i].elementNodesArray[eN*4+2] = n3;
        //                 multilevelMesh.meshArray[i].elementNodesArray[eN*4+3] = n2;
        //               }
        //             n0 = multilevelMesh.meshArray[i].elementNodesArray[eN*4+0];
        //             n1 = multilevelMesh.meshArray[i].elementNodesArray[eN*4+1];
        //             n2 = multilevelMesh.meshArray[i].elementNodesArray[eN*4+2];
        //             n3 = multilevelMesh.meshArray[i].elementNodesArray[eN*4+3];
            
        //             t[0][0] = multilevelMesh.meshArray[i].nodeArray[n1*3+0] - multilevelMesh.meshArray[i].nodeArray[n0*3+0];
        //             t[0][1] = multilevelMesh.meshArray[i].nodeArray[n1*3+1] - multilevelMesh.meshArray[i].nodeArray[n0*3+1];
        //             t[0][2] = multilevelMesh.meshArray[i].nodeArray[n1*3+2] - multilevelMesh.meshArray[i].nodeArray[n0*3+2];
            
        //             t[1][0] = multilevelMesh.meshArray[i].nodeArray[n2*3+0] - multilevelMesh.meshArray[i].nodeArray[n0*3+0];
        //             t[1][1] = multilevelMesh.meshArray[i].nodeArray[n2*3+1] - multilevelMesh.meshArray[i].nodeArray[n0*3+1];
        //             t[1][2] = multilevelMesh.meshArray[i].nodeArray[n2*3+2] - multilevelMesh.meshArray[i].nodeArray[n0*3+2];
            
        //             t[2][0] = multilevelMesh.meshArray[i].nodeArray[n3*3+0] - multilevelMesh.meshArray[i].nodeArray[n0*3+0];
        //             t[2][1] = multilevelMesh.meshArray[i].nodeArray[n3*3+1] - multilevelMesh.meshArray[i].nodeArray[n0*3+1];
        //             t[2][2] = multilevelMesh.meshArray[i].nodeArray[n3*3+2] - multilevelMesh.meshArray[i].nodeArray[n0*3+2];
            
        //             det = fabs(t[0][0]*(t[1][1]*t[2][2] - t[1][2]*t[2][1]) -   
        //                        t[0][1]*(t[1][0]*t[2][2] - t[1][2]*t[2][0]) +   
        //                        t[0][2]*(t[1][0]*t[2][1] - t[1][1]*t[2][0]));
        //             //cout<<"det "<<det<<endl;
            
        //           }
        //mwftodo need to come up with convention for assigning new node ids
        multilevelMesh.meshArray[i].nodeMaterialTypes = new int[multilevelMesh.meshArray[i].nNodes_global];
        std::copy(multilevelMesh.meshArray[i-1].nodeMaterialTypes,
                  multilevelMesh.meshArray[i-1].nodeMaterialTypes+multilevelMesh.meshArray[i-1].nNodes_global,
                  multilevelMesh.meshArray[i].nodeMaterialTypes);
        //new nodes get default material type, should be set on interior and 
        //boundary in constructElementBoundaryElementsArray_*
        if (multilevelMesh.meshArray[i].nNodes_global > multilevelMesh.meshArray[i-1].nNodes_global)
          memset(multilevelMesh.meshArray[i].nodeMaterialTypes+multilevelMesh.meshArray[i-1].nNodes_global,
                 DEFAULT_NODE_MATERIAL,
                 (multilevelMesh.meshArray[i].nNodes_global-multilevelMesh.meshArray[i-1].nNodes_global)*sizeof(int));
 
        for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
          {
            multilevelMesh.meshArray[i].nodeMaterialTypes[nodeItr->nN] = nodeItr->flag;
          }
      }
    return 0;
  }
  
  int globallyRefineHexahedralMesh(const int& nLevels, Mesh& mesh, MultilevelMesh& multilevelMesh, bool averageNewNodeFlags)
  {
    using namespace std;
    multilevelMesh.nLevels = nLevels;
    multilevelMesh.meshArray = new Mesh[nLevels];
    for(int i=0;i<nLevels;i++)
      initializeMesh(multilevelMesh.meshArray[i]);
    multilevelMesh.elementChildrenArray = new int*[nLevels];
    multilevelMesh.elementChildrenOffsets = new int*[nLevels];
    multilevelMesh.elementParentsArray = new int*[nLevels];
    multilevelMesh.meshArray[0] = mesh; //shallow copy
    for(int i=1;i<nLevels;i++)
      {
        std::cout<<"Hexahedron refinement not implemented"<<std::endl;
        assert(false);
      }
    return 0;
  }
 
  int assignElementBoundaryMaterialTypesFromParent(Mesh& parentMesh, Mesh& childMesh, const int* levelElementParentsArray,
                                                   const int& nSpace_global)
  {
    int failed = 0;
    for (int ebNI = 0; ebNI < childMesh.nInteriorElementBoundaries_global; ebNI++)
      {
        int ebN = childMesh.interiorElementBoundariesArray[ebNI];
        int eN_left  = childMesh.elementBoundaryElementsArray[ebN*2+0];
        int eN_right = childMesh.elementBoundaryElementsArray[ebN*2+1];
        int eN_left_parent = levelElementParentsArray[eN_left];
        int eN_right_parent= levelElementParentsArray[eN_right];
        if (eN_left_parent == eN_right_parent) //interior to element on coarser level
          childMesh.elementBoundaryMaterialTypes[ebN] = INTERIOR_ELEMENT_BOUNDARY_MATERIAL;
        else
          {
            //find local element boundary id on left parent
            int left_parent_ebN = -1;
            for (int ebN_left_parent_element = 0; ebN_left_parent_element < parentMesh.nElementBoundaries_element; 
                 ebN_left_parent_element++)
              {
                if (parentMesh.elementNeighborsArray[eN_left_parent*parentMesh.nElementBoundaries_element +
                                                     ebN_left_parent_element] 
                    == eN_right_parent)
                  {
                    left_parent_ebN = ebN_left_parent_element;
                    break;
                  }
              }
            assert(0 <= left_parent_ebN < parentMesh.nElementBoundaries_element);
            int ebN_parent = parentMesh.elementBoundariesArray[eN_left_parent*parentMesh.nElementBoundaries_element + 
                                                               left_parent_ebN];       
            childMesh.elementBoundaryMaterialTypes[ebN] = parentMesh.elementBoundaryMaterialTypes[ebN_parent];
          }
      }//interior element boundaries
    for (int ebNE = 0; ebNE < childMesh.nExteriorElementBoundaries_global; ebNE++)
      {
        int ebN = childMesh.exteriorElementBoundariesArray[ebNE];
        int eN  = childMesh.elementBoundaryElementsArray[ebN*2+0];
        int eN_parent = levelElementParentsArray[eN];
        int lastExteriorElementBoundaryOnParent=-1;
        int nExteriorElementBoundariesOnParent =0;
        for (int ebN_element = 0; ebN_element < parentMesh.nElementBoundaries_element; ebN_element++)
          {
            if (parentMesh.elementNeighborsArray[eN_parent*parentMesh.nElementBoundaries_element + ebN_element] < 0)
              {
                lastExteriorElementBoundaryOnParent = ebN_element;
                nExteriorElementBoundariesOnParent++;
              }
          }
        assert(nExteriorElementBoundariesOnParent > 0); assert(lastExteriorElementBoundaryOnParent >= 0);
        if (nExteriorElementBoundariesOnParent > 1)
          {
            //more than 1 face of parent is on exterior boundary, have to figure out which one
            //holds the current face on the child
            if (nSpace_global == 1)
              {
                const int ebN_element_child = childMesh.elementBoundaryLocalElementBoundariesArray[ebN*2+0];
                //across from ebN
                const int nN_child0         = childMesh.elementNodesArray[eN*childMesh.nNodes_element+ebN_element_child];
                //on ebN
                const int nN_child1         = childMesh.elementNodesArray[eN*childMesh.nNodes_element+((ebN_element_child+1)%2)];     
                const double n_child  = (childMesh.nodeArray[nN_child1*3+0]-childMesh.nodeArray[nN_child0*3+0])/
                  fabs(childMesh.nodeArray[nN_child1*3+0]-childMesh.nodeArray[nN_child0*3+0]);
 
                const int nN_parent0 = parentMesh.elementNodesArray[eN_parent*childMesh.nNodes_element + 0];
                const int nN_parent1 = parentMesh.elementNodesArray[eN_parent*childMesh.nNodes_element + 1];
                const double n_parent0  = (parentMesh.nodeArray[nN_parent1*3+0]-parentMesh.nodeArray[nN_parent0*3+0])/
                  fabs(parentMesh.nodeArray[nN_parent1*3+0]-parentMesh.nodeArray[nN_parent0*3+0]);
                const double n_parent1  = (parentMesh.nodeArray[nN_parent0*3+0]-parentMesh.nodeArray[nN_parent1*3+0])/
                  fabs(parentMesh.nodeArray[nN_parent0*3+0]-parentMesh.nodeArray[nN_parent1*3+0]);
                
                int ebN_parent = -1;
                if (fabs(n_child-n_parent1) < 1.0e-8)
                  ebN_parent = parentMesh.elementBoundariesArray[eN_parent*parentMesh.nElementBoundaries_element + 1];
                else if (fabs(n_child-n_parent0) < 1.0e-8)
                  ebN_parent = parentMesh.elementBoundariesArray[eN_parent*parentMesh.nElementBoundaries_element + 0];
                assert(ebN_parent >= 0);
                childMesh.elementBoundaryMaterialTypes[ebN] = parentMesh.elementBoundaryMaterialTypes[ebN_parent];
              }//1d
            else if (nSpace_global == 2)
              {
                //mwftodo assumes relavant coordinates are first two
                const int nN_ebN_child[2] = {childMesh.elementBoundaryNodesArray[2*ebN+0],childMesh.elementBoundaryNodesArray[2*ebN+1]};
                double n_child[2]   = {childMesh.nodeArray[nN_ebN_child[1]*3+1]-childMesh.nodeArray[nN_ebN_child[0]*3+1],
                                       childMesh.nodeArray[nN_ebN_child[0]*3+0]-childMesh.nodeArray[nN_ebN_child[1]*3+0]};
                double tmp = sqrt(n_child[0]*n_child[0] + n_child[1]*n_child[1]);
                assert(tmp > 0.0);
                n_child[0] /= tmp; n_child[1] /= tmp;
 
                int ebN_parent = -1;
                for (int ebN_element_parent = 0; ebN_element_parent < parentMesh.nElementBoundaries_element; ebN_element_parent++)
                  {
                    int ebN_cur = parentMesh.elementBoundariesArray[eN_parent*parentMesh.nElementBoundaries_element + ebN_element_parent];
                    const int nN_ebN_parent[2] = {parentMesh.elementBoundaryNodesArray[2*ebN_cur+0],
                                                  parentMesh.elementBoundaryNodesArray[2*ebN_cur+1]};
                    
                    double n_parent[2]   = {parentMesh.nodeArray[nN_ebN_parent[1]*3+1]-parentMesh.nodeArray[nN_ebN_parent[0]*3+1],
                                            parentMesh.nodeArray[nN_ebN_parent[0]*3+0]-parentMesh.nodeArray[nN_ebN_parent[1]*3+0]};
                    tmp = sqrt(n_parent[0]*n_parent[0] + n_parent[1]*n_parent[1]);
                    assert(tmp > 0.0);
                    n_parent[0] /= tmp; n_parent[1] /= tmp;
                    //looking for face on parent with unit normal in same direction as child face
                    //but don't enforce same outside inside convention on mesh orderings
                    tmp = n_parent[0]*n_child[0] + n_parent[1]*n_child[1];
                               
                    if (fabs(sqrt(tmp*tmp) - 1.0) < 1.0e-8)
                      {
                        ebN_parent = ebN_cur;
                        break;
                      }
                  }
                assert(ebN_parent >= 0);
                childMesh.elementBoundaryMaterialTypes[ebN] = parentMesh.elementBoundaryMaterialTypes[ebN_parent];
              }//2d
            else
              {
                const double * nN_ebN_child_0_x = childMesh.nodeArray + 3*childMesh.elementBoundaryNodesArray[3*ebN+0];
                const double * nN_ebN_child_1_x = childMesh.nodeArray + 3*childMesh.elementBoundaryNodesArray[3*ebN+1];
                const double * nN_ebN_child_2_x = childMesh.nodeArray + 3*childMesh.elementBoundaryNodesArray[3*ebN+2];
                double n_child[3] = {0.,0.,0.};
                //(nN_1-nN_0) x (nN_2-nN_0)
                n_child[0] = (nN_ebN_child_1_x[1]-nN_ebN_child_0_x[1])*(nN_ebN_child_2_x[2]-nN_ebN_child_0_x[2])
                  -(nN_ebN_child_1_x[2]-nN_ebN_child_0_x[2])*(nN_ebN_child_2_x[1]-nN_ebN_child_0_x[1]);
 
                n_child[1] = (nN_ebN_child_1_x[2]-nN_ebN_child_0_x[2])*(nN_ebN_child_2_x[0]-nN_ebN_child_0_x[0])
                  -(nN_ebN_child_1_x[0]-nN_ebN_child_0_x[0])*(nN_ebN_child_2_x[2]-nN_ebN_child_0_x[2]);
 
                n_child[2] = (nN_ebN_child_1_x[0]-nN_ebN_child_0_x[0])*(nN_ebN_child_2_x[1]-nN_ebN_child_0_x[1])
                  -(nN_ebN_child_1_x[1]-nN_ebN_child_0_x[1])*(nN_ebN_child_2_x[0]-nN_ebN_child_0_x[0]);
 
                double tmp = sqrt(n_child[0]*n_child[0] + n_child[1]*n_child[1] + n_child[2]*n_child[2]);
                assert(tmp > 0.0);
                n_child[0] /= tmp; n_child[1] /= tmp; n_child[2] /= tmp;
 
                int ebN_parent = -1;
                double n_parent[3] = {0.,0.,0.};
                for (int ebN_element_parent = 0; ebN_element_parent < parentMesh.nElementBoundaries_element; ebN_element_parent++)
                  {
                    int ebN_cur = parentMesh.elementBoundariesArray[eN_parent*parentMesh.nElementBoundaries_element + ebN_element_parent];
 
                    const double * nN_ebN_parent_0_x = parentMesh.nodeArray + 3*parentMesh.elementBoundaryNodesArray[3*ebN_cur+0];
                    const double * nN_ebN_parent_1_x = parentMesh.nodeArray + 3*parentMesh.elementBoundaryNodesArray[3*ebN_cur+1];
                    const double * nN_ebN_parent_2_x = parentMesh.nodeArray + 3*parentMesh.elementBoundaryNodesArray[3*ebN_cur+2];
                    //(nN_1-nN_0) x (nN_2-nN_0)
                    n_parent[0] = (nN_ebN_parent_1_x[1]-nN_ebN_parent_0_x[1])*(nN_ebN_parent_2_x[2]-nN_ebN_parent_0_x[2])
                      -(nN_ebN_parent_1_x[2]-nN_ebN_parent_0_x[2])*(nN_ebN_parent_2_x[1]-nN_ebN_parent_0_x[1]);
                    
                    n_parent[1] = (nN_ebN_parent_1_x[2]-nN_ebN_parent_0_x[2])*(nN_ebN_parent_2_x[0]-nN_ebN_parent_0_x[0])
                      -(nN_ebN_parent_1_x[0]-nN_ebN_parent_0_x[0])*(nN_ebN_parent_2_x[2]-nN_ebN_parent_0_x[2]);
                    
                    n_parent[2] = (nN_ebN_parent_1_x[0]-nN_ebN_parent_0_x[0])*(nN_ebN_parent_2_x[1]-nN_ebN_parent_0_x[1])
                      -(nN_ebN_parent_1_x[1]-nN_ebN_parent_0_x[1])*(nN_ebN_parent_2_x[0]-nN_ebN_parent_0_x[0]);
 
                    tmp = sqrt(n_parent[0]*n_parent[0] + n_parent[1]*n_parent[1] + n_parent[2]*n_parent[2]);
                    assert(tmp > 0.0);
                    n_parent[0] /= tmp; n_parent[1] /= tmp; n_parent[2] /= tmp;
                
                    //looking for face on parent with unit normal in same direction as child face
                    //but don't enforce same outside inside convention on mesh orderings
                    tmp = n_parent[0]*n_child[0] + n_parent[1]*n_child[1] + n_parent[2]*n_child[2];
                               
                    if (fabs(sqrt(tmp*tmp) - 1.0) < 1.0e-8)
                      {
                        ebN_parent = ebN_cur;
                        break;
                      }
                    
                  }
                assert(ebN_parent >= 0);
                childMesh.elementBoundaryMaterialTypes[ebN] = parentMesh.elementBoundaryMaterialTypes[ebN_parent];
              }//3d
          }//nExteriorBoundaries > 1
        else 
          {
            //only 1 exterior boundary on parent so must be same as child's face
            int ebN_parent = parentMesh.elementBoundariesArray[eN_parent*parentMesh.nElementBoundaries_element + 
                                                               lastExteriorElementBoundaryOnParent];
            childMesh.elementBoundaryMaterialTypes[ebN] = parentMesh.elementBoundaryMaterialTypes[ebN_parent];
            
          }
 
      }//exterior element boundaries
    return failed;
  }
}
 
int readElements(std::istream& meshFile, Mesh& mesh)
{
  using namespace std;
  assert(meshFile);
 
  string word,elementType;
  //read in the mesh file. I just read in each token in the order it
  //appeards in the .3dm files. This will break if there is
  //non-whitespace trash in the file.
  meshFile>>word;
  if(word == "MESH1D")
    {
      elementType = "E2E";
      mesh.nNodes_element = 2;
      //cout<<"Reading 1D edge mesh"<<endl;
    }
  else if(word == "MESH2D")
    {
      elementType = "E3T";
      mesh.nNodes_element = 3;
      //cout<<"Reading 2D triangular mesh"<<endl;
    }
  else if (word == "MESH3D")
    {
      elementType = "E4T";
      mesh.nNodes_element = 4;
      //cout<<"Reading 3D tetrahedral mesh"<<endl;
    }
  else
    {
      cerr<<"Unrecognized mesh type"<<endl;
      return 1;
    }
  vector<vector<int> > elementNodesVector;
  vector<int> nodes(mesh.nNodes_element);
  vector<int> materialTypes;
  int type;
  meshFile>>word;
  while(word == elementType)
    {
      meshFile>>word; //discard element numbering
      for(int nN=0;nN<mesh.nNodes_element;nN++)
        {
          meshFile>>nodes[nN];
          nodes[nN]-=1;
        }
      elementNodesVector.push_back(nodes);
      meshFile>>type;
      materialTypes.push_back(type);
      meshFile>>word;
    }
  //cout<<"Number of  elements = "<<elementNodesVector.size()<<endl;
  mesh.nElements_global = elementNodesVector.size();
  mesh.elementNodesArray = new int[mesh.nElements_global*mesh.nNodes_element];
  mesh.elementMaterialTypes = new int[mesh.nElements_global];
  for (int eN=0;eN<mesh.nElements_global;eN++)
    {
      mesh.elementMaterialTypes[eN] = materialTypes[eN];
      for (int nN=0;nN<mesh.nNodes_element;nN++)
        mesh.elementNodesArray[eN*mesh.nNodes_element+nN]=elementNodesVector[eN][nN];
    }
  return 1;
}
 
int writeNodes(std::ostream& meshFile, const Mesh& mesh)
{
  using namespace std;
  for(int nN=0;nN<mesh.nNodes_global;nN++)
    meshFile<<"ND"<<setw(7)<<nN+1
            <<scientific<<setprecision(8)<<setw(16)<<mesh.nodeArray[3*nN+0]
            <<scientific<<setprecision(8)<<setw(16)<<mesh.nodeArray[3*nN+1]
            <<scientific<<setprecision(8)<<setw(16)<<mesh.nodeArray[3*nN+2]
            <<endl;
  return 0;
}
 
int writeElements(std::ostream& meshFile, const Mesh& mesh)
{
  //std::cout<<"printing mesh "<<mesh.nElements_global<<"\t"<<mesh.nNodes_element<<std::endl;
  meshFile<<"Try to write something"<<std::endl;
  using namespace std;
  string elementType;
  int width;
  if(mesh.nNodes_element == 2)
    {
      width=7;
      elementType = "E2E";
      meshFile<<"MESH1D"<<endl;
    }
  else if(mesh.nNodes_element == 3)
    {
      width=7;
      elementType = "E3T";
      meshFile<<"MESH2D"<<endl;
    }
  else if (mesh.nNodes_element == 4)
    {
      width=8;
      elementType = "E4T";
      meshFile<<"MESH3D"<<endl;
    }
  else
    {
      cerr<<"Unknown element type"<<endl;
      return 1;
    }
  for (int eN=0;eN<mesh.nElements_global;eN++)
    {
      meshFile<<elementType;
      meshFile<<setw(width)<<eN+1;
      for (int nN=0;nN<mesh.nNodes_element;nN++)
        meshFile<<setw(width)<<(mesh.elementNodesArray[eN*mesh.nNodes_element+nN]+1);
      //mwftodo decide if have convention about material types and base 0
      meshFile<<setw(width)<<mesh.elementMaterialTypes[eN]+1;
      meshFile<<endl;
    }
  return 0;
}
 
 
int setFromTriangleElements(triangulateio* trimesh, Mesh& mesh, int base)
{
  int failed = 0;
  assert(trimesh); assert(trimesh->trianglelist);
  //assume mesh hasn't been allocated
  assert(mesh.nElements_global == 0);
 
  mesh.nNodes_element = 3;
  mesh.nElements_global = trimesh->numberoftriangles;
  mesh.elementNodesArray = new int[mesh.nElements_global*mesh.nNodes_element];
  mesh.elementMaterialTypes = new int[mesh.nElements_global];
  //copy material types now, copy first attribute of triangle
  if (trimesh->triangleattributelist)
    {
      for (int eN = 0; eN < trimesh->numberoftriangles; eN++)
        mesh.elementMaterialTypes[eN] = trimesh->triangleattributelist[eN*trimesh->numberoftriangleattributes+0]; 
    }
  else
    memset(mesh.elementMaterialTypes,DEFAULT_ELEMENT_MATERIAL,mesh.nElements_global*sizeof(int));
 
  for (int eN = 0; eN < mesh.nElements_global; eN++)
    {
      //copy onl vertices even if triangle stores "nonlinear" (6pt) triangle 
      for (int ebN = 0; ebN < mesh.nNodes_element; ebN++)
        {
          mesh.elementNodesArray[eN*mesh.nNodes_element+ebN] = 
            trimesh->trianglelist[eN*trimesh->numberofcorners+ebN]-base;
        }    
     
    }
 
  return failed;
}
 
int setFromTriangleNodes(triangulateio* trimesh, Mesh& mesh, int base)
{
  int failed = 0;
  assert(trimesh); assert(trimesh->pointlist);
  //assume mesh hasn't been allocated