readInput.cpp

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <vector>
#include <math.h>
#include <cstdlib>
#include <string> 
#include "input.h"
#include "readInputFunctions.h"
#include "generatingPoints.h"

short stopCondition;

unsigned int nPoly;

double domainSize[3];

double h;

float radiiListIncrease;

bool disableFram;

bool visualizationMode;

bool tripleIntersections;

bool boundaryFaces[6];   

bool keepOnlyLargestCluster;

bool printRejectReasons;

bool outputFinalRadiiPerFamily;

bool outputAcceptedRadiiPerFamily;

bool *ebetaDistribution;

bool *rbetaDistribution;

bool insertUserRectanglesFirst;

bool forceLargeFractures;

unsigned int seed; 

float domainSizeIncrease[3];


int nFamRect;

int nFamEll;

float *famProb;

float *famProbOriginal; 

int *edistr;

float *easpect;

unsigned int *enumPoints;

bool  eAngleOption;

float *etheta;

float *ephi;

float *ebeta; 

float *ekappa;

float *eLogMean;

float *esd;

float *eExpMean;

float *rLogMin;

float *rLogMax;

float *rExpMin;

float *rExpMax; 

float *eLogMin;

float *eLogMax;

float *eExpMin;

float *eExpMax;

float *econst;

float *emin;

float *emax; 

float *ealpha;

float *e_p32Targets;

unsigned int *rdistr;

float *raspect;
 
bool rAngleOption;

float removeFracturesLessThan;

float *rtheta; 

float *rphi; 

float *rbeta;

float *rkappa;  

float *rLogMean; 

float *rsd;

float *rmin;

float *rmax;

float *ralpha;

float *r_p32Targets;

float *rExpMean;

float *rconst;

bool userEllipsesOnOff;

int nUserEll;

bool ueAngleOption;

float *ueRadii;

float *ueBeta;

float *ueaspect;

double *uetranslation;

double *uenormal; 

unsigned int *uenumPoints;

bool userRectanglesOnOff;

bool userRecByCoord;

bool userEllByCoord;

bool permOption;

bool outputAllRadii;

int nUserRect;

float *urRadii;

bool urAngleOption;

float *urBeta;

float *uraspect;

double *urtranslation; 

double *urnormal; 

unsigned int nRectByCoord;

unsigned int nEllByCoord;

unsigned int nEllNodes;

double *userRectCoordVertices;

double *userEllCoordVertices;

float meanAperture;

float stdAperture; 

int aperture;  

float apertureFromTransmissivity[2];

double constantAperture; 

double lengthCorrelatedAperture[2];

double constantPermeability;

int rejectsPerFracture; 

float *layers;

float *layerVol;

/*  Defines which domain, or layer, the family belings to.
    Layer 0 is the entire domain ('domainSize').
    Layers numbered > 0 coorespond to layers defined above (see 'Layers:').
    1 corresponts to the first layer listed, 2 is the next layer listed, etc*/
int *rLayer;

int *eLayer;

bool *p32Status;

bool ignoreBoundaryFaces;

int numOfLayers;

/******************************************************************/
void getInput(char* input, std::vector<Shape> &shapeFamily){

    std::string tempstring;
    char ch;
   
    std::ifstream inputFile;
    std::cout<<"DFN Generator Input File: "<< input <<"\n\n"; 

    // Open input file and initialize variables

    inputFile.open(input, std::ifstream::in);
    checkIfOpen(inputFile, tempstring.c_str());

    searchVar(inputFile, "stopCondition:");
    inputFile >> stopCondition;

    searchVar(inputFile, "printRejectReasons:");
    inputFile >> printRejectReasons;

    searchVar(inputFile, "domainSize:");
    inputFile >> ch >> domainSize[0] >> ch >> domainSize[1] >> ch >> domainSize[2];

    searchVar(inputFile, "numOfLayers:");
    inputFile >> numOfLayers;                 
                    

    if (numOfLayers > 0 ) {
        layers = new float[numOfLayers*2]; // Multiply by 2 for +z and -z for each layer
        layerVol = new float[numOfLayers];
        searchVar(inputFile, "layers:");

        std::cout << "Number of Layers: " << numOfLayers << "\n";
        for (int i = 0; i < numOfLayers; i++){
            int idx = i*2;
            inputFile >> ch >> layers[idx] >> ch >> layers[idx+1] >> ch;
            std::cout<<"    Layer "<<i+1<<"{-z,+z}: {"<< layers[idx] << ","<<layers[idx+1] << "}, Volume: ";
            layerVol[i] = domainSize[0]*domainSize[1]* (std::abs(layers[idx+1]-layers[idx]));
            std::cout << layerVol[i] << "m^3\n";
        }
        std::cout<<"\n";
    }

    searchVar(inputFile, "h:");
    inputFile >> h;

    searchVar(inputFile, "disableFram:");
    inputFile >> disableFram;

    if (disableFram == true) {
        std::cout<<"\nFRAM IS DISABLED\n";
    }

    searchVar(inputFile, "tripleIntersections:");
    inputFile >> tripleIntersections;

    searchVar(inputFile, "forceLargeFractures:");
    inputFile >> forceLargeFractures;

    searchVar(inputFile, "visualizationMode:");
    inputFile >> visualizationMode;

    if (disableFram == true){
        visualizationMode = 1;
    }
    
    searchVar(inputFile, "outputAllRadii:");
    inputFile >> outputAllRadii;

    searchVar(inputFile, "outputFinalRadiiPerFamily:");
    inputFile >> outputFinalRadiiPerFamily;

    searchVar(inputFile, "outputAcceptedRadiiPerFamily:");
    inputFile >> outputAcceptedRadiiPerFamily;

    searchVar(inputFile, "seed:");
    inputFile >> seed;

    searchVar(inputFile, "domainSizeIncrease:");
    inputFile >> ch >> domainSizeIncrease[0] >> ch >> domainSizeIncrease[1] >> ch >> domainSizeIncrease[2];

    searchVar(inputFile, "keepOnlyLargestCluster:");
    inputFile >> keepOnlyLargestCluster;

    searchVar(inputFile, "ignoreBoundaryFaces:");
    inputFile >> ignoreBoundaryFaces;
    
    searchVar(inputFile, "boundaryFaces:");
    inputFile >>ch>>boundaryFaces[0]>>ch>>boundaryFaces[1]>>ch>>boundaryFaces[2]>>ch>>boundaryFaces[3]>>ch>>boundaryFaces[4]>>ch>>boundaryFaces[5];

    searchVar(inputFile, "rejectsPerFracture:");
    inputFile>>rejectsPerFracture;

    searchVar(inputFile, "nFamRect:");
    inputFile >> nFamRect;

    searchVar(inputFile, "nFamEll:");
    inputFile >> nFamEll;

    searchVar(inputFile, "removeFracturesLessThan:");
    inputFile >> removeFracturesLessThan;

    if (nFamEll > 0 || nFamRect >0) {
        searchVar(inputFile, "famProb:");
        famProb = new float[(nFamEll+nFamRect)];
        getInputAry(inputFile, famProb, (nFamEll+nFamRect));

        searchVar(inputFile, "famProb:");
        famProbOriginal = new float[(nFamEll+nFamRect)];
        getInputAry(inputFile, famProbOriginal, (nFamEll+nFamRect));
            
        searchVar(inputFile, "radiiListIncrease:");
        inputFile >> radiiListIncrease;

    }


    if (nFamEll > 0) {

        searchVar(inputFile, "ebetaDistribution:");
        ebetaDistribution = new bool[nFamEll];
        getInputAry(inputFile, ebetaDistribution, nFamEll);

        searchVar(inputFile, "eLayer:");
        eLayer = new int[nFamEll];
        getInputAry(inputFile, eLayer, nFamEll);

        searchVar(inputFile, "edistr:");    
        edistr = new int[nFamEll];
        getInputAry(inputFile, edistr, nFamEll);

        searchVar(inputFile, "easpect:");
        easpect = new float[nFamEll];
        getInputAry(inputFile, easpect, nFamEll);

        searchVar(inputFile, "enumPoints:");
        enumPoints = new unsigned int[nFamEll];
        getInputAry(inputFile,enumPoints, nFamEll);

        searchVar(inputFile, "eAngleOption:");
        inputFile >> eAngleOption;

        searchVar(inputFile, "etheta:");
        etheta = new float[nFamEll];
        getInputAry(inputFile, etheta, nFamEll);

        searchVar(inputFile, "ephi:");
        ephi = new float[nFamEll];
        getInputAry(inputFile, ephi, nFamEll);

        searchVar(inputFile, "ebeta:");
        ebeta = new float[nFamEll];
        getInputAry(inputFile, ebeta, nFamEll);

        searchVar(inputFile, "ekappa:");
        ekappa = new float[nFamEll];
        getInputAry(inputFile, ekappa, nFamEll);

        searchVar(inputFile, "eLogMean:");
        eLogMean = new float[nFamEll];
        getInputAry(inputFile, eLogMean, nFamEll);

        searchVar(inputFile, "esd:");
        esd = new float[nFamEll];
        getInputAry(inputFile, esd, nFamEll);

        searchVar(inputFile, "eExpMean:");
        eExpMean = new float[nFamEll];
        getInputAry(inputFile, eExpMean, nFamEll);
        
        searchVar(inputFile, "emin:");
        emin = new float[nFamEll];
        getInputAry(inputFile, emin, nFamEll);

        searchVar(inputFile, "emax:");
        emax = new float[nFamEll];
        getInputAry(inputFile, emax, nFamEll);

        searchVar(inputFile, "ealpha:");
        ealpha = new float[nFamEll];
        getInputAry(inputFile, ealpha, nFamEll);

        searchVar(inputFile, "econst:");
        econst = new float[nFamEll];
        getInputAry(inputFile, econst, nFamEll);
        
        searchVar(inputFile, "eLogMin:");
        eLogMin = new float[nFamEll];
        getInputAry(inputFile, eLogMin, nFamEll);
        
        searchVar(inputFile, "eLogMax:");
        eLogMax = new float[nFamEll];
        getInputAry(inputFile, eLogMax, nFamEll);
        
        searchVar(inputFile, "eExpMin:");
        eExpMin = new float[nFamEll];
        getInputAry(inputFile, eExpMin, nFamEll);
 
        searchVar(inputFile, "eExpMax:");
        eExpMax = new float[nFamEll];
        getInputAry(inputFile, eExpMax, nFamEll);
               
        
        if (stopCondition == 1) {
            // Get temp array for ellipse p32 targets 
            // Used to simplify initialization of shape family structures below
            searchVar(inputFile, "e_p32Targets:");
            e_p32Targets = new float[nFamEll];
            getInputAry(inputFile, e_p32Targets, nFamEll);
        }
    }

    
        // Distribution counters
        int shape1 = 0; //longnormal dist counter
        int shape2 = 0; //truncated power-law dist counter
        int shape3 = 0; //exponential dist counter
        int shape4 = 0; //constant dist counter
        int betaCount = 0;

        // Create shape strucutres from data gathered above
        for (int i = 0; i < nFamEll; i++){
        struct Shape newShapeFam;

        newShapeFam.shapeFamily = 0; // shapFam = 0 = ellipse, 1 = rect
        newShapeFam.distributionType = edistr[i];
        newShapeFam.numPoints = enumPoints[i];
        newShapeFam.aspectRatio = easpect[i];
        newShapeFam.theta = etheta[i];
        newShapeFam.phi = ephi[i];
        newShapeFam.kappa = ekappa[i];
        newShapeFam.angleOption = eAngleOption;
        newShapeFam.layer = eLayer[i];
        
        generateTheta(newShapeFam.thetaList, newShapeFam.aspectRatio, newShapeFam.numPoints);
        

        if (ebetaDistribution[i] == 1 ){ // If constant user defined beta option
            newShapeFam.betaDistribution = 1;
            newShapeFam.beta = ebeta[betaCount];
            betaCount++; 
        }
        else {
            newShapeFam.betaDistribution = 0;
        }

        // dist options:1 = lognormal, 2= truncated power-law, 3= exponential, 4=constant
        switch (edistr[i]) {
            case 1: // Lognormal
                newShapeFam.mean = eLogMean[shape1];
                newShapeFam.sd = esd[shape1];
                newShapeFam.logMin = eLogMin[shape1];
                newShapeFam.logMax = eLogMax[shape1];
                shape1++;
                break;
            case 2: // Truncated power-law
                newShapeFam.min = emin[shape2];
                newShapeFam.max = emax[shape2];
                newShapeFam.alpha = ealpha[shape2];
                shape2++;
                break;
            case 3:

                newShapeFam.expMean = eExpMean[shape3];
                newShapeFam.expLambda = 1/eExpMean[shape3];
                newShapeFam.expMin = eExpMin[shape3];
                newShapeFam.expMax = eExpMax[shape3];
                shape3++;
                break;
            case 4:
                newShapeFam.constRadi = econst[shape4];
                shape4++;
                break;
        }
        
        if (stopCondition == 1){    
            newShapeFam.p32Target = e_p32Targets[i];
        }
        
        // Save shape family to perminant array
        shapeFamily.push_back(newShapeFam);

    }
    
    if (nFamEll > 0 ) {
        // Can now delete/free the memory for these arrays
        delete[] ebetaDistribution;
        delete[] edistr;
        delete[] easpect;
        delete[] enumPoints;
        delete[] etheta;
        delete[] ephi;
        delete[] ebeta;
        delete[] ekappa;
        delete[] eLogMean;
        delete[] esd;
        delete[] eExpMean;
        delete[] emin;
        delete[] emax;
        delete[] ealpha;
        delete[] econst;
        delete[] eLayer;
        delete[] eLogMin;
        delete[] eLogMax;
        delete[] eExpMin;
        delete[] eExpMax;
        
        if (stopCondition == 1){ 
            delete[] e_p32Targets;
        }
    }


    if (nFamRect >0) {

        searchVar(inputFile, "rbetaDistribution:");
        rbetaDistribution = new bool[nFamRect];
        getInputAry(inputFile, rbetaDistribution, nFamRect); 
         
        searchVar(inputFile, "rdistr:");
        rdistr = new unsigned int[nFamRect];
        getInputAry(inputFile, rdistr, nFamRect);

        searchVar(inputFile, "rLayer:");
        rLayer = new int[nFamRect];
        getInputAry(inputFile, rLayer, nFamRect);

        searchVar(inputFile, "raspect:");
        raspect = new float[nFamRect];
        getInputAry(inputFile, raspect, nFamRect);

        searchVar(inputFile, "rAngleOption:");
        inputFile >> rAngleOption;

        searchVar(inputFile, "rtheta:");
        rtheta = new float[nFamRect];
        getInputAry(inputFile, rtheta, nFamRect);

        searchVar(inputFile, "rphi:");
        rphi = new float[nFamRect];
        getInputAry(inputFile, rphi, nFamRect);

        searchVar(inputFile, "rbeta:");
        rbeta = new float[nFamRect];
        getInputAry(inputFile, rbeta, nFamRect);

        searchVar(inputFile, "rkappa:");
        rkappa = new float[nFamRect];
        getInputAry(inputFile, rkappa, nFamRect);

        searchVar(inputFile, "rLogMean:");
        rLogMean = new float[nFamRect];
        getInputAry(inputFile, rLogMean, nFamRect);

        searchVar(inputFile, "rsd:");
        rsd = new float[nFamRect];
        getInputAry(inputFile, rsd, nFamRect);

        searchVar(inputFile, "rmin:");
        rmin = new float[nFamRect];
        getInputAry(inputFile, rmin, nFamRect);

        searchVar(inputFile, "rmax:");
        rmax = new float[nFamRect];
        getInputAry(inputFile, rmax, nFamRect);

        searchVar(inputFile, "ralpha:");
        ralpha = new float[nFamRect];
        getInputAry(inputFile, ralpha, nFamRect);

        searchVar(inputFile, "rExpMean:");
        rExpMean = new float[nFamRect];
        getInputAry(inputFile, rExpMean, nFamRect);
               
        searchVar(inputFile, "rconst:");
        rconst = new float[nFamRect];
        getInputAry(inputFile, rconst, nFamRect);

        searchVar(inputFile, "rLogMin:");
        rLogMin = new float[nFamRect];
        getInputAry(inputFile, rLogMin, nFamRect);

        searchVar(inputFile, "rLogMax:");
        rLogMax = new float[nFamRect];
        getInputAry(inputFile, rLogMax, nFamRect);
        
        searchVar(inputFile, "rExpMin:");
        rExpMin = new float[nFamRect];
        getInputAry(inputFile, rExpMin, nFamRect);
        
        searchVar(inputFile, "rExpMax:");
        rExpMax = new float[nFamRect];
        getInputAry(inputFile, rExpMax, nFamRect);

        if (stopCondition == 1) {    
            // Get temp array for rectangle p32 targets, 
            // Used to simplify initialization of shape family structures below
            searchVar(inputFile, "r_p32Targets:");
            r_p32Targets = new float[nFamRect];
            getInputAry(inputFile, r_p32Targets, nFamRect);
        }
    }


    // Set stop condition variables
    if (nFamRect > 0 || nFamEll > 0) {

        if (stopCondition == 0){ // npoly option
            searchVar(inputFile, "nPoly:");
            inputFile >> nPoly;
        }
        else if (stopCondition == 1) { // Stop program when p32 conditions per fracture family are met
            p32Status = new bool [nFamRect + nFamEll]; // Staus array for whether or not the family has reached its p32 requirement 

            for (int i = 0; i < (nFamRect + nFamEll); i++){
                // Zero the status flags
                p32Status[i] = 0;
            }        
        }
    }

         
        // Counters, used to place variable into correct array index
        // distribution counters
        shape1 = 0; // Longnormal dist counter
        shape2 = 0; // Truncated power-law dist counter
        shape3 = 0; // Exponential dist counter
        shape4 = 0; // Constant dist counter
        betaCount = 0;
    

    // Create shape strucutres from data gathered above
    for (int i = 0; i < nFamRect; i++) {
    
        struct Shape newShapeFam;
       
        newShapeFam.shapeFamily = 1; // shapFam = 0 = ellipse, 1 = rect
        newShapeFam.distributionType = rdistr[i];
        newShapeFam.numPoints = 4; // Rectangle
        newShapeFam.aspectRatio = raspect[i];
        newShapeFam.theta = rtheta[i];
        newShapeFam.phi = rphi[i];
        newShapeFam.kappa = rkappa[i];
        newShapeFam.angleOption = rAngleOption;
        newShapeFam.layer = rLayer[i];
    
        if (rbetaDistribution[i] == 1 ) { // If constant beta option
            newShapeFam.betaDistribution = 1;
            newShapeFam.beta = rbeta[betaCount];
            betaCount++; 

        }
        else {
            newShapeFam.betaDistribution = 0;
        }
        
        // dist options:1 = lognormal, 2= truncated power-law, 3= exponential, 4=constant
        switch (rdistr[i]) {
            case 1: // Lognormal
                newShapeFam.mean = rLogMean[shape1];
                newShapeFam.sd = rsd[shape1];
                newShapeFam.logMin = rLogMin[shape1];
                newShapeFam.logMax = rLogMax[shape1];
                shape1++;
                break;
            case 2: // Truncated power-law
                newShapeFam.min = rmin[shape2];
                newShapeFam.max = rmax[shape2];
                newShapeFam.alpha = ralpha[shape2];
                shape2++;
                break;
            case 3:
                newShapeFam.expMean = rExpMean[shape3];
                newShapeFam.expLambda = 1/rExpMean[shape3];
                newShapeFam.expMin = rExpMin[shape3];
                newShapeFam.expMax = rExpMax[shape3];
                shape3++;
                break;
            case 4:
                newShapeFam.constRadi = rconst[shape4];
                shape4++;
                break;
        }
        
        if (stopCondition == 1){    
            newShapeFam.p32Target = r_p32Targets[i];
        }    
        
        
        // Save family to perminant array
        shapeFamily.push_back(newShapeFam);

    }

    if (nFamRect > 0 ) {
        // Can now delete/free the memory for these arrays
        delete[] rbetaDistribution;
        delete[] rdistr;
        delete[] raspect;
        delete[] rtheta;
        delete[] rphi;
        delete[] rbeta;
        delete[] rkappa;
        delete[] rLogMean;
        delete[] rsd;
        delete[] rExpMean;
        delete[] rmin;
        delete[] rmax;
        delete[] ralpha;
        delete[] rconst;
        delete[] rLayer;
        delete[] rLogMin;
        delete[] rLogMax;
        delete[] rExpMin;
        delete[] rExpMax;
        
        if (stopCondition == 1) { 
            delete[] r_p32Targets;
        }
    }


    searchVar(inputFile, "userEllipsesOnOff:");
    inputFile >> userEllipsesOnOff;

    if (userEllipsesOnOff != 0) {

        searchVar(inputFile, "UserEll_Input_File_Path:");
        inputFile >> tempstring;
        std::ifstream uEllFile;
        uEllFile.open(tempstring.c_str(),std::ifstream::in);
        checkIfOpen(uEllFile, tempstring);
        std::cout << "User Defined Ellipses File: " << tempstring << std::endl;

        searchVar(uEllFile, "nUserEll:");
        uEllFile >> nUserEll;

        searchVar(uEllFile, "Radii:");
        ueRadii = new float[nUserEll];
        getElements(uEllFile, ueRadii, nUserEll);

        searchVar(uEllFile, "Aspect_Ratio:");
        ueaspect = new float[nUserEll];
        getElements(uEllFile, ueaspect, nUserEll);

        searchVar(uEllFile, "AngleOption:");
        uEllFile >> ueAngleOption;

        searchVar(uEllFile, "Beta:");
        ueBeta = new float[nUserEll];
        getElements(uEllFile, ueBeta, nUserEll);

        searchVar(uEllFile, "Translation:");
        uetranslation = new double[3*nUserEll];
        get2dAry(uEllFile, uetranslation, nUserEll);

        searchVar(uEllFile, "Normal:");
        uenormal = new double[3*nUserEll];
        get2dAry(uEllFile, uenormal, nUserEll);

        searchVar(uEllFile, "Number_of_Vertices:");
        uenumPoints = new unsigned int[nUserEll];
        getElements(uEllFile, uenumPoints, nUserEll);

        uEllFile.close();
    }


    searchVar(inputFile, "userRectanglesOnOff:");
    inputFile >> userRectanglesOnOff;


    if (userRectanglesOnOff != 0) {

        searchVar(inputFile, "UserRect_Input_File_Path:");
        
        inputFile >> tempstring;
        std::ifstream uRectFile;
        uRectFile.open(tempstring.c_str(),std::ifstream::in);
        checkIfOpen(uRectFile, tempstring);
        std::cout << "User Defined Rectangles File: " << tempstring << std::endl;

        searchVar(uRectFile,"nUserRect:");
        uRectFile >> nUserRect;

        searchVar(uRectFile, "Radii:");
        urRadii = new float[nUserRect];
        getElements(uRectFile, urRadii, nUserRect);


        searchVar(uRectFile, "AngleOption:");
        uRectFile >> urAngleOption;


        searchVar(uRectFile, "Beta:");
        urBeta = new float[nUserRect];
        getElements(uRectFile, urBeta, nUserRect);


        searchVar(uRectFile, "Aspect_Ratio:");
        uraspect = new float[nUserRect];
        getElements(uRectFile, uraspect, nUserRect);


        searchVar(uRectFile, "Translation:");
        urtranslation = new double[3*nUserRect];
        get2dAry(uRectFile, urtranslation, nUserRect);
        
        searchVar(uRectFile, "Normal:");
        urnormal = new double[3*nUserRect];
        get2dAry(uRectFile, urnormal, nUserRect);



        uRectFile.close();
    }
    
    searchVar(inputFile, "userEllByCoord:");
    inputFile >> userEllByCoord;

    searchVar(inputFile, "userRecByCoord:");
    inputFile >> userRecByCoord;



    if ((userRectanglesOnOff == 1 || userRecByCoord == 1) && (userEllipsesOnOff == 1 || userEllByCoord == 1)) {
        searchVar(inputFile, "insertUserRectanglesFirst:");
        inputFile >> insertUserRectanglesFirst;
    }
    else {
        insertUserRectanglesFirst = 0;
    }



    if (userEllByCoord != 0) { 
        searchVar(inputFile, "EllByCoord_Input_File_Path:");
        inputFile >> tempstring;
        std::ifstream file;
        file.open(tempstring.c_str(),std::ifstream::in);
        checkIfOpen(file, tempstring);
        std::cout << "User Defined Ellipses by Coordinates File: " << tempstring << std::endl;
        
        searchVar(file, "nEllipses:");
        file >> nEllByCoord;

        searchVar(file, "nNodes:");
        file >> nEllNodes;

        userEllCoordVertices = new double[3 * nEllNodes * nEllByCoord];
        searchVar(file, "Coordinates:");

        getCords(file, userEllCoordVertices, nEllByCoord, nEllNodes);
        
        file.close();   
    }


    if (userRecByCoord != 0) {
        searchVar(inputFile, "RectByCood_Input_File_Path:");
        inputFile >> tempstring;
        std::ifstream uCoordFile;
        uCoordFile.open(tempstring.c_str(),std::ifstream::in);
        checkIfOpen(uCoordFile, tempstring);
        std::cout << "User Defined Rectangles by Coordinates File: " << tempstring << std::endl;

        searchVar(uCoordFile, "nRectangles:");
        uCoordFile >> nRectByCoord;

        searchVar(uCoordFile, "Coordinates:");
        //4 vertices, 12 elements x,y,z per rectangle: 
        userRectCoordVertices = new double[12 * nRectByCoord];
        getRectCoords(uCoordFile, userRectCoordVertices, nRectByCoord);

        uCoordFile.close();
    } 

    searchVar(inputFile, "aperture:");
    inputFile >> aperture;

    if (aperture == 1) {
    
        searchVar(inputFile, "meanAperture:");
        inputFile >> meanAperture;

        searchVar(inputFile, "stdAperture:");
        inputFile >> stdAperture;

    }
    else if (aperture == 2) {
        searchVar(inputFile,"apertureFromTransmissivity:");
        getInputAry(inputFile, apertureFromTransmissivity, 2);
    }
    else if (aperture == 3) {
        searchVar(inputFile, "constantAperture:");
        inputFile >> constantAperture;
    }  
    else if (aperture == 4) {
        searchVar(inputFile, "lengthCorrelatedAperture:");
        getInputAry(inputFile, lengthCorrelatedAperture, 2);
    }
    else {
        std::cerr<<"\nERROR: Aperture option not recognised\n";
        exit(1);
    }

    searchVar(inputFile, "permOption:");
    inputFile >> permOption;
    if (permOption != 0) {
        searchVar(inputFile, "constantPermeability:");
        inputFile >> constantPermeability;
    }

        // Error check on stopping parameter
    if (nFamEll + nFamRect == 0 && stopCondition != 0) { // If no stochastic shapes, use nPoly option with npoly = number of user polygons
        std::cout<<"\nWARNING: You have defined stopCondition = 1 (P32 program stopping condition) but have no stochastic shape families defined. Automatically setting stopCondition to 0 for use with user defined polygons and nPoly.\n\n";
        stopCondition = 0;
        
        if (userEllipsesOnOff == 0 && userRectanglesOnOff == 0 && userRecByCoord == 0 ) {
            std::cout << "ERROR: All polygon generating options are off or undefined, please check input file for errors.\n\n";
            exit(1);
        }
        
        int count = 0; // Count of user defined polygons
        if (userEllipsesOnOff == 1){
            count += nUserEll;
        }
        
        if (userRectanglesOnOff == 1){
            count += nUserRect;
        }
        
        if (userRecByCoord == 1) {
            count += nRectByCoord;
        }
        
        // Set nPoly to the amount of user defined polygons
        nPoly = count;
        
    }
    
    inputFile.close();

    // Convert angles to rad if necessary, all functions and code require radians
    for (unsigned int i = 0; i < shapeFamily.size(); i ++) {
        if (shapeFamily[i].angleOption == 1 ) { // Convert deg to rad
            double temp = M_PI/180;
            shapeFamily[i].beta *= temp;
            shapeFamily[i].theta *= temp;
            shapeFamily[i].phi *= temp;
            shapeFamily[i].angleOption = 0; // Angles now in radians
        }
    }


}

/****************************************************************************************/
// Depreciated Function
// Prints all input variables, useful for debugging.
// NOTE: Needs to be updated. There bave been changes to input variables
void printInputVars() {
    std::cout << "npoly = " << nPoly << std::endl;
    std::cout<< "tripleIntersections = " << tripleIntersections << std::endl;
    std::cout << "domainsize = {" << domainSize[0] << ", " << domainSize[1] << ", " << domainSize[2] << "}\n";
    std::cout << "h = " << h << std::endl;
    std::cout << "visualizationMode = " << visualizationMode << std::endl;
    std::cout << "seed = " << seed << std::endl;
    std::cout << "keepOnlyLargestCluster = "<<keepOnlyLargestCluster << std::endl;
    std::cout << "domainSizeIncrease = {" << domainSizeIncrease[0]<<","<<domainSizeIncrease[1]<<","<<domainSizeIncrease[2]<< "}\n";
    std::cout<< "boundaryFaces = {"<<boundaryFaces[0]<<","<<boundaryFaces[1]<<","<<boundaryFaces[2]<<","<<boundaryFaces[3]<<","<<boundaryFaces[4]<<","<<boundaryFaces[5]<<"}\n";
    std::cout<< "nFamRect = " <<nFamRect<<std::endl;
    std::cout << "nFamEll = " << nFamEll << std::endl;
    printAry(famProb, "famProb", (nFamEll+nFamRect));
    printAry(edistr, "edistr", nFamEll);
    printAry(enumPoints, "enumPoints", nFamEll);
    std::cout << "eAngleOption = " << eAngleOption << std::endl;
    printAry(easpect, "easpect", nFamEll);
    printAry(etheta, "etheta", nFamEll);
    printAry(ephi, "ephi", nFamEll);
    printAry(ebeta, "ebeta", nFamEll);
    printAry(ekappa, "ekappa", nFamEll);
    printAry(eLogMean, "eLogMean", nFamEll);
    printAry(esd, "esd", nFamEll);
    printAry(eExpMean, "eExpMean", nFamEll);
    printAry(econst, "econst", nFamEll);
    printAry(emin, "emin", nFamEll);
    printAry(emax, "emax", nFamEll);
    printAry(ealpha, "ealpha", nFamEll);
    printAry(rdistr, "rdistr", nFamRect);
    printAry(raspect, "raspect", nFamRect);
    std::cout << "rAngleOption = " << rAngleOption << std::endl;
    printAry(rtheta, "rtheta", nFamRect);
    printAry(rphi, "rphi", nFamRect);
    printAry(rbeta, "rbeta", nFamRect);
    printAry(rLogMean, "rLogMean", nFamRect);
    printAry(rsd, "rsd", nFamRect);
    printAry(rmin, "rmin", nFamRect);
    printAry(rmax, "rmax", nFamRect);
    printAry(ralpha, "ralpha", nFamRect);
    printAry(rkappa, "rkappa", nFamRect);
    printAry(rExpMean, "rExpMean", nFamRect);
    printAry(rconst, "rconst", nFamRect);
    std::cout << "userEllipsesOnOff = " << userEllipsesOnOff << std::endl;

    if (userEllipsesOnOff != 0) {
        std::cout<<"nUserEll = " << nUserEll << "\n";
        printAry(ueRadii, "ueRadii", nUserEll);  
        std::cout << "urAngleOption = " << ueAngleOption <<std::endl;
        printAry(ueBeta, "ueBeta", nUserEll);
        printAry(ueaspect, "ueaspect", nUserEll);
        print2dAry(uetranslation, "uetranslation", nUserEll);
        print2dAry(uenormal, "uenormal", nUserEll);
        printAry(uenumPoints, "uenumPoints", nUserEll);
    }

    std::cout << "userRectanglesOnOff = " << userRectanglesOnOff << std::endl;

    if (userRectanglesOnOff != 0) {

        std::cout<<"nUserRect = "<< nUserRect <<"\n";
        printAry(urRadii, "urRadii", nUserRect);
        std::cout<< "urAngleOption = " << urAngleOption << std::endl;
        printAry(urBeta, "urBeta", nUserRect);
        printAry(urRadii, "urRadii", nUserRect);
        printAry(uraspect, "uraspect", nUserRect);
        print2dAry(urtranslation, "urtranslation", nUserRect);
        print2dAry(urnormal, "urnormal", nUserRect);
    }
      
    std::cout << "userRecByCoord = " << userRecByCoord << std::endl;
    if (userRecByCoord != 0) {
        std::cout << "nRectByCoord = " << nRectByCoord << std::endl;
        printRectCoords(userRectCoordVertices, "userRectCoordVertices" , nRectByCoord);  
    }

    std::cout << "aperture option: " << aperture << std::endl;
    if (aperture == 1) {
        std::cout << "meanAperture = " << meanAperture << std::endl;
        std::cout << "stdAperture = " << stdAperture << std::endl;
      
    }
    else if (aperture == 2) {
        printAry(apertureFromTransmissivity, "apertureFromTransmissivity", 2);
    }
    else if (aperture == 3) {
        std::cout <<  "constantAperture = " << constantAperture << std::endl;
    }
    else if (aperture == 4) {
        printAry(lengthCorrelatedAperture, "lengthCorrelatedAperture", 2);
    }
 
      
    if (permOption == 0) {
        std::cout << "Permeability: Function of aperture\n";
    }
    else{
        std::cout << "ConstantPermeability = " << constantPermeability << std::endl;
    }  
}