importExport.cpp

#include "CResults.h"
#include "importExport.h"
#include "DSimpleMessage.h"

#include <armadillo>
#include <cassert>
#include <string_view>

using namespace arma;

void showErrorMessage(const string& message);
void showWarningMessage(const string& message);
void showInfoMessage(const string& message);
void showLineErrorMessage(
        const string& message,
        const uInt& lineIndex);


static void setAOA(CMesh *pmesh, double AOA)
{
    const double origin[3] = {0, 0, 0},
                 tip[3] = {0, 1, 0};
    setUInt nodeSelection, elemSelection;
    nodeSelection.clear(); elemSelection.clear();

    pmesh->rotate(AOA / pi * 180, origin, tip, nodeSelection, elemSelection);
}

static void setroll(CMesh *pmesh, double roll)
{
    const double origin[3] = {0, 0, 0},
                 tip[3] = {1, 0, 0};
    setUInt nodeSelection, elemSelection;
    nodeSelection.clear(); elemSelection.clear();

    pmesh->rotate(roll / pi * 180, origin, tip, nodeSelection, elemSelection);
}

static void setsweep(CMesh *pmesh,  double sweep)
{
    const double origin[3] = {0, 0, 0},
                 tip[3] = {0, 0, 1};
    setUInt nodeSelection, elemSelection;
    nodeSelection.clear(); elemSelection.clear();

    pmesh->rotate(sweep / pi * 180, origin, tip, nodeSelection, elemSelection);
}

static void setydist(CMesh *pmesh, double ydist)
{
    const double distance[3] = {0, ydist, 0};
    setUInt nodeSelection;
    nodeSelection.clear();

    pmesh->translate(distance, nodeSelection);
}

//#define DEBUG_FIND
int findright(
        const bool isright,
        vecUInt *found,
        const bool isspanwise,
        const uInt id,
        const vecUInt &elemids,
        const vecUInt &types,
        const vecVecUInt &nodeIDs,
        CMesh &mesh) {

    vecVecDbl cv;
    uInt elemid;

    if (isright) {
#ifdef DEBUG_FIND
        elemid = elemids[id]; //TODO: maybe i
        mesh.getElemCoords(elemid, cv);
        cerr << "appended " << id << " elemid:" << elemid <<
                " x: " << cv[0][0] << " y: " << cv[0][1] << " z: " << cv[0][2] << endl;
#endif
        found->emplace_back(id);
    }
    // get element where id ==0 and a[3] == b[2] && a[1] == b[4] maybe
    for (uInt i = 0; i < elemids.size(); i++) {
        if (i == id) continue;
#if 0
        if ((i % elemids.size() == 0 && (id + 1) % elemids.size() == 0) ||
                (id % elemids.size() == 0 && (i + 1) % elemids.size() == 0)) {
            cerr << " skip wrap " << endl;
            continue;
        }
#endif
        if (types[i] != 1) continue;
        //cerr << "i:" << i << " type:" << types[i] << endl;
        // get element where id ==0 and a[3] == b[2] && a[1] == b[4] maybe
        //if (id == 0) cerr << "node 0:" << nodeIDs[id][0] << " " << nodeIDs[id][1] << " " << nodeIDs[id][2] << " " << nodeIDs[id][3] << endl;
        //if (i == 1) cerr << "node 1:" << nodeIDs[i][0] << " " << nodeIDs[i][1] << " " << nodeIDs[i][2] << " " << nodeIDs[i][3] << endl;
        bool ismatch = false;
        if (isright) {
            if ((isspanwise && (nodeIDs[id][2] == nodeIDs[i][1] &&
                                 nodeIDs[id][3] == nodeIDs[i][0])) ||
                (!isspanwise && (nodeIDs[id][2] == nodeIDs[i][3] &&
                                nodeIDs[id][1] == nodeIDs[i][0]))) {
                ismatch = true;
            }
        } else {
            if ((isspanwise && (nodeIDs[id][3] == nodeIDs[i][2] &&
                                 nodeIDs[id][0] == nodeIDs[i][1])) ||
                (!isspanwise && (nodeIDs[id][0] == nodeIDs[i][3] &&
                                nodeIDs[id][1] == nodeIDs[i][2]))) {
                ismatch = true;
            }
        }
        if (ismatch) {
#ifdef DEBUG_FIND_END
            cerr << "FOUND isspanwise? " << isspanwise << " id= " << id << " i= " << i;
#endif
            elemid = elemids[i]; //TODO: maybe i
            mesh.getElemCoords(elemid, cv);
#ifdef DEBUG_FIND_END
            //cerr << "x: " << cv[0][0] << " y: " << cv[0][1] << endl;
            //cerr << "x: " << cv[1][0] << " y: " << cv[1][1] << endl;
            //cerr << "x: " << cv[2][0] << " y: " << cv[2][1] << endl;
            //cerr << "x: " << cv[3][0] << " y: " << cv[3][1] << endl;
#ifdef DEBUG_FIND_END
            cerr << " ELEM " << elemid << endl;
#endif
            //cerr << "id:" << id << " nodes: " <<nodeIDs[id][0] << " " << nodeIDs[id][1] << " " << nodeIDs[id][2] << " " << nodeIDs[id][3] << endl;
            //cerr << "i :" << i << " nodes: " <<nodeIDs[i][0] << " " << nodeIDs[i][1] << " " << nodeIDs[i][2] << " " << nodeIDs[i][3] << endl;
#endif
            //TODO hack - must prevent wrap correctly
            for (auto &it: *found) {
                if (it == i) {
#ifdef DEBUG_FIND_END
                    cerr << "skipping, already found: " << it << endl;
#endif
                    return 0;
                }
            }
            if (!isright) {
                continue;  // doesn't need to work in this direction
                for (auto &it: *found) {
#ifdef DEBUG_FIND_END
                    cerr << "it: " << it << " i: " << i << endl;
#endif
                    if (it == i) {
                        continue;
#ifdef DEBUG_FIND_END
                        cerr << "skipping, already found: " << it << endl;
#endif
                    }
                }
                found->emplace(found->begin(), i);
#ifdef DEBUG_FIND
                elemid = elemids[i]; //TODO: maybe i
                mesh.getElemCoords(elemid, cv);
                cerr << "prepended " << i << " elemid:" << elemid <<
                        " x: " << cv[0][0] << " y: " << cv[0][1] << " z: " << cv[0][2] << endl;
#endif
    }
            int ret = findright(isright, found, isspanwise, i, elemids, types, nodeIDs, mesh);
            return ret;
        }
    }
    return -1;
}

void findgrid(CMesh &mesh,
          vecUInt &spanwise,
          vector<vector<uInt>*> &chord,
          const vecUInt &elemids,
          const vecUInt &types,
          const vecVecUInt &nodeIDs) {

    // first, get first spanwise
    // then, for each spanwise, go chordwise
    //TODO still problem going left from origin.
    // for wing.bdf, element 81 should match 1218
    // 81: 1 2 337 80  1218: 1173 189 1 80  <---match 0=2  3=3 ???
#ifdef DEBUG_FIND
    cerr << "looking right" << endl;
#endif
    int fr = findright(true, &spanwise, true, 0, elemids, types, nodeIDs, mesh);
#ifdef DEBUG_FIND
    cerr << "looking left" << endl;
#endif
    fr = findright(false, &spanwise, true, 0, elemids, types, nodeIDs, mesh);
#ifdef DEBUG_FIND
    cerr << "spanwise: " << " size:" << spanwise.size() << endl;
#endif
    // then for each in 'spanwise', do another find for columns
    //cerr << "fr:" << fr << "spanwise:" << " size:" << spanwise.size() << "list: ";
    for (int i = 0; i < spanwise.size(); i++) {
        chord.push_back(new vector<uInt>);
#ifdef DEBUG_FIND
        cerr << "chord for " << spanwise[i] << " " << endl;;
        cerr << "looking down" << endl;
#endif
        int fr = findright(true, chord[i], false, spanwise[i], elemids, types, nodeIDs, mesh);
#if 0
#ifdef DEBUG_FIND
        cerr << "looking up" << endl;
#endif
        fr = findright(false, chord[i], true, spanwise[i], elemids, types, nodeIDs, mesh);
#ifdef DEBUG_FIND
        cerr << "spanwise: " << " size:" << chord[i]->size() << endl;
#endif
#endif
    }
}

void read_NASTRAN(const string meshfile, cube &XYZ, uInt &nnodes, uInt &nelems, int &m, int &n)
{
    CMesh mesh;
    mesh.clearData();

    const long shiftNodeIds = 0, shiftElemIds = 0;

    readNastranFile(&mesh, meshfile, shiftNodeIds, shiftElemIds);
    //cerr << mesh.getNumBodyQuad() << endl;
    // type=1 QUAD 2 TRIANGLE
    vecUInt elemids;
    vecUInt types;
    vecVecUInt nodeIDs;
    vecVecVecDbl coordinates;

    nnodes = mesh.getNumNode();
    nelems = mesh.getNumBodyQuad();

    mesh.getElemsData(elemids, types, nodeIDs);
    //cerr << "elemids (count" << elemids.size() << "):";
    //for (auto &eid: elemids) {
        //cerr << " " << eid;
    //}
    //cerr << endl;
    vecUInt spanwise;
    vector<vector<uInt>*> chord;
    findgrid(mesh, spanwise, chord, elemids, types, nodeIDs);
#if 0
    // don't do this here
    setAOA(&mesh, AOA);
    setroll(&mesh, roll);
    setsweep(&mesh, sweep);
    cerr << "not calling setydist!!!" << endl; // setydist(&mesh, ydist);
#endif

    //load_grid(blade_grid, importXYZ);
    //TODO warn/faill if chords don't have consistent size
    n = spanwise.size();
    int m2 = chord[0]->size();
    m = m2 / 2;
    XYZ.set_size(3, 2 * m + 2, n + 1);
    cerr << "NASTRAN parse found m=" << m << " n=" << n << endl;
    int ie = 0;
    vector<uInt> *el;
    XYZ.fill(0.0);
    //XYZ.fill(0.9999);
    // TODO add top right point [3] horizontally and [2] vertically for entire row
    // start view from bottom with trailing edge up
    // --------------------------
    // 0 3
    // 1 2
    int j;
    int abcd = 3;
    vecVecDbl cv;
    for (j = 0; j < chord.size(); j++) {
        el = chord[j];
        int i;
        abcd = 0;
        for (i = 0; i < (*el).size(); i++ ) {
        //for (uInt eid; *el) {
            //cerr << "ELEM index " << i << " ";
            uInt elemid = elemids[(*el)[i]]; //TODO: maybe i
            //cerr << "ELEM " << elemid;
            mesh.getElemCoords(elemid, cv);
            //cerr << "x: " << cv[0][0] << " y: " << cv[0][1] << " z: " << cv[0,2] << endl;
            //cerr << "x: " << cv[1][0] << " y: " << cv[1][1];
    //cerr << " i: " << i << " j: " << j <<" x: " << cv[2][0] << " y: " << cv[2][1] << " z: " << cv[2][2] << endl;
            //cerr << "x: " << cv[3][0] << " y: " << cv[3][1] << endl;
            XYZ(0, i, j) = cv[abcd][0];
            XYZ(1, i, j) = cv[abcd][1];
            XYZ(2, i, j) = cv[abcd][2];
            //if (i == 0) XYZ.print(cerr, "XYZ at 1=");
            ie++;
        }
        abcd = 1;
        XYZ(0, i, j) = cv[abcd][0];
        XYZ(1, i, j) = cv[abcd][1];
        XYZ(2, i, j) = cv[abcd][2];
        //cerr << "x: " << cv[3][0] << " y: " << cv[3][1] << " z: " << cv[3][2] << endl;
        ie++;
        //cerr << "last in chord: x: " << cv[3][0] << " y: " << cv[3][1] << endl;
    }
    // add right row, top right points
    el = chord[chord.size() - 1];
    int i = 0;
    abcd = 3;
    for (auto ii = (*el).begin(); ii != (*el).end(); ii++) {
        //cerr << "ELEM index " << *ii << " ";
        uInt elemid = elemids[*ii]; //TODO: maybe i
        //cerr << "ELEM " << elemid;
        mesh.getElemCoords(elemid, cv);
        //cerr << " i: " << i << " n:" << n << " x: " << cv[1][0] << " y: " << cv[1][1] << endl;
        XYZ(0, i, n) = cv[abcd][0];
        XYZ(1, i, n) = cv[abcd][1];
        XYZ(2, i, n) = cv[abcd][2];
        ie++;
        i++;
    }
    // add bottom right point
    el = chord[chord.size() - 1];
    uInt elemid = elemids[el->at(el->size() - 1)];
    mesh.getElemCoords(elemid, cv);
    abcd = 2;
    //cerr << "x: " << cv[2][0] << " y: " << cv[2][1] << endl;
    //cerr << "sizes:" << XYZ.n_cols << " " << XYZ.n_slices <<  " j: " << j << endl;
    //cerr << "m: " << m << " n:" << n << " x: " << cv[0][0] << " y: " << cv[0][1] << endl;
    XYZ(0, m2, n) = cv[abcd][0];
    XYZ(1, m2, n) = cv[abcd][1];
    XYZ(2, m2, n) = cv[abcd][2];
    ie++;
#if 0
    for (j = 0; j < n + 1; j++) {
        //cerr << "span: " << j << endl;
        for (i = 0; i < m2 + 2; i++) {
            cerr << XYZ(0, i, j) << " " <<
                    XYZ(1, i, j) << " " <<
                    XYZ(2, i, j) << " " << endl;
        }
        //cerr << endl;
    }
#endif
    //cerr << "last in last: x: " << cv[2][0] << " y: " << cv[2][1] << endl;

    //cerr << "ie:" << ie << " m+1 " <<  (m+1) << " n+1 " << (n+1) << endl;
    assert(ie == (m2+1) * (n+1));
}

//-----------------------------------------------------------------------/
// reads Nastran input file into mesh
//-----------------------------------------------------------------------/
void readNastranFile(
        CMesh* pMesh,
        const string& fileName,
        const long& shiftNodeIDs,
        const long& shiftElemIDs)
{
    // open file
    ifstream inputFile(fileName.c_str());
    if (!inputFile.is_open()){
        showErrorMessage("Cannot open NASTRAN file");
        return;
    }

    // read lines with stripped comments
    vecStr fileLines;
    string line;
    while (getLine(inputFile, line))
        fileLines.push_back(line.substr(0, line.find('$')));
    inputFile.close();

    // read nodes first
    bool nodefound = false;
    for (uInt i=0; i<fileLines.size(); i++){
        // check if this is node line
        if (convert2UpperCase(removeLeadingSpaces(fileLines[i])).substr(0,4) == "GRID"){
            nodefound = true;

            // declare node definition strings
            string sNodeID;
            vecStr coords(3);

            // check if free format: GRID, ID, , X, Y, Z, $ this is a comment
            if (fileLines[i].find(',') != string::npos){
                // split free format line
                vecStr lineParts = splitString(fileLines[i], ',');

                // check number of line parts
                if (lineParts.size() < 6){
                    showLineErrorMessage("Improper free format grid definition", i+1);
                    return;
                }
                sNodeID = lineParts[1];
                coords[0] = lineParts[3];
                coords[1] = lineParts[4];
                coords[2] = lineParts[5];
            }
            // check if small format (by searching of "*" in first field)
            else if (fileLines[i].substr(0,8).find('*') == string::npos){
                // check line size
                if (fileLines[i].size() <= 40){
                    showLineErrorMessage("Line too short for node definition", i+1);
                    return;
                }
                sNodeID = fileLines[i].substr(8,8);
                coords[0] = fileLines[i].substr(24,8);
                coords[1] = fileLines[i].substr(32,8);
                coords[2] = fileLines[i].substr(40,8);
            }
            // large format is what remains
            else{
                // check line size
                if (fileLines[i].size() <= 56){
                    showLineErrorMessage("Line too short for node definition", i+1);
                    return;
                }
                sNodeID = fileLines[i].substr(8,16);
                coords[0] = fileLines[i].substr(40,16);
                coords[1] = fileLines[i].substr(56,16);
                // skip to next line for additional Z coordinate
                if (fileLines.size() == ++i){
                    showLineErrorMessage("Premature end of file", i+1);
                    return;
                }
                if (fileLines[i].size() <= 8){
                    showLineErrorMessage("Line too short for node definition", i+1);
                    return;
                }
                coords[2] = fileLines[i].substr(8,16);
            }

            // extract data from strings
            //cerr << " sNodeID:" << sNodeID << " shiftNodeIDs:" << shiftNodeIDs << endl;
            uInt nodeID = atoi(sNodeID.c_str())+shiftNodeIDs;
            if (nodeID <= 0){
                showLineErrorMessage("Node ID should be greater than 0", i+1);
                cerr << "nodeID=" << nodeID << " sNodeID:" << sNodeID << " shiftNodeIDs:" << shiftNodeIDs << endl;
                return;
            }
            double coordinates[3];
            coordinates[0] = readNastranSci(coords[0]);
            coordinates[1] = readNastranSci(coords[1]);
            coordinates[2] = readNastranSci(coords[2]);

            // insert node into mesh
            int retCode = pMesh->insertNodeWithIDCoords(nodeID, coordinates);
            if (retCode == -1){
                showLineErrorMessage("Node ID already exists", i+1);
                return;
            }
            else if (retCode != 0){
                showLineErrorMessage("Unknown error reading node", i+1);
                return;
            }
        }
    }

    // check if something was read
    if (!nodefound)
        showWarningMessage("Zero nodes found");

    // read elements
    bool elemfound = false;
    for (uInt i=0; i<fileLines.size(); i++){
        // check if this is (supported) element line
        if (convert2UpperCase(removeLeadingSpaces(fileLines[i])).substr(0,6) == "CTRIA3" ||
            convert2UpperCase(removeLeadingSpaces(fileLines[i])).substr(0,6) == "CQUAD4" ){
            elemfound = true;

            // declare element definition strings
            string sElemID;
            vecStr sNodeIDs;

            // check if free format
            if (fileLines[i].find(',') != string::npos){
                // split free format line
                vecStr lineParts = splitString(fileLines[i], ',');

                // add data common to both element types
                if (lineParts.size() > 5){
                    sElemID = lineParts[1];
                    sNodeIDs.resize(3);
                    sNodeIDs[0] = lineParts[3];
                    sNodeIDs[1] = lineParts[4];
                    sNodeIDs[2] = lineParts[5];
                }

                // check if quad type and split line size adequate
                if (convert2UpperCase(removeLeadingSpaces(fileLines[i])).substr(0,6) == "CQUAD4"){
                    if (lineParts.size() > 6){
                        // add additional node
                        sNodeIDs.resize(4);
                        sNodeIDs[3] = lineParts[6];
                    }
                    else{
                        showLineErrorMessage("Improper free format element definition", i+1);
                        return;
                    }
                }
            }
            // check if small format (by searching of "*" in first field)
            else if (fileLines[i].substr(0,8).find('*') == string::npos){
                // check line size
                if (fileLines[i].size() <= 40){
                    showLineErrorMessage("Line too short for element definition", i+1);
                    return;
                }
                sElemID = fileLines[i].substr(8,8);
                sNodeIDs.resize(3);
                sNodeIDs[0] = fileLines[i].substr(24,8);
                sNodeIDs[1] = fileLines[i].substr(32,8);
                sNodeIDs[2] = fileLines[i].substr(40,8);

                // check if quad type and split line size adequate
                if (convert2UpperCase(removeLeadingSpaces(fileLines[i])).substr(0,6)=="CQUAD4" && fileLines[i].size()>48){
                    // add additional node
                    sNodeIDs.resize(4);
                    sNodeIDs[3] = fileLines[i].substr(48,8);
                }
            }
            // large format is what remains
            else{
                // check line size
                if (fileLines[i].size() <= 56){
                    showLineErrorMessage("Line too short for element definition", i+1);
                    return;
                }
                sElemID = fileLines[i].substr(8,16);
                sNodeIDs.resize(3);
                sNodeIDs[0] = fileLines[i].substr(40,16);
                sNodeIDs[1] = fileLines[i].substr(56,16);
                // skip to next line for additional node ID
                if (fileLines.size() == ++i){   // LINE INCREMENT HERE
                    showLineErrorMessage("Premature end of file", i+1);
                    return;
                }
                if (fileLines[i].size() <= 8){
                    showLineErrorMessage("Line too short for element definition", i+1);
                    return;
                }
                sNodeIDs[2] = fileLines[i].substr(8,16);
                if (convert2UpperCase(removeLeadingSpaces(fileLines[i])).substr(0,6)=="CQUAD4" && fileLines[i].size()>32){
                    sNodeIDs.resize(4);
                    sNodeIDs[3] = fileLines[i].substr(32,16);
                }
                else{
                    showLineErrorMessage("Line too short for element definition", i+1);
                    return;
                }
            }

            // convert data from strings
            vecUInt nodeIDs(sNodeIDs.size());
            for (uInt j=0; j<sNodeIDs.size(); j++)
                nodeIDs[j] = str2Int(sNodeIDs[j])+shiftNodeIDs;

            // insert element into mesh
            int retVal = pMesh->insertElemWithIDNodeIDs(str2Int(sElemID)+shiftElemIDs, nodeIDs);

            // process return value
            string message;
            if      (retVal == -1) message = "Duplicate element ID detected";
            else if (retVal == -2) message = "Element ID 0 detected";
            else if (retVal == -4) message = "Duplicate node ID(s) given in element definition";
            else if (retVal == -5) message = "Element refers to non-existing node ID(s)";
            else if (retVal == -6) message = "Zero surface element detected";
            if (retVal != 0){
                showLineErrorMessage(message, i+1);
                return;
            }
        }
    }

    // check if something was read
    if (!elemfound)
        showWarningMessage("Zero elements found");
}

//-----------------------------------------------------------------------/
// converts Nastran scientific format (missing "e") into double
// search is based on finding '-' or '+' sign without preoccuring character
//-----------------------------------------------------------------------/
double readNastranSci(const string& nastranNumber)
{
    // determine if scientific format
    bool scientificFormat = false;
    bool signExists = false;
    int exponentLocation = 0;
    for (int i=nastranNumber.size()-1; i>=0; i--){
        // check if sign exists
        if (nastranNumber[i]=='-' || nastranNumber[i]=='+')
            signExists = true;
        else if (signExists && nastranNumber[i]!=' '){
            // scientific format detected
            scientificFormat = true;
            exponentLocation = i+1;
            break;
        }
    }

    // read data into double
    if (scientificFormat){
        double coef = atof(nastranNumber.substr(0, exponentLocation).c_str());
        double expo = atof(nastranNumber.substr(   exponentLocation).c_str());
        return coef*pow(10, expo);
    }
    else
        return atof(nastranNumber.c_str());
}

//-----------------------------------------------------------------------/
// reads VTK file into mesh
//-----------------------------------------------------------------------/
void readVTKFile(
        CMesh* pMesh,
        const string& fileName,
        const long& shiftNodeIDs,
        const long& shiftElemIDs)
{
    // open file
    ifstream inputFile(fileName.c_str());
    if (!inputFile.is_open()){
        showErrorMessage("Cannot open VTK file");
        return;
    }

    // read lines with stripped comments
    vecStr fileLines;
    string line;
    while (getLine(inputFile, line))
        fileLines.push_back(line.substr(0, line.find("**")));
    inputFile.close();

    // read nodes
    uInt nodeCounter = 0;
    bool nodefound = false;
    for (uInt lineIndex=0; lineIndex<fileLines.size(); lineIndex++){
        vecStr vecStartLine = splitString(change2SingleSpace(fileLines[lineIndex]), ' ');
        if (vecStartLine.size()>1 && convert2UpperCase(vecStartLine[0])=="POINTS"){
            // node definition table found so read nodes
            uInt nodeNum = atoi(vecStartLine[1].c_str());
            if (checkEndOfFile(nodeNum, lineIndex, fileLines.size()))
                break;

            for (uInt nodeIndex=0; nodeIndex<nodeNum; nodeIndex++){
                lineIndex++;
                vecStr vecLine = splitString(change2SingleSpace(fileLines[lineIndex]), ' ');
                if (vecLine.size() != 3){
                    cerr << fileLines[lineIndex] << endl;

                    showLineErrorMessage(fileLines[lineIndex], lineIndex);
                    showLineErrorMessage("Incorrect node definition", lineIndex);
                    return;
                }

                uInt nodeID = nodeCounter+1+shiftNodeIDs;
                nodeCounter++;

                // read node coordinates
                double coordinates[3];
                for (int i=0; i<3; i++)
                    coordinates[i] = atof(vecLine[i].c_str());

                int retCode = pMesh->insertNodeWithIDCoords(nodeID, coordinates);
                if (retCode == -2){
                    showLineErrorMessage("Node ID already exists", lineIndex);
                    return;
                }
                else if (retCode == -3){
                    showLineErrorMessage("Node ID is equal 0", lineIndex);
                    return;
                }
                nodefound = true;
            }
        }
    }

    // check if something was read
    if (!nodefound)
        showWarningMessage("Zero nodes found");

    // read elements
    uInt elemCounter = 0;
    bool elemfound = false;
    for (uInt lineIndex=0; lineIndex<fileLines.size(); lineIndex++){
        vecStr vecStartLine = splitString(change2SingleSpace(fileLines[lineIndex]), ' ');
        if (vecStartLine.size()>1 && convert2UpperCase(vecStartLine[0])=="CELLS"){
            // element definition table found so read elements
            uInt cellNum = atoi(vecStartLine[1].c_str());
            if (checkEndOfFile(cellNum, lineIndex, fileLines.size()))
                break;

            for (uInt cellIndex=0; cellIndex<cellNum; cellIndex++){
                lineIndex++;
                vecStr vecLine = splitString(change2SingleSpace(fileLines[lineIndex]), ' ');
                // check element type
                uInt elemType;
                uInt elemId;
                vecUInt nodeIDs;
                if (vecLine.size()==5 && vecLine[0]=="4"){
                    // read linear quad plate elements
                    elemType = 1;

                    // read element ID
                    elemId = elemCounter+1+shiftElemIDs;
                    elemCounter++;
                    if (elemId <= 0){
                        showLineErrorMessage("Element ID should be greater than 0", lineIndex);
                        return;
                    }

                    // read node IDs
                    nodeIDs.resize(4);
                    for (uInt i=0; i<4; i++)
                        nodeIDs[i] = atoi(vecLine[i+1].c_str())+1+shiftNodeIDs;
                }
                else if (vecLine.size()==4 && vecLine[0]=="3"){
                    // read linear tria plate elements
                    elemType = 2;

                    // read element ID
                    elemId = elemCounter+1+shiftElemIDs;
                    elemCounter++;
                    if (elemId <= 0){
                        showLineErrorMessage("Element ID should be greater than 0", lineIndex);
                        return;
                    }

                    // read node IDs
                    nodeIDs.resize(3);
                    for (uInt i=0; i<3; i++)
                        nodeIDs[i] = atoi(vecLine[i+1].c_str())+1+shiftNodeIDs;
                }

                // insert element into mesh
                int retVal = pMesh->insertElemWithIDTypeNodeIDs(elemId, elemType, nodeIDs);

                // process return value
                string message = "Error in element definition";
                if      (retVal == -2) message = "Duplicate element ID detected";
                else if (retVal == -3) message = "Element ID 0 detected";
                else if (retVal == -4) message = "Duplicate node ID(s) given in element definition";
                else if (retVal == -5) message = "Zero surface element detected";
                else if (retVal == -6) message = "Element refers to non-existing node ID";
                if (retVal != 0){
                    showLineErrorMessage(message, lineIndex);
                    return;
                }

                elemfound = true;
            }
        }
    }

    // check if something was read
    if (!elemfound)
        showWarningMessage("Zero elements found");
}

//-----------------------------------------------------------------------/
// exports Nastran file
//-----------------------------------------------------------------------/
int exportNastranFile(
        const string& fileName,
        const vector<string>& comment,
        CMesh* pMesh,
        CResults* pResults,
        const bool& writeMesh,
        const bool& writeProperty,
        const bool& writeLoads,
        const string& modulus,
        const string& poisson,
        const string& thickness,
        const string& selectedData,
        const string& selectedCase)
{
    // open Nastran file for writing
    ofstream nastranFile(fileName.c_str());
    if (!nastranFile.is_open()){
        showErrorMessage("Cannot open NASTRAN file for writing");
        return -1;
    }

    nastranFile.precision(15);

    // write title
    nastranFile << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$" << endl;
    nastranFile << "$  NASTRAN (.bdf) mesh file generated by IT CAS for demonstration purposes" << endl;
    nastranFile << "$    illustrating a typical mesh to be imported by the PM software." << endl;
    nastranFile << "$  The PM user must eventually provide a similar mesh file representing some kind of wing, blade or foil." << endl;
    nastranFile << "$" << endl;
    for (auto &c: comment) nastranFile << "$  " << c << endl;
    nastranFile << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$" << endl;

    // write load
    if (writeLoads){
        // write step data
        nastranFile << "INIT MASTER(S)" << endl;
        nastranFile << "NASTRAN SYSTEM(319)=1" << endl;
        nastranFile << "ID Apame,Apame" << endl;
        nastranFile << "SOL SESTATIC" << endl;
        nastranFile << "TIME 10000" << endl;
        nastranFile << "CEND" << endl;
        nastranFile << "TITLE = Static Load from Apame 3D Panel Method" << endl;
        nastranFile << "ECHO = NONE" << endl;
        nastranFile << "DISPLACEMENT(PLOT) = ALL" << endl;
        nastranFile << "STRESS(PLOT) = ALL" << endl;
        nastranFile << "LOAD = 1" << endl;
        nastranFile << "BEGIN BULK" << endl;
        nastranFile << "$" << endl;

        // get index of selected data in results data structure
        int dataIndex = -1;
        for (uInt i=0; i<pResults->surfDatas.size(); i++){
            if (pResults->surfDatas[i].dataName == selectedData){
                dataIndex = i;
                break;
            }
        }
        if (dataIndex == -1){
            showErrorMessage("Selected field not found in results, contact support!");
            nastranFile.close();
            return -1;
        }

        // get index of case in results data structure
        int caseIndex = -1;
        for (uInt i=0; i<pResults->surfDatas[dataIndex].caseTitles.size(); i++){
            if (pResults->surfDatas[dataIndex].caseTitles[i] == selectedCase){
                caseIndex = i;
                break;
            }
        }
        if (caseIndex == -1){
            showErrorMessage("Selected case not found in results, contact support!");
            nastranFile.close();
            return -1;
        }

        // check data sizes
        if (pResults->surfDatas[dataIndex].data[caseIndex].size() != pResults->surfDatas[dataIndex].elemIDs.size()){
            showErrorMessage("Field data size not matching element number, contact support!");
            nastranFile.close();
            return -1;
        }

        // write distributed load data
        nastranFile << "*DLOAD, OP=NEW" << endl;
        for (uInt elemIndex=0; elemIndex<pResults->surfDatas[dataIndex].elemIDs.size(); elemIndex++)
            nastranFile << "PLOAD4,1," << pResults->surfDatas[dataIndex].elemIDs[elemIndex] << "," << -pResults->surfDatas[dataIndex].data[caseIndex][elemIndex] << endl;
        nastranFile << "$" << endl;
    }

    // write mesh
    if (writeMesh){
        // get node data
        vecUInt nodeIDs;
        vecVecDbl coordinates;
        pMesh->getNodeData(nodeIDs, coordinates);

        // write nodes
        if (!nodeIDs.empty()){
            for (uInt nodeIndex=0; nodeIndex<nodeIDs.size(); nodeIndex++)
                nastranFile << "GRID," << nodeIDs    [nodeIndex]    << ",0," <<
                                          coordinates[nodeIndex][0] << ","  <<
                                          coordinates[nodeIndex][1] << ","  <<
                                          coordinates[nodeIndex][2] << ",0" << endl;
            nastranFile << "$" << endl;
        }

        // get element data
        vecUInt elemIDs;
        vecUInt elemTypes;
        vecVecUInt elemNodes;
        pMesh->getElemsData(elemIDs, elemTypes, elemNodes);

        // write quad elements
        if (pMesh->getNumBodyQuad()+pMesh->getNumDummQuad() > 0){
            for (uInt elemIndex=0; elemIndex<elemIDs.size(); elemIndex++){
                if (elemTypes[elemIndex]==1 || elemTypes[elemIndex]==20)
                    nastranFile << "CQUAD4," << elemIDs  [elemIndex]    << ",1," <<
                                                elemNodes[elemIndex][0] << ","  <<
                                                elemNodes[elemIndex][1] << ","  <<
                                                elemNodes[elemIndex][2] << ","  <<
                                                elemNodes[elemIndex][3] << endl;
            }
            nastranFile << "$" << endl;
        }

        // write tria elements
        if (pMesh->getNumBodyTria()+pMesh->getNumDummTria() > 0){
            for (uInt elemIndex=0; elemIndex<elemIDs.size(); elemIndex++){
                if (elemTypes[elemIndex]==2 || elemTypes[elemIndex]==21){
                    nastranFile << "CTRIA3," << elemIDs  [elemIndex]    << ",1," <<
                                                elemNodes[elemIndex][0] << ","  <<
                                                elemNodes[elemIndex][1] << ","  <<
                                                elemNodes[elemIndex][2] << endl;
                }
            }
            nastranFile << "$" << endl;
        }
    }

    // write property data
    if (writeProperty){
        // write plate property
        nastranFile << "PSHELL,1,1," << thickness << ",1,1,0." << endl;
        nastranFile << "$" << endl;
        // write material data
        nastranFile << "MAT1,1," << modulus << "," << poisson << ",0.,0.,0." << endl;
        nastranFile << "$" << endl;
    }

    nastranFile << "ENDDATA" << endl;

    nastranFile.close();

    return 0;
}

auto read_file(string_view path) -> std::string {
    constexpr auto read_size = std::size_t{4096};
    auto stream = std::ifstream{path.data()};
    stream.exceptions(std::ios_base::badbit);

    auto out = std::string{};
    auto buf = std::string(read_size, '\0');
    while (stream.read(& buf[0], read_size)) {
        out.append(buf, 0, stream.gcount());
    }
    out.append(buf, 0, stream.gcount());
    return out;
}