output.cpp

Go to the documentation of this file.

#include "output.h"
#include <fstream>
#include <stdlib.h> // remove dir
#include <algorithm>
#include "input.h"
#include <iostream>
#include "structures.h"
#include "vectorFunctions.h"
#include "insertShape.h" // getFamilyNum()
#include "mathFunctions.h"
#include "generatingPoints.h"
#include "computationalGeometry.h" // Rotation matrix functions
#include <cmath>
#include <sys/stat.h> // Checking if output directory exists DIR_exists(dir)
#include <iomanip> // std::setprecision()
#include <sys/stat.h> // mkdir system call 
#include "readInputFunctions.h" // error check for file open checkIfOpen()

//NOTE: do not use std::endl for new lines when writing to files. This will flush the output buffer. Use '\n'

/* void writeOutput() ************************************************************************/
void writeOutput(char* outputFolder, std::vector<Poly> &acceptedPoly, std::vector<IntPoints> &intPts, std::vector<Point> &triplePoints, struct Stats &pstats, std::vector<unsigned int> &finalFractures, std::vector<Shape> &shapeFamilies) {
    
    std::string output = outputFolder;
    
    // Define Output Files:
    std::string permOutputFile = output + "/perm.dat";
    std::string aperture = output + "/aperture.dat";
    std::string intersectionFolder = output + "/intersections";
    std::string radiiFolder = output + "/radii/"; 

    // Create intersections folder
    //makeDIR(intersectionFolder.c_str()); 
    // Write intersection files (must be first file written, rotates polys to x-y plane)
    writeIntersectionFiles(finalFractures, acceptedPoly, intPts, triplePoints, intersectionFolder, pstats);
    // Write polys.inp
    writePolysInp(finalFractures, acceptedPoly, output);
    // Write params.txt
    writeParamsFile(finalFractures, acceptedPoly, shapeFamilies, pstats, triplePoints, output);
    // Write aperture file
    writeApertureFile(finalFractures, acceptedPoly, output);
    // Write permability file
    writePermFile(finalFractures, acceptedPoly, output);
    // Write radii file
    writeRadiiFile(finalFractures, acceptedPoly, output);
    // Write rejection stats file
    writeRejectionStats(pstats, output);
    // Write families to output Files
    writeShapeFams(shapeFamilies, output);
    // Write fracture translations file
    writeFractureTranslations(finalFractures, acceptedPoly, output);
    // Write fracture connectivity (edge graph) file
    writeConnectivity(finalFractures, acceptedPoly, intPts, output);
    // Write rotation data
    writeRotationData(acceptedPoly, finalFractures, shapeFamilies, output);
    // Write rejects per fracture insertion attempt data
    writeRejectsPerAttempt(pstats, output);
    // Write all accepted radii, for Nataliia
    writeFinalPolyRadii(finalFractures, acceptedPoly, output);

    if (outputAcceptedRadiiPerFamily) {
        std::cout << "Writing Accepted Radii Files Per Family\n";
        // Creates radii files per family, before isolated fracture removal.
        int size = shapeFamilies.size();
        for (int i = 0; i < size; i++) {
           writeAllAcceptedRadii_OfFamily(i, acceptedPoly, radiiFolder);
        }

        if (userRectanglesOnOff) {
            // Fractures are marked -2 for user rects
            writeAllAcceptedRadii_OfFamily(-2, acceptedPoly, radiiFolder);
        }

        if (userEllipsesOnOff) {
            // Fractures are marked -1 for user ellipses
            writeAllAcceptedRadii_OfFamily(-1, acceptedPoly, radiiFolder);
        }
    }


    if (outputFinalRadiiPerFamily) {
        
        std::cout << "Writing Final Radii Files Per Family\n";    
    
        int size = shapeFamilies.size();
        for (int i = 0; i < size; i++) {
           writeFinalRadii_OfFamily(finalFractures, i, acceptedPoly, radiiFolder);
        }

        if (userRectanglesOnOff) { 
            writeFinalRadii_OfFamily(finalFractures, -2, acceptedPoly, radiiFolder);
        }

        if (userEllipsesOnOff) {
            writeFinalRadii_OfFamily(finalFractures, -1, acceptedPoly, radiiFolder);
        }
    }

    // If triple intersections are on, write triple intersection points file
    if (tripleIntersections) {
        std::cout << "Writing Triple Intersection Points File\n";
        writeTriplePts(triplePoints, finalFractures, acceptedPoly, intPts, output);
    }

} // End writeOutput()


/*=================================   OUTPUT.CPP FUNCTIONS   ================================*/
/*===========================================================================================*/


/* void writePoints() ************************************************************************/
void writePoints(std::ostream &output, std::vector<Point> &points, int start, unsigned int &count){
    int n=points.size();
    for (int i = start; i < n; i++) {
        output << std::setprecision(13) << count << " " << points[i].x 
               << " " << points[i].y << " " << points[i].z<< "\n";
        count++;
    }
}

/* finishWritingIntFile() ********************************************************************/
void finishWritingIntFile(std::ostream &fractIntFile, int fract1, int numPoints, int numIntersections, 
     std::vector<unsigned int> &intStart, std::vector<unsigned int> &intersectingFractures){
    
    unsigned int count=1;
    int idx = 0;
    
    // Lines
    for (int i = 0; i<numPoints-numIntersections; i++){
        if(intStart[idx] == count+1){
            count++;
            idx++;
        }
        fractIntFile<< i+1 << " " <<fract1<< " line " <<count<< " " <<count+1<< "\n";
        count++;
    }    
    
    fractIntFile<< "2 1 1\n";
    fractIntFile<< "a_b, integer\n";
    fractIntFile<< "b_a, integer\n";
    idx=0;
    for (int i=0; i<numPoints; i++){
        if(intStart[idx] == (unsigned) i+1){
            idx++;
        }
    fractIntFile<<i+1<< " " <<fract1<< " " <<intersectingFractures[idx]<< "\n";
    }
    
    // Move to header postion
    fractIntFile.seekp(0);
    fractIntFile<<numPoints<< " " <<numPoints - numIntersections<< " " <<2<< " " << "0 0";
}



/* bool DIR_exists() *************************************************************************/
bool DIR_exists(const char *path){ 
    struct stat sb;

    if ((stat(path, &sb) == 0 && S_ISDIR(sb.st_mode))) {
        return 1;
    }
return 0;
}


/* adjustIntFractIds() **********************************************************************/
void adjustIntFractIDs(std::vector<unsigned int> &finalFractures, std::vector<Poly> &allPolys, std::vector<IntPoints> &intPts){
    //go through all final fractures
    for (unsigned int i = 0; i < finalFractures.size(); i++){
        //go through each final fractures intersections
        for (unsigned int j = 0; j < allPolys[finalFractures[i]].intersectionIndex.size(); j++){
            //change matching fracture numbers to new number (order of finalFractures list) for output
            unsigned int intIdx = allPolys[finalFractures[i]].intersectionIndex[j];
            if (intPts[intIdx].fract1 ==  finalFractures[i]){
                intPts[intIdx].fract1 =  -( (long int) i+1);
            }
            else if (intPts[intIdx].fract2 == finalFractures[i]){
                intPts[intIdx].fract2 = -( (long int) i+1);
            }
        }
    }
}


/* writeIntersectionFiles() ******************************************************************/
void writeIntersectionFiles(std::vector<unsigned int> &finalFractures, std::vector<Poly> &acceptedPoly, std::vector<IntPoints> &intPts, std::vector<Point> &triplePoints, std::string intersectionFolder, struct Stats &pstats) { 
    
    Point tempPoint1, tempPoint2; // Keeps track of current un-rotated points we are working with

    std::cout << "\nWriting Intersection Files\n";    

    adjustIntFractIDs(finalFractures,acceptedPoly, intPts);
 
    std::ofstream fractIntFile; 

    // Go through finalFractures. Rotate poly, intersections, and triple intersection points 
    // to XY plane. Discretize and write to file
    for (unsigned int i = 0; i < finalFractures.size(); i++){
        
        // Starting positions for each intersection. Lets us know how to make the line connections
        std::vector<unsigned int> intStart;         
        unsigned int count = 1;

        // Counter for intersection header  number of points
        unsigned int numIntPts = 0; 
        // Used in wiring which nodes belong to which two fractures in finishWritingOutput()
        std::vector<unsigned int> intersectingFractures; 
    
        // Make new intersection file in intersections folder
        std::string file = intersectionFolder + "/intersections_" + std::to_string(i+1) + ".inp";
        fractIntFile.open(file.c_str(), std::ofstream::out | std::ofstream::trunc);

        checkIfOpen(fractIntFile, file);

        // Buffer first line to leave space to write header  
        fractIntFile<< "                                                               \n";    
                                    
        // Go through each final fracture's intersections and write to output
        unsigned int size = acceptedPoly[finalFractures[i]].intersectionIndex.size();
        for (unsigned int j = 0; j < size; j++){        

            // tempTripPts holds rotated triple points for an intersection. Triple pts must be rotated 3 different 
            // ways so we cannot change the original data 
            std::vector<Point> tempTripPts; 

            // Used to measure current length before rotation (used to calculate number of points 
            // to discretize) based on original intersection. This fixes any precision errors we 
            // when calculating length after rotation, both rotations for the same intersection will
            // always have the same step size and same number of discretized points.
            double curLength = 0;
        
            unsigned int polyIntIdx = acceptedPoly[finalFractures[i]].intersectionIndex[j];
                        
            // Similarly to above, the intersection must be rotated two different ways, 
            // one for each intersecting poly. We can't change the original data so we must use temp data    
            IntPoints tempIntersection = polyAndIntersection_RotationToXY(intPts[polyIntIdx], 
                                         acceptedPoly[finalFractures[i]], triplePoints, tempTripPts);
            // poly and intersection now rotated
        
            int triplePtsSize = tempTripPts.size();
                
            // fracture 1 is i
            // fracture 2 is the other intersecting fracture
            unsigned int fract2;
            if (intPts[polyIntIdx].fract1 == -(i+1)){
                fract2 = -intPts[polyIntIdx].fract2;
                intersectingFractures.push_back(fract2);    
            }
            else {
                fract2 = -intPts[polyIntIdx].fract1;
                intersectingFractures.push_back(fract2);
            }
            // If triple points exist on intersection, discretize from endpoint to closest triple point,
            // from triple to next triple point, and finally to other end point
            if (triplePtsSize != 0) {
                
                // Keep track of number of triple points which will be in the 
                // DFN (this is after isolated fracture removal)
                // NOTE: This will need to be divided by six to get correct value. 
                // Division by six was determined through testing.
                pstats.tripleNodeCount += triplePtsSize;
                 
                // Order the triple points by distances to know to discretize from point to next closest point
                double *distances = new double[triplePtsSize]; 
            
                double pt1[3] = {tempIntersection.x1, tempIntersection.y1, tempIntersection.z1};
                tempPoint1.x = intPts[polyIntIdx].x1;
                tempPoint1.y = intPts[polyIntIdx].y1;
                tempPoint1.z = intPts[polyIntIdx].z1;
                // Create array of distances first end point to triple points
                for (int k = 0; k < triplePtsSize; k++) { //loop through triple points on  intersection i
                    double point[3] = {tempTripPts[k].x, tempTripPts[k].y, tempTripPts[k].z};//triple pt
                    distances[k] = euclideanDistance(pt1, point);//create array of distances 
                }
        
                // Order the indices of the distances array shortest to largest distance
                // this lets us know which point to discritize to next
                int *s = sortedIndex(distances, triplePtsSize);
            
                // Discretize from end point1 to first triple pt
                // pt1 already = enpoint1
                double pt2[3] = {tempTripPts[s[0]].x, tempTripPts[s[0]].y, tempTripPts[s[0]].z};
                tempPoint2 = triplePoints[intPts[polyIntIdx].triplePointsIdx[s[0]]];
                curLength = euclideanDistance(tempPoint1, tempPoint2);
                std::vector<Point> points = discretizeLineOfIntersection(pt1, pt2, curLength);
            
                // Write points to file
                numIntPts += points.size();
                writePoints(fractIntFile, points, 0, count);

                // If one trip pt, set up points to discretize from only triple pt to other end point
                if (triplePtsSize == 1){ 
                    pt1[0] = pt2[0];    pt1[1]=pt2[1];     pt1[2]=pt2[2];
                    pt2[0]=tempIntersection.x2;     pt2[1]=tempIntersection.y2;     pt2[2]=tempIntersection.z2; 

                    tempPoint1 = tempPoint2;
                    tempPoint2.x = intPts[polyIntIdx].x2;
                    tempPoint2.y = intPts[polyIntIdx].y2;
                    tempPoint2.z = intPts[polyIntIdx].z2;
                }
                else { // More than 1 triple point
                    for (int jj=0; jj < (triplePtsSize-1); jj++){
                    
                        pt1[0] = tempTripPts[s[jj]].x;   pt1[1]=tempTripPts[s[jj]].y; pt1[2]=tempTripPts[s[jj]].z;
                        pt2[0] = tempTripPts[s[jj+1]].x; pt2[1]=tempTripPts[s[jj+1]].y; pt2[2]=tempTripPts[s[jj+1]].z;

                        tempPoint1 = triplePoints[intPts[polyIntIdx].triplePointsIdx[s[jj]]];
                        tempPoint2 = triplePoints[intPts[polyIntIdx].triplePointsIdx[s[jj+1]]];
                        curLength = euclideanDistance(tempPoint1, tempPoint2);
                        points = discretizeLineOfIntersection(pt1, pt2, curLength);
                    
                        // Write points for first fracture to file, save second set of points to temp
                        numIntPts += points.size()-1;
                        writePoints(fractIntFile, points, 1, count);
                    }    
                    // Set up points to go from last triple point to last endpoint
                    pt1[0]=pt2[0];     pt1[1]=pt2[1];     pt1[2]=pt2[2];
                    pt2[0]=tempIntersection.x2;     pt2[1]=tempIntersection.y2;     pt2[2]=tempIntersection.z2;
                    
                    tempPoint1 = tempPoint2;
                    tempPoint2.x = intPts[polyIntIdx].x2;
                    tempPoint2.y = intPts[polyIntIdx].y2;
                    tempPoint2.z = intPts[polyIntIdx].z2;

                }
                curLength = euclideanDistance(tempPoint1, tempPoint2);
                points = discretizeLineOfIntersection(pt1, pt2, curLength);

                numIntPts += points.size()-1;
                writePoints(fractIntFile, points, 1, count);
            
                delete[] s; // Need to delete these manually. created with new[]
                delete[] distances;    
            }
            else { // No triple intersection points on intersection line
                double pt1[3] = {tempIntersection.x1, tempIntersection.y1, tempIntersection.z1};
                double pt2[3] = {tempIntersection.x2, tempIntersection.y2, tempIntersection.z2};
                tempPoint1.x = intPts[polyIntIdx].x1;
                tempPoint1.y = intPts[polyIntIdx].y1;
                tempPoint1.z = intPts[polyIntIdx].z1;

                tempPoint2.x = intPts[polyIntIdx].x2;
                tempPoint2.y = intPts[polyIntIdx].y2;
                tempPoint2.z = intPts[polyIntIdx].z2;
               
                curLength = euclideanDistance(tempPoint1, tempPoint2);
                std::vector<Point> points = discretizeLineOfIntersection(pt1, pt2, curLength);
                numIntPts += points.size();
                writePoints(fractIntFile, points, 0, count);
            }    
        
            intStart.push_back(count);
        } 
        // Done with fracture and intersections

        pstats.intersectionNodeCount += numIntPts;    

            
        // Write line connectivity and header
        finishWritingIntFile(fractIntFile, i+1, numIntPts, size, intStart, intersectingFractures);
        intersectingFractures.clear();
        intStart.clear();
        fractIntFile.close();
    }

    //Divide by 6 to remove the duplicate counts
    pstats.tripleNodeCount /= 6;
}

/* rotateFractures() **************************************************************************/
void rotateFractures(std::vector<unsigned int> &finalFractures, std::vector<Poly> &acceptedPoly) {
   
    for (unsigned int i = 0; i < finalFractures.size(); i++) {
            
        if (acceptedPoly[finalFractures[i]].XYPlane == true) {
            continue; // Go to next interation of loop
        }    

        acceptedPoly[finalFractures[i]].XYPlane = true;

        double normalB[3] = {0,0,1};
        applyRotation3D(acceptedPoly[finalFractures[i]], normalB); 
    }
}

/* writePolysInp() ****************************************************************************/
void writePolysInp_old(std::vector<unsigned int> &finalFractures, std::vector<Poly> &acceptedPoly, std::string &output){

    std::ofstream polyOutput;
    std::string polyOutputFile = output+"/polys.inp";
    polyOutput.open(polyOutputFile.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(polyOutput, polyOutputFile);   
    std::cout << "Writing " << polyOutputFile << "\n";

    unsigned long long int vertexCount = 0;

    //HEADER
    for (unsigned int j = 0; j < finalFractures.size(); j++){ // Count vertices
        vertexCount += acceptedPoly[finalFractures[j]].numberOfNodes;
    }

    polyOutput << vertexCount << " " << vertexCount - finalFractures.size() << " 0" << " 0" << " 0" << "\n";
    
    int count = 1;
    int polyCount = finalFractures.size();
    for (int j = 0; j<polyCount; j++){    

        // Write vertices    
        for (int i = 0; i < acceptedPoly[finalFractures[j]].numberOfNodes; i++) {
                int idx = i*3;        
                polyOutput << std::setprecision(15) << count << " " 
                           << acceptedPoly[finalFractures[j]].vertices[idx] << " " 
                           << acceptedPoly[finalFractures[j]].vertices[idx+1] << " " 
                           << acceptedPoly[finalFractures[j]].vertices[idx+2] << "\n"; 
                count++;
            }
    }    
        
    // Write line connectivity
    count = 1;// Counter for node numbers
    int count2 = 1;// Counter for lines written
    int j;
    for (j = 0; j<polyCount; j++){        
        for (int i = 0; i < acceptedPoly[finalFractures[j]].numberOfNodes-1; i++) {
            polyOutput << count2<< " " << j+1 << " line " << count << " " <<count+1<< "\n"; 
            count++;
            count2++;
        }
    count++;
    }
        
    polyOutput.close(); // Done with polygons inp file
}

/* writePolysInp() ****************************************************************************/
void writePolysInp(std::vector<unsigned int> &finalFractures, std::vector<Poly> &acceptedPoly, std::string &output) {

    std::ofstream polyOutput;
    std::cout << "Writing poly inp files\n";

    int polyCount = finalFractures.size();
  
    for (int j = 0; j < polyCount; j++) {    
        std::string polyOutputFile = output + "/polys/poly_" + std::to_string(j+1) + ".inp";
        polyOutput.open(polyOutputFile.c_str(), std::ofstream::out | std::ofstream::trunc);

        // Write header
        polyOutput << acceptedPoly[finalFractures[j]].numberOfNodes << " " 
                   << acceptedPoly[finalFractures[j]].numberOfNodes-1 << " 0" << " 0" << " 0 " << "\n";

        // Write vertices   
        int numberOfNodes = acceptedPoly[finalFractures[j]].numberOfNodes; 
        for (int i = 0; i < numberOfNodes; i++) {
            int idx = i*3;        
            polyOutput << std::setprecision(15) << i+1 << " " 
                       << acceptedPoly[finalFractures[j]].vertices[idx] << " " 
                       << acceptedPoly[finalFractures[j]].vertices[idx+1] << " " 
                       << acceptedPoly[finalFractures[j]].vertices[idx+2] << "\n"; 
        }

        // Write line connectivity
        for (int i = 1; i < numberOfNodes; i++) {
            polyOutput << i << " " << j+1  << " line " << i  << " " << i+1 << "\n";
        }
        polyOutput.close();
    }    
        
        
}

/* writeParamsFile() **************************************************************************/
void writeParamsFile(std::vector<unsigned int> &finalFractures, std::vector<Poly> &acceptedPoly, std::vector<Shape> &shapeFamilies, Stats &pstats, std::vector<Point> &triplePoints, std::string &output) {

    std::ofstream params;
    std::string paramsOutputFile = output+"/params.txt";
    params.open(paramsOutputFile.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(params, paramsOutputFile);
    std::cout << "Writing " << paramsOutputFile << "\n";
    params << finalFractures.size() << "\n";
    params << h << "\n";

    params << visualizationMode << "\n"; // Production mode
    params << pstats.intersectionNodeCount/2 - pstats.tripleNodeCount << "\n";
    
    params.close();
}


/* writeApertureFile() ************************************************************************/
void writeApertureFile(std::vector<unsigned int> &finalFractures, std::vector<Poly> &acceptedPoly, std::string &output) {
    std::string file = output+"/aperture.dat";
    std::ofstream ap;
    ap.open(file.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(ap, file);
    std::cout << "Writing aperture.dat\n";
    ap<< "aperture.dat" << "\n";
   
    int size = finalFractures.size();
    for (int i = 0; i < size; i++) {
        ap << -(7+i) << " 0 0 " << std::setprecision(10) << acceptedPoly[finalFractures[i]].aperture << "\n";
    }
    
    ap.close();
}

/* writePermFile() ****************************************************************************/
void writePermFile(std::vector<unsigned int> &finalFractures, std::vector<Poly> &acceptedPoly, std::string &output) {
    std::string file = output+"/perm.dat"; 
    std::ofstream perm;
    perm.open(file.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(perm, file);
    std::cout << "Writing perm.dat\n";
    perm<< "permeabiliy" << "\n";
    
    int size = finalFractures.size();
    for (int i = 0; i < size; i++) {
       perm << -(7+i) << " 0 0 " <<  std::setprecision(10) << acceptedPoly[finalFractures[i]].permeability << " " <<acceptedPoly[finalFractures[i]].permeability<< " " <<acceptedPoly[finalFractures[i]].permeability<< "\n";
                                            
       }
    
    perm.close();
}


/* writeRadiiFile() ***************************************************************************/
void writeRadiiFile(std::vector<unsigned int> &finalFractures, std::vector<Poly> &acceptedPoly, std::string &output) {

    std::cout << "Writing Radii File (radii.dat)\n";

    std::string file = output+"/radii.dat";
    std::ofstream radii;
    radii.open(file.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(radii, file);
    radii << "Format: xRadius yRadius Family# Removed (-1 = userRectangle, 0 = userEllipse, > 0 is family in order of famProb)\n";

    unsigned int finalFractLimit = finalFractures.size() - 1;
    unsigned int size = acceptedPoly.size();
    unsigned int curFinalIdx = 0;

    if (finalFractures.size() <= 0) return;

    for (unsigned int i = 0; i < size; i++) {
        radii <<  std::setprecision(8) << acceptedPoly[i].xradius << " " << acceptedPoly[i].yradius
                   << " " << acceptedPoly[i].familyNum+1;

        // If poly is not in finalFractures list
        // mark that is was removed 
        if (i != finalFractures[curFinalIdx]) {
           radii << " R\n";
        }
        else {
            if (curFinalIdx < finalFractLimit) {
                curFinalIdx++;
            }
            radii << "\n";
        }
    
    }
   radii.close();
}


/* writeFractureTranslations() ****************************************************************/
void writeFractureTranslations(std::vector<unsigned int> &finalFractures, std::vector<Poly> &acceptedPoly, std::string &output) {

    std::cout << "Writing Fracture Translations File (translations.dat)\n";
    std::string filePath = output + "/translations.dat";
    std::ofstream file;
    file.open(filePath.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(file, filePath);
    file << "Format: x y z  (R = removed from domain due to fracture isolation)\n";

    unsigned int finalFractLimit = finalFractures.size() - 1;
    unsigned int size = acceptedPoly.size();
    unsigned int curFinalIdx = 0;

    if (finalFractures.size() <= 0) return;

    for (unsigned int i = 0; i < size; i++) {
        file << std::setprecision(10) << acceptedPoly[i].translation[0] << " " 
             << acceptedPoly[i].translation[1] << " " << acceptedPoly[i].translation[2];

        // If poly is not in finalFractures list
        // mark that is was removed 
        if (i != finalFractures[curFinalIdx]) {
           file << " R\n";
        }
        else {
            if (curFinalIdx < finalFractLimit) {
                curFinalIdx++;
            }
            file << "\n";
        }
    }
   file.close();
}

/* writeFinalPolyRadii() **********************************************************************/
void writeFinalPolyRadii(std::vector<unsigned int> &finalFractures, std::vector<Poly> &acceptedPoly, std::string &output){
    std::string file = output + "/radii_Final.dat";
    std::ofstream radiiFinal;
    radiiFinal.open(file.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(radiiFinal, file);
    radiiFinal << "Fracture Radii List After Isolated Fracture and Cluster Removal\n";
    radiiFinal << "Format: xRadius yRadius Family# (-1 = userRectangle, 0 = userEllipse, > 0 is family in order of famProb)\n";

    int size =   finalFractures.size();
    for (int i = 0; i < size; i++) {
        radiiFinal << acceptedPoly[finalFractures[i]].xradius << " " 
                   << acceptedPoly[finalFractures[i]].yradius << " " 
                   << acceptedPoly[finalFractures[i]].familyNum + 1 << "\n";
    }
    radiiFinal.close();
}


/* writeAllAcceptedRadii() ********************************************************************/
void writeAllAcceptedRadii(std::vector<Poly> &acceptedPoly, std::string &output){
    std::string file = output + "/radii_AllAccepted.dat";
    std::ofstream radiiAcpt;
    radiiAcpt.open(file.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(radiiAcpt, file);
    radiiAcpt << "Fracture Radii List Before Isolated Fracture and Cluster Removal\n";
    radiiAcpt << "Format: xRadius yRadius Distrubution# (-1 = userRectangle, 0 = userEllipse, > 0 is family in order of famProb)\n";

    int size = acceptedPoly.size();
    for (int i = 0; i < size; i++) {
        radiiAcpt << acceptedPoly[i].xradius << " " << acceptedPoly[i].yradius << " " <<acceptedPoly[i].familyNum+1 << "\n";
    }
    radiiAcpt.close();
}


/* writeAllAcceptedRadii_OfFamily() ***********************************************************/
void writeAllAcceptedRadii_OfFamily(int familyNum, std::vector<Poly> &acceptedPoly, std::string &output){
    std::string fileName = output + "/radii_AllAccepted_Fam_" + std::to_string(familyNum + 1) + ".dat";
    std::ofstream file;
    file.open(fileName.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(file, fileName);
    file << "Fracture Radii List Before Isolated Fracture and Cluster Removal (Family " << familyNum+1 << ")\n";
    file << "Format: xRadius yRadius Distrubution# (-1 = userRectangle, 0 = userEllipse, > 0 is family in order of famProb)\n";

    int size = acceptedPoly.size();
    for (int i = 0; i < size; i++) {
        if (acceptedPoly[i].familyNum == familyNum) {
            file << acceptedPoly[i].xradius << " " << acceptedPoly[i].yradius << " " <<acceptedPoly[i].familyNum+1 << "\n";
        }
    }
    file.close();
}

/* writeAllAcceptedRadii_OfFamily() ***********************************************************/
void writeFinalRadii_OfFamily(std::vector<unsigned int> &finalFractures, int familyNum, std::vector<Poly> &acceptedPoly, std::string &output){
    std::string fileName = output+"/radii_Final_Fam_" + std::to_string(familyNum + 1) + ".dat";
    std::ofstream file;
    file.open(fileName.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(file, fileName);
    file << "Fracture Radii List After Isolated Fracture and Cluster Removal (Family " << familyNum+1 << ")\n";
    file << "Format: xRadius yRadius Distrubution# (-1 = userRectangle, 0 = userEllipse, > 0 is family in order of famProb)\n";

    int size = finalFractures.size();
    for (int i = 0; i < size; i++) {
        if (acceptedPoly[finalFractures[i]].familyNum == familyNum) {
            file << acceptedPoly[finalFractures[i]].xradius << " " << acceptedPoly[finalFractures[i]].yradius << " " <<acceptedPoly[finalFractures[i]].familyNum+1 << "\n";
        }
    }
    file.close();
}


/* writeAllAcceptedRadii_OfFamily() ***********************************************************/
void writeTriplePts(std::vector<Point> &triplePoints, std::vector<unsigned int> &finalFractures, std::vector<Poly> &acceptedPoly, std::vector<IntPoints> &intPts, std::string &output) {
 

    std::string fileName = output + "/triple_points.dat";
    std::ofstream file;
    file.open(fileName.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(file, fileName);
   
    // Save triple point indices from final list of fractures to temp array
    // There will be duplicates  which need to be removed
    // before writing to file
    std::vector<unsigned int> triplePtsList;

    int size = finalFractures.size();
    for (int i = 0; i < size; i++) {
        int intersectCount = acceptedPoly[finalFractures[i]].intersectionIndex.size();
        for (int j = 0; j < intersectCount; j++) {                                        
            int tripleSize = intPts[acceptedPoly[finalFractures[i]].intersectionIndex[j]].triplePointsIdx.size();
            for (int k = 0; k < tripleSize; k++) {
                int triplePtIdx = intPts[acceptedPoly[finalFractures[i]].intersectionIndex[j]].triplePointsIdx[k];
                triplePtsList.push_back(triplePtIdx);
            }
        }
    }
 
    size = triplePtsList.size();

    if (size > 0) {

        //remove duplicates
        std::sort(triplePtsList.begin(), triplePtsList.end());
        std::vector<unsigned int> finalPts;
        unsigned int  prevPt = triplePtsList[0];
        finalPts.push_back(prevPt);
        for (int i = 1; i < size; i++) {
            unsigned int curPt = triplePtsList[i];
            if (curPt != prevPt) {
                finalPts.push_back(curPt);
            }
            prevPt = curPt;
        }
        
        triplePtsList.clear();


        size = finalPts.size();
        for (int i = 0; i < size; i++) {
        file << std::setprecision(17) 
             << triplePoints[finalPts[i]].x  
             << " " << triplePoints[finalPts[i]].y 
             << " " << triplePoints[finalPts[i]].z 
             << "\n";
        }
    }
    file.close();
}


/* writeRejectionStats() **********************************************************************/
void writeRejectionStats(Stats &pstats, std::string &output) {
    std::cout << "Writing Rejection Statistics File (rejections.dat)\n";
    std::string fileName = output+"/rejections.dat";
    std::ofstream file;
    file.open(fileName.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(file, fileName);

    file << pstats.rejectionReasons.shortIntersection << " Short Intersections \n";
    file << pstats.rejectionReasons.closeToNode << " Close to Node\n";
    file << pstats.rejectionReasons.closeToEdge << " Close to Edge\n";
    file << pstats.rejectionReasons.closePointToEdge << " Vertice Close to Edge\n";
    file << pstats.rejectionReasons.outside << " Outside of Domain\n";
    file << pstats.rejectionReasons.triple << " Triple intersection Rejections\n";
    file << pstats.rejectionReasons.interCloseToInter << " Intersections Close to Other Intersections\n";
    
    file.close();
}


/* writeShapeFams() ***************************************************************************/
void writeShapeFams(std::vector<Shape> &shapeFamilies, std::string &output){
    
    double radToDeg = 180 / M_PI;

    std::cout << "Writing Family Definitions File (families.dat)\n";
    std::string fileName = output+"/families.dat";
    std::ofstream file;
    file.open(fileName.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(file, fileName);

    using namespace std;

    for(unsigned int i = 0; i < shapeFamilies.size(); i++) {

        //name(rect or ell) and number of family
        file << shapeType(shapeFamilies[i]) << " Family "
             << getFamilyNumber(i, shapeFamilies[i].shapeFamily) << ":\n";

        file << "Global Family " << i + 1 << "\n";

        // Print vertice number  
        if (shapeFamilies[i].shapeFamily == 0) {  // If ellipse family
            file << "Number of Vertices: " << shapeFamilies[i].numPoints << endl;
        }
        else {
            file << "Number of Vertices: 4" << endl;
        }

        // aspect ratio
        file << "Aspect Ratio: " << shapeFamilies[i].aspectRatio << endl;

        // p32 target 
        if (stopCondition == 1) {
            file << "P32 (Fracture Intensity) Target: "
                 << shapeFamilies[i].p32Target << endl;
        }

        // beta distribution, rotation around normal vector
        if (shapeFamilies[i].betaDistribution == 0) {
            file << "Beta Distribution (Rotation Around Normal Vector): [0, 2PI)" << endl;
        }
        else {
            file << "Beta (Rotation Around Normal Vector): "
                 << shapeFamilies[i].beta << " rad, " 
                 << shapeFamilies[i].beta * radToDeg << " deg" << endl;
        }

        // Theta (angle normal makes with z axis
        file << "Theta: " << shapeFamilies[i].theta << " rad, " 
             << shapeFamilies[i].theta * radToDeg << " deg" << endl;

        // Phi (angle the projection of normal onto x-y plane  makes with +x axis
        file << "Phi: " << shapeFamilies[i].phi << " rad, " 
             << shapeFamilies[i].phi * radToDeg << " deg" << endl;

        // kappa
        file << "Kappa: " << shapeFamilies[i].kappa << endl;


        // Print layer family belongs to    
        if (shapeFamilies[i].layer == 0) {
            file << "Layer: Entire domain" << endl;
        }
        else {
            int idx = (shapeFamilies[i].layer - 1) * 2;
            file << "Layer: " << shapeFamilies[i].layer << " {" << layers[idx]
                 << ", " << layers[idx+1]
                 << "}" << endl;
        }

        // Print distribution data        
        switch (shapeFamilies[i].distributionType) {
            case 1: // lognormal
                file << "Distribution: Lognormal\n";
                file << "Mean: " << shapeFamilies[i].mean << endl;
                file << "Standard Deviation: " << shapeFamilies[i].sd <<  endl;
                file << "Minimum Radius: " << shapeFamilies[i].logMin << "m" << endl;
                file << "Maximum Radius: " << shapeFamilies[i].logMax << "m" << endl;
                break;

            case 2: // power-law
                file << "Distribution: Truncated Power-Law\n";
                file << "Alpha: " << shapeFamilies[i].alpha << endl;
                file << "Minimum Radius: " << shapeFamilies[i].min << "m" << endl;
                file << "Maximum Radius: " << shapeFamilies[i].max << "m" << endl;
                break;

            case 3: // exponential
                file << "Distribution: Exponential\n";
                file << "Mean: " << shapeFamilies[i].expMean << endl;
                file << "Lambda: " << shapeFamilies[i].expLambda << endl;
                file << "Minimum Radius: " << shapeFamilies[i].expMin << "m" << endl;
                file << "Maximum Radius: " << shapeFamilies[i].expMax << "m" << endl;
                break;

            case 4: // constant
                file << "Distribution: Constant\n";
                file << "Radius: " << shapeFamilies[i].constRadi << "m" << endl;
        }
        file << "Family Insertion Probability: " << famProbOriginal[i] << "\n\n";
    }
    file.close();
}


/* makeDIR() **********************************************************************************/
void makeDIR(const char *dir) {

    // If dir already exists, remove it
    if (DIR_exists(dir)) {
        std::string tempStr = "rm -r ";
        tempStr += dir;
        if (system(tempStr.c_str())) {
            std::cout << "ERROR: Problem executing system " 
                      << "command: " << tempStr << "\n";
            exit(1);
        }
    }

    // If directory doesn't exist, create it
    if (!DIR_exists(dir)) {
        int dir_err = mkdir(dir, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
        if (-1 == dir_err)  {
        std::cout << "\nError creating directory " << dir << "\n";
        exit(1);
        }
    }
}


/* writeConnectivity() **********************************************************************************/
void writeConnectivity(std::vector<unsigned int> &finalFractures, std::vector<Poly> &acceptedPoly, std::vector<IntPoints> &intPts, std::string &output) {
    
    std::cout << "Writing Connectivity Data (connectivity.dat)\n";
    std::string fileName = output + "/connectivity.dat";
    std::ofstream file;
    file.open(fileName.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(file, fileName);

    // Format for connectivity: Line number = poly number, integers on line are intersecting fractures
    //  eg: 
    // 2 3 
    // 1 3
    // 1 2
    // Fracture 1 intersects frac 2 and 3, fract 2 intersects 1 and 3 etc

    for (unsigned int i = 0; i < finalFractures.size(); i++) {
    
        unsigned int size = acceptedPoly[finalFractures[i]].intersectionIndex.size();
        for (unsigned int j = 0; j < size; j++) {
            long int intIdx = acceptedPoly[finalFractures[i]].intersectionIndex[j];
            // Don't write fractures own fracture ID 

            // NOTE: At this stage in the program, the fracture ID's have been changed to 
            // negative. See 'adjustIntFractIds()' for more details.
            // Use '-' to make them positive again.
       
            if (finalFractures[i]+1 != -intPts[intIdx].fract1) {
                file << -intPts[intIdx].fract1 << " ";
            }
            else {
                file << -intPts[intIdx].fract2 << " ";
            }
        }
        file << "\n";
    }  
    file.close();
}


/* writeRotationData() ******************************************************************************/
void writeRotationData(std::vector<Poly> &acceptedPoly, std::vector<unsigned int> &finalFractures, std::vector<Shape> &shapeFamilies, std::string output) {

    std::ofstream file;
    std::string fileOutputFile = output+"/poly_info.dat";
    file.open(fileOutputFile.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(file, fileOutputFile);
    std::cout << "Writing Rotation Data File (poly_info.dat)\n";

        double maxDomainSize = domainSize[0];
        if (maxDomainSize < domainSize[1]) {
            maxDomainSize = domainSize[1];
        }
        if (maxDomainSize < domainSize[2]) {
            maxDomainSize = domainSize[2];
        } 
        
        maxDomainSize *= 10;        

    for (unsigned int i = 0; i < finalFractures.size(); i++ ) {

        // poly's normal is already normalized at this point
        double normal[3] = {acceptedPoly[finalFractures[i]].normal[0], acceptedPoly[finalFractures[i]].normal[1], acceptedPoly[finalFractures[i]].normal[2]};
        double e3[3] = {0,0,1};
        
        // Rotation angle in radians
        double theta = std::acos(dotProduct(normal, e3));
        // rad to deg
        theta = theta * (180/M_PI);
    
        // Rotation into xy plane
        double *v = crossProduct(e3, normal);
        if (!(std::abs(v[0]) < eps && std::abs(v[1]) < eps && std::abs(v[2]) < eps)){ //if not zero vector
            normalize(v);
        }
        
                
        double x0 = (-maxDomainSize*v[0]);
        double y0 = (-maxDomainSize*v[1]);
        double z0 = (-maxDomainSize*v[2]);
        
        double x1 = maxDomainSize*v[0];
        double y1 = maxDomainSize*v[1];
        double z1 = maxDomainSize*v[2];    
    
        
        // The last number is the shape family number. Throughout this program,
        // -1 and -2 are used to denote user defined ell. and rectangles.
        // -1 and -2 are not good numbers to use as material IDs in Lagrit.
        // If the family number is -1 or -2 we change these numbers to be
        // number of stochastic families + 1 and number of stochastic families + 2
        int famNum;
        if (acceptedPoly[finalFractures[i]].familyNum == -1) {
            famNum = shapeFamilies.size() + 1;
        }
        else if (acceptedPoly[finalFractures[i]].familyNum == -2) {
            famNum = shapeFamilies.size() + 2;
        }
        else {
            famNum = acceptedPoly[finalFractures[i]].familyNum + 1;
        }

        // Format: fracture#, x0, y0, z0, x1, y1, z1, family#
        file << (i+1) << " " << famNum << std::setprecision(15) << " " << theta << " " << x0 
               << " " << y0 << " " << z0 << " " << x1 << " " << y1 << " " << z1 << "\n";

        delete[] v;
    }    
    file.close();
}


/* writeRejectsPerAttempt()**************************************************************************/
void writeRejectsPerAttempt(Stats &pstats, std::string &output) {

    std::ofstream file;
    std::string fileOutputFile = output+"/rejectsPerAttempt.dat";
    file.open(fileOutputFile.c_str(), std::ofstream::out | std::ofstream::trunc);
    checkIfOpen(file, fileOutputFile);
    std::cout << "Writing Rotation Data File (rejectsPerAttempt.dat)\n";
 
    for (unsigned int i = 0; i < pstats.rejectsPerAttempt.size(); i++) {

        file << pstats.rejectsPerAttempt[i] << "\n";
    }

    file.close();
}