правила создания библиотеки gcdfes, использующей библиотеку comsdk (4-5 этапы)

README.md

@@ -10,4 +10,11 @@
 В данной версии nla3d улучшен алгоритм вывода данных для VTK-файлов. Скорость вывода увеличена в разы, особенно хорошо это заметно на больших сетках.
 Для успешной компиляции необходимо наличие исходников библиотеки линейной алгебры eigen.
 В ходе компиляции будет большое количество warning'ов, это нормально и связано с eigen.
+4. С помощью CMAKE откомпилировать библиотеку comsdk (https://sa2systems.ru:881/com/comsdk).
+5. Создать проект библиотеки gcdfes, состоящий из четырёх файлов - gcdfes.cpp (https://sa2systems.ru:881/gcd/gcdfes/blob/master/dev/gcdfes.cpp), gcdfes.h (https://sa2systems.ru:881/gcd/gcdfes/blob/master/dev/gcdfes.h) и FEReaders.cpp/FEReaders.h (внутри папки libgcdfes).
+При создании библиотеки указать, какие решатели будем использовать - paralution (-DNLA3D_USE_PARALUTION), mkl (-DNLA3D_USE_MKL) или оба сразу.
+Пример использования решателей mkl внутри кода: math::PARDISO_equationSolver eqSolver = math::PARDISO_equationSolver();
+Пример использования решателей paralution внутри кода: math::PARALUTION_equationSolver eqSolver = math::PARALUTION_equationSolver();
+При создании библиотеки gcdfes для линковщика подключить следующие библиотеки в следующем порядке (для варианта, когда используются как решатели mkl, так и решатели paralution): -lcomsdk -lnla3d -lmath -lparalution -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl -lboost_system
+Пример проекта CodeBlocks (libgcdfes.cbp) и make-файла находятся в папке libgcdfes (для корректной работы нужно добавить два файла - gcdfes.cpp и gcdfes.h).
 

libgcdfes/FEReaders.cpp

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d
+
+//#define NOMINMAX
+
+#include "FEReaders.h"
+#include <unordered_map>
+#include <algorithm>
+
+using std::ifstream;
+using std::vector;
+using std::string;
+using std::map;
+using std::cout;
+using std::set;
+using std::endl;
+
+const string space_chars = " \t";
+
+
+void MeshData::clear() {
+    cellNumbers.clear();
+    cellNodes.clear();
+    cellDoubleData.clear();
+    cellIntData.clear();
+    cellStringData.clear();
+    nodesNumbers.clear();
+    nodesPos.clear();
+    loadBcs.clear();
+    fixBcs.clear();
+    mpcs.clear();
+    feComps.clear();
+}
+
+
+void MeshData::compressNumbers() {
+  // numbers should start from 1
+  if (nodesNumbers.size() == *std::max_element(nodesNumbers.begin(), nodesNumbers.end())) {
+    // nothing to do
+    return;
+  }
+  std::map<uint32, uint32> old2new;
+  uint32 nextNumber = 1;
+  for (uint32 i = 0; i < nodesNumbers.size(); i++) {
+    old2new[nodesNumbers[i]] = nextNumber;
+    nodesNumbers[i] = nextNumber;
+    nextNumber++;
+  }
+  assert(nextNumber - 1 == nodesNumbers.size());
+
+  // go through elements and change element's node numbers
+  for (auto e : cellNodes) {
+    for (uint16 i = 0; i < e.size(); i++) {
+      e[i] = old2new[e[i]];
+    }
+  }
+
+  // go through loadBCs
+  for (auto bc : loadBcs) {
+    bc.node = old2new[bc.node];
+  }
+
+  // go through fixBCs
+  for (auto bc : fixBcs) {
+    bc.node = old2new[bc.node];
+  }
+
+  // go through Mpcs
+  for (auto mpc : mpcs) {
+    for (auto term : mpc->eq) {
+      term.node = old2new[term.node];
+    }
+  }
+
+  // go through FE components
+  for (auto pair : feComps) {
+    auto comp = pair.second;
+    if (comp.type == FEComponent::NODES) {
+      for (auto n : comp.list) {
+        n = old2new[n];
+      }
+    }
+  }
+}
+
+
+std::vector<uint32> MeshData::getDegeneratedCells() {
+  std::vector<uint32> res;
+  for (size_t i = 0; i < cellNumbers.size(); i++) {
+    std::set<uint32> s(cellNodes[i].begin(), cellNodes[i].end());
+    if (s.size() != cellNodes[i].size()) {
+      res.push_back(i);
+    }
+  }
+  return res;
+}
+
+
+std::vector<uint32> MeshData::getCellsByAttribute(std::string atr_name, uint32 atr_val) {
+  std::vector<uint32> res;
+  for (size_t i = 0; i < cellNumbers.size(); i++) {
+    if (cellIntData[atr_name][i] == atr_val) {
+      res.push_back(i);
+    }
+  }
+  return res;
+}
+
+
+string&& strim(string&& str) {
+  size_t st = str.find_first_not_of(space_chars, 0);
+  str.erase(0, st);
+  size_t en = str.find_last_not_of(space_chars);
+  str.erase(en+1);
+  return std::move(str);
+}
+
+
+string& strim(string& str) {
+  size_t st = str.find_first_not_of(space_chars, 0);
+  str.erase(0, st);
+  size_t en = str.find_last_not_of(space_chars);
+  str.erase(en+1);
+  return str;
+}
+
+
+string& stolower(string& str) {
+  std::transform(str.begin(), str.end(), str.begin(), ::tolower);
+  return str;
+}
+
+
+string& stoupper(string& str) {
+  std::transform(str.begin(), str.end(), str.begin(), ::toupper);
+  return str;
+}
+
+
+char sfirstNotBlank(const string& str) {
+  size_t st = str.find_first_not_of(space_chars, 0);
+  return str[st];
+}
+
+
+std::vector<std::string> ssplit(const std::string& line, const std::vector<int>& widths, bool strict) {
+  int ind = 0;
+  vector<string> vv;
+  for (int w : widths) {
+    if (line.length() - ind < w) {
+      // no room for another field
+      if (strict) {
+        LOG(FATAL) << "ssplit: incomplete line: \"" << line << "\"";
+      } else {
+        if (line.length() - ind != 0) {
+          LOG(FATAL) << "ssplit: incomplete field: \"" << line << "\"";
+        }
+        return std::move(vv);
+      }
+    }
+    vv.push_back(line.substr(ind, w));
+    ind += w;
+  }
+  return vv;
+}
+
+
+bool iequals(const string& a, const string& b) {
+    unsigned int sz = a.size();
+    if (b.size() != sz)
+        return false;
+    for (unsigned int i = 0; i < sz; ++i)
+        if (tolower(a[i]) != tolower(b[i]))
+            return false;
+    return true;
+}
+
+// taken from
+// http://stackoverflow.com/questions/6089231/getting-std-ifstream-to-handle-lf-cr-and-crlf/6089413#6089413
+std::istream& getLine(std::istream& is, std::string& t) {
+  t.clear();
+
+  // The characters in the stream are read one-by-one using a std::streambuf.
+  // That is faster than reading them one-by-one using the std::istream.
+  // Code that uses streambuf this way must be guarded by a sentry object.
+  // The sentry object performs various tasks,
+  // such as thread synchronization and updating the stream state.
+
+  std::istream::sentry se(is, true);
+  std::streambuf* sb = is.rdbuf();
+
+  for (;;) {
+    int c = sb->sbumpc();
+    switch (c) {
+      case '\n':
+        return is;
+      case '\r':
+        if (sb->sgetc() == '\n') {
+          sb->sbumpc();
+        }
+        return is;
+      case EOF:
+        // Also handle the case when the last line has no line ending
+        if (t.empty()) {
+          is.setstate(std::ios::eofbit);
+        }
+        return is;
+      default:
+        t += (char)c;
+    }
+  }
+}
+
+
+int Tokenizer::tokenize(const string& line) {
+  tokens.clear();
+  dtokens.clear();
+
+  string _delimiters(delimiters.begin(), delimiters.end());
+
+  int found = 0;
+  size_t st = 0;
+
+  while (1) {
+    size_t ind = line.find_first_of(_delimiters, st);
+    found++;
+    if (ind == string::npos) {
+      tokens.push_back(strim(line.substr(st, string::npos)));
+      break;
+    } else {
+      tokens.push_back(strim(line.substr(st, ind - st)));
+      // if we want to store which exactly delimeter was used
+      dtokens.push_back(line[ind]);
+    }
+    st = ind+1;
+  }
+
+  if (tolower) {
+    for (int i = 0; i < tokens.size(); i++) {
+      stolower(tokens[i]);
+    }
+  }
+  return found;
+}
+
+
+int Tokenizer::tokenInt(size_t ind) {
+  return std::stoi(tokens[ind]);
+}
+
+
+double Tokenizer::tokenDouble(size_t ind) {
+  return std::stod(tokens[ind]);
+}
+
+bool readNeuFile(std::string filename, MeshData& md)
+{
+  uint32 n_number, en;
+  ifstream file(filename);
+  if (!file) {
+    LOG(WARNING) << "Can't open neu file " << filename;
+    return false;
+  }
+  md.clear();
+
+  // read nodes section
+  int nodesCount;
+  file >> nodesCount;
+  md.nodesPos.resize(nodesCount);
+  md.nodesNumbers.reserve(nodesCount);
+
+  for (int node_num(0); node_num < nodesCount; ++node_num) {
+    for (int col(0); col < 3; ++col) {
+      file >> md.nodesPos[node_num][col];
+    }
+    md.nodesNumbers.push_back(node_num+1);
+  }
+  md.cellIntData.insert(std::make_pair("MAT", std::vector<uint32>()));
+
+  // read volumes section
+  int volumesCount;
+  file >> volumesCount;
+  md.cellNumbers.reserve(volumesCount);
+  for (int volume_num(0); volume_num < volumesCount; ++volume_num) {
+    int mat_num;
+    file >> mat_num;
+    md.cellIntData["MAT"].push_back(mat_num);
+    vector<uint32> enodes;
+    enodes.resize(4);
+    for (int node_num(0); node_num < 4; ++node_num) {
+      file >> enodes[node_num];
+    }
+    for (int node_num(4); node_num < 8; ++node_num) {
+        enodes.push_back(enodes[3]);
+    }
+    md.cellNodes.push_back(enodes);
+    md.cellNumbers.push_back(volume_num + 1);
+  }
+
+  // read surfaces section
+  int surfaceCount;
+  file >> surfaceCount;
+  for (int surfNum(0); surfNum < surfaceCount; ++surfNum) {
+    int bc;
+    file >> bc;
+    string compname = std::to_string(bc);
+    if (md.feComps.find(compname) == md.feComps.end()) {
+      FEComponent comp;
+      comp.name = compname;
+      comp.type = FEComponent::ELEMENTS;
+      md.feComps.insert(std::make_pair(comp.name, comp));
+    }
+    for (int node_num(0); node_num < 3; ++node_num) {
+      int node;
+      file >> node;
+      md.feComps[compname].list.push_back(node);
+    }
+  }
+  file.close();
+  return true;
+}
+
+
+
+bool readCdbFile(std::string filename, MeshData& md) {
+  uint32 n_number, en;
+  ifstream file(filename);
+  if (!file) {
+    LOG(WARNING) << "Can't open cdb file " << filename;
+    return false;
+  }
+  md.clear();
+  // init element addition data which is specific to cdb file format
+  md.cellIntData.insert(std::make_pair("MAT", std::vector<uint32>()));
+  md.cellIntData.insert(std::make_pair("REAL", std::vector<uint32>()));
+  md.cellIntData.insert(std::make_pair("SECT", std::vector<uint32>()));
+  md.cellIntData.insert(std::make_pair("CS", std::vector<uint32>()));
+  md.cellIntData.insert(std::make_pair("SHAPE", std::vector<uint32>()));
+  md.cellIntData.insert(std::make_pair("TYPE", std::vector<uint32>()));
+
+  std::unordered_map<std::string, double> apdlParameters;
+
+  // current element type number. We keep the current element type number while proceed the apld
+  // file. Currently, this info is not used, but in the past it was used to determine element type
+  // while parsing EBLOCK blocks.
+  uint32 type = 0;
+
+  Tokenizer t;
+  // do not convert tokens to lower symbols
+  t.tolower = false;
+  t.delimiters.insert(',');
+
+  string line;
+  while (!getLine(file, line).eof()) {
+    // split line with delimiters ','
+    t.tokenize(line);
+    if (t.tokens[0] == "") continue;
+    if (t.tokens[0][0] == '/') continue;
+    if (t.tokens[0][0] == '!') continue;
+
+    if (iequals(t.tokens[0], "NBLOCK")) {
+      // example from cdb Ansys APDL 12.1:
+      //NBLOCK,6,SOLID,     9355,     9355
+      //(3i9,6e20.13)
+      //        1        0        0 7.0785325971794E+01 6.5691449317818E+01-3.6714639015390E+01
+      //
+      // from Ansys WB 17.1:
+      // nblock,3
+      // (1i9,3e20.9e3)
+      //         1    2.000000000E+000    6.000000000E+000    0.000000000E+000
+      //         2    0.000000000E+000    6.000000000E+000    0.000000000E+000
+      //
+      // from Ansys WB APDL 17.1:
+      //NBLOCK,6,SOLID,       341,       341
+      //(3i9,6e21.13e3)
+      //        1        0        0 0.0000000000000E+000 4.0000000000000E+000
+      //        2        0        0 0.0000000000000E+000 6.0000000000000E+000
+      // parse format string
+      getLine(file, line);
+      stolower(strim(line));
+      std::regex re(R"(\((\d+)i(\d+),(\d+)e(\d+)\.[0-9e]+\))");
+      std::smatch match;
+      int int_num;
+      int int_field;
+      int float_num;
+      int float_field;
+      if (std::regex_match(line, match, re)) {
+        int_num = std::stoi(match[1]);
+        int_field = std::stoi(match[2]);
+        float_num = std::stoi(match[3]);
+        float_field = std::stoi(match[4]);
+      } else {
+        LOG(FATAL) << "Don't understand nblock format string \"" << line << "\"";
+      }
+
+      // prepare array of field widths
+      std::vector<int> widths;
+      for (int i = 0; i < int_num; i++) {
+        widths.push_back(int_field);
+      }
+      for (int i = 0; i < float_num; i++) {
+        widths.push_back(float_field);
+      }
+
+      // read nodes info line by line
+      while(!getLine(file, line).eof()) {
+        // check for the end of nodal table
+        char ch = sfirstNotBlank(line);
+        if (ch == '-' || ch == 'N') break;
+
+        auto v = ssplit(line, widths, false);
+        // column no. 1 - node number
+        uint32 number = std::stoi(v[0]);
+        // Vec is initialized with zeros by default
+        Vec<3> pos;
+        // cdb file could skip last components of positions if they are equal to zero
+        if (v.size() > int_num) pos[0] = std::stod(v[int_num]);
+        if (v.size() > int_num + 1) pos[1] = std::stod(v[int_num + 1]);
+        if (v.size() > int_num + 2) pos[2] = std::stod(v[int_num + 2]);
+        md.nodesNumbers.push_back(number);
+        md.nodesPos.push_back(pos);
+      }
+    }//NBLOCK
+    else if (iequals(t.tokens[0], "EBLOCK")) {
+      // Example of record:
+      //EBLOCK,19,SOLID,      7024,      7024
+      //(19i9)
+      // ...
+      // Info from Ansys 17 help pages:
+      /*The format of the element "block" is as follows for the SOLID format:
+      Field 1 - The material number.
+      Field 2 - The element type number.
+      Field 3 - The real constant number.
+      Field 4 - The section ID attribute (beam section) number.
+      Field 5 - The element coordinate system number.
+      Field 6 - The birth/death flag.
+      Field 7 - The solid model reference number.
+      Field 8 - The element shape flag.
+      Field 9 - The number of nodes defining this element if Solkey = SOLID; otherwise, Field 9 = 0.
+      Field 10 - Not used.
+      Field 11 - The element number.
+      Fields 12-19 - The node numbers. The next line will have the additional node numbers if there
+      are more than eight.
+
+      The format without the SOLID keyword is:
+      Field 1 - The element number.
+      Field 2 - The type of section ID.
+      Field 3 - The real constant number.
+      Field 4 - The material number.
+      Field 5 - The element coordinate system number.
+      Fields 6-15 - The node numbers. The next line will have the additional node numbers if there
+      are more than ten.
+      The final line of the block will be a -1 in field 1. */
+
+      //TODO: It seems that getline keeps windows line ending
+      bool isSolid = iequals(t.tokens[2], "SOLID");
+      // parse format string
+      getLine(file, line);
+      stolower(strim(line));
+      std::regex re(R"(\((\d+)i(\d+)\))");
+      std::smatch match;
+      int int_num;
+      int int_field;
+      if (std::regex_match(line, match, re)) {
+        int_num = std::stoi(match[1]);
+        int_field = std::stoi(match[2]);
+      } else {
+        LOG(FATAL) << "Don't understand eblock format string \"" << line << "\"";
+      }
+
+      // prepare array of field widths
+      std::vector<int> widths;
+      for (int i = 0; i < int_num; i++) {
+        widths.push_back(int_field);
+      }
+      // read element data line by line
+      while(!getLine(file, line).eof()) {
+        // check for the end of element table
+        char ch = sfirstNotBlank(line);
+        if (ch == '-' || ch == 'N') break;
+
+        auto v = ssplit(line, widths, false);
+
+        // initialize values that we read from element table
+        uint32 MAT = 0;
+        uint32 ETYPE = 0;
+        uint32 REAL = 0;
+        uint32 SECT = 0;
+        uint32 CS = 0;
+        uint32 SHAPE = 0;
+        uint32 nNodes = 0;
+        uint32 elNumber = 0;
+        size_t st = 0;
+        if (isSolid) {
+          MAT = static_cast<uint32>(std::stoi(v[0]));
+          ETYPE = static_cast<uint32>(std::stoi(v[1]));
+          REAL = static_cast<uint32>(std::stoi(v[2]));
+          SECT = static_cast<uint32>(std::stoi(v[3]));
+          CS = static_cast<uint32>(std::stoi(v[4]));
+          SHAPE = static_cast<uint32>(std::stoi(v[7]));
+          nNodes = static_cast<uint32>(std::stoi(v[8]));
+          elNumber = static_cast<uint32>(std::stoi(v[10]));
+          st = 11;
+        } else {
+          elNumber = static_cast<uint32>(std::stoi(v[0]));
+          // based on my experiments on *.apdl files this field is about etype, but not about secnum
+          ETYPE = static_cast<uint32>(std::stoi(v[1]));
+          REAL = static_cast<uint32>(std::stoi(v[2]));
+          MAT = static_cast<uint32>(std::stoi(v[3]));
+          CS = static_cast<uint32>(std::stoi(v[4]));
+
+          // in this case we need to determine number of nodes by ourselfs
+          nNodes = v.size() - 5;
+          st = 5;
+        }
+        vector<uint32> enodes;
+
+        if (v.size() != st + nNodes) {
+          LOG(FATAL) << "Not enought fields to read element nodes";
+        }
+        // parse element nodes
+        for (int i = 0; i < nNodes; i++) {
+          enodes.push_back(static_cast<uint32>(std::stoi(v[st + i])));
+        }
+
+        md.cellNumbers.push_back(elNumber);
+        md.cellNodes.push_back(enodes);
+        md.cellIntData["TYPE"].push_back(ETYPE);
+        md.cellIntData["MAT"].push_back(MAT);
+        md.cellIntData["REAL"].push_back(REAL);
+        md.cellIntData["SECT"].push_back(SECT);
+        md.cellIntData["CS"].push_back(CS);
+        md.cellIntData["SHAPE"].push_back(SHAPE);
+      } // end of element table
+    }//EBLOCK
+    else if (iequals(t.tokens[0], "D")) {
+      // Fixed dof record
+      fixBC bnd;
+      bool isNumeric = false;
+      try {
+        bnd.node = t.tokenInt(1);
+        bnd.value = t.tokenDouble(3);
+        isNumeric = true;
+      } catch (const std::invalid_argument& ia) {
+        isNumeric = false;
+      }
+      if (isNumeric) {
+        // node and value field are number, we can deal with in. In other cases these can be 'all',
+        bnd.node_dof = Dof::label2dofType(t.tokens[2]);
+        md.fixBcs.push_back(bnd);
+      }
+    }//D
+    else if (iequals(t.tokens[0], "CE")) {
+      // CE - is a command to declare MPC equation
+      // How MPC looks like this in cdb file:
+      //CE,R5.0,DEFI,       2,       1,  0.00000000
+      //CE,R5.0,NODE,      1700,UX  ,  1.00000000    ,      1700,UZ  ,  1.00000000
+      // TODO: Mpc is dynamicaly allocated, but MeshData won't free this memory (this should be done
+      // in FEStorage). This is potential memory leak
+      Mpc* mpc = new Mpc();
+      mpc->b = t.tokenDouble(5);
+      int n_terms = t.tokenInt(3);
+      // read MPC terms. They are stored by 2 in a row
+      while (n_terms > 0) {
+        getLine(file, line);
+        t.tokenize(line);
+        uint16 place = 6;
+        for (int i = 0; i < std::max(n_terms, 2); i++) {
+          uint32 node = t.tokenInt(3 + 3 * i + 0);
+          Dof::dofType dof = Dof::label2dofType(t.tokens[3 + 3 * i + 1]);
+          double coef = t.tokenDouble(3 + 3 * i + 2);
+          mpc->eq.push_back(MpcTerm(node,dof,coef));
+          n_terms--;
+        }
+      }
+      md.mpcs.push_back(mpc);
+    }//CE (MPC)
+    else if (iequals(t.tokens[0], "CMBLOCK")) {
+      // Example of CMBLOCK command:
+      //CMBLOCK,BOTTOM_SIDE,NODE,      17  ! users node component definition
+      //(8i10)
+      //      5037     -5341      6330     -6352      6355      6357      6433     -6456
+      //      6459      6470     -6473      6537     -6556      6566     -6569      6633
+      //     -6652
+      if (t.tokens[1][0] != '_') {
+        // we work only with component names started not with underscore '_'
+        FEComponent comp;
+        comp.name = t.tokens[1];
+        // NODE or ELEMENT component?
+        comp.type = FEComponent::typeFromString(stoupper(t.tokens[2]));
+
+        size_t numRanges = t.tokenInt(3);
+        getLine(file, line);
+        // we need to take a format of columns "(8i10)"
+        // read how many records in a single row
+        stolower(strim(line));
+        std::regex re(R"(\((\d+)i(\d+)\))");
+        std::smatch match;
+        int int_num;
+        int int_field;
+        if (std::regex_match(line, match, re)) {
+          int_num = std::stoi(match[1]);
+          int_field = std::stoi(match[2]);
+        } else {
+          LOG(FATAL) << "Don't understand CMBLOCK format string \"" << line << "\"";
+        }
+        // prepare array of field widths
+        std::vector<int> widths;
+        for (int i = 0; i < int_num; i++) {
+          widths.push_back(int_field);
+        }
+
+        uint16 in_row = 0;
+        uint16 all = 0;
+
+        // read all ranges and calculate overall number of entities in component list
+        vector<int32> rangesVec;
+        rangesVec.reserve(numRanges);
+        getLine(file, line);
+        auto v = ssplit(line, widths, false);
+        size_t numEntity = 0;
+        while (all < numRanges) {
+           if (in_row == int_num) {
+              getLine(file, line);
+              v = ssplit(line, widths, false);
+              in_row = 0;
+           }
+           rangesVec.push_back(std::stoi(v[in_row]));
+           if (rangesVec[all] > 0) {
+              numEntity++;
+           } else {
+              // if next index is negative that means a range with step 1
+              // for example: 4 -9 --> 4 5 6 7 8 9
+              assert(all > 0);
+              assert(rangesVec[all-1] > 0);
+              numEntity += -rangesVec[all] - rangesVec[all-1];
+           }
+           in_row ++;
+           all ++;
+        }
+        // fill list of entities numbers
+        comp.list.reserve(numEntity);
+        all = 0;
+        while (all < numRanges) {
+          if (rangesVec[all] > 0) {
+            comp.list.push_back(rangesVec[all]);
+          } else {
+            for (uint32 i = rangesVec[all-1] + 1; i < static_cast<uint32> (-rangesVec[all]+1); i++) {
+              // TODO: need assert for overflow
+              comp.list.push_back(i);
+            }
+          }
+          all++;
+        }
+        md.feComps.insert(std::make_pair(comp.name, comp));
+      }
+    } //CMBLOCK
+    else if (iequals(t.tokens[0], "ET") || iequals(t.tokens[0], "ETYPE")) {
+      // Keep a track on last ET command.. We need this info to know element type.
+      // If we can convert ET argument into int, then we try to find the APDL parameter in
+      // apdlParameters dictionary. If it's not there - adpParameters will return new zero value.
+      // That means, that if parameter is not declared(somehow..) we just put type = 0
+      // NOTE: in current implementation type value is not used, but it could be used again in the
+      // future.
+      try {
+        type = t.tokenInt(1);
+      } catch (const std::invalid_argument& ia) {
+        type = static_cast<uint32>(apdlParameters[t.tokens[1]]);
+      }
+    } // ETYPE
+    else if (iequals(t.tokens[0], "TYPE")) {
+      // Keep the track on last TYPE command.. We need this info to know element type.
+      // NOTE: in current implementation type value is not used, but it could be used again in the
+      // future.
+      try {
+        type = t.tokenInt(1);
+      } catch (const std::invalid_argument& ia) {
+        type = 0;
+      }
+    } // ETYPE
+    else if (iequals(t.tokens[0], "F")) {
+      // Dof force record
+      //f,42,heat,800000.
+      loadBC bnd;
+      bool isNumeric = false;
+      try {
+        bnd.node = t.tokenInt(1);
+        bnd.value = t.tokenDouble(3);
+        isNumeric = true;
+      } catch (const std::invalid_argument& ia) {
+        isNumeric = false;
+      }
+      if (isNumeric) {
+        // node and value field are number, we can deal with in. In other cases these can be 'all',
+        // parameter or reference to table '%my_apdl_table%. We can't work with these complex cases.
+        if (iequals(t.tokens[2], "HEAT")) {
+          bnd.node_dof = Dof::TEMP;
+        } else if (iequals(t.tokens[2], "FX")) {
+          bnd.node_dof = Dof::UX;
+        } else if (iequals(t.tokens[2], "FY")) {
+          bnd.node_dof = Dof::UY;
+        } else if (iequals(t.tokens[2], "FZ")) {
+          bnd.node_dof = Dof::UZ;
+        } else {
+          bnd.node_dof = Dof::label2dofType(t.tokens[2]);
+        }
+        md.loadBcs.push_back(bnd);
+      }
+    }//F
+    else if (iequals(t.tokens[0], "*set")) {
+      // This is declaration of APDL parameters. We track all APDL numeric variables which is
+      // defined in APDL file.
+      // NOTE: currently this information is not used, but it could be useful later
+      bool isNumeric = false;
+      double value;
+      std::string parName;
+      // try to read parameter's name and it's value as double
+      // if fail - skip this expression
+      try {
+        parName = t.tokens[1];
+        value = t.tokenDouble(2);
+        isNumeric = true;
+      } catch (const std::invalid_argument& ia) {
+        isNumeric = false;
+      }
+      if(isNumeric) {
+        apdlParameters[parName] = value;
+      }
+    }
+  }
+  file.close();
+  return true;
+}

libgcdfes/FEReaders.h

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#pragma once
+#include <vector>
+#include <set>
+#include <string>
+#include <iostream>
+#include <algorithm>
+#include <regex>
+
+#include "sys.h"
+#include "Node.h"
+#include "elements/element.h"
+#include "FEStorage.h"
+#include "FESolver.h"
+
+
+using namespace nla3d;
+
+
+// neutral data structure to keep FE data.
+class MeshData {
+  public:
+    std::vector<uint32> cellNumbers;
+    std::vector<std::vector<uint32>> cellNodes;
+
+    // custom element data of different types
+    // NOTE: every element in MeshData should have the same set of custom data
+    std::map<std::string, std::vector<double>> cellDoubleData;
+    std::map<std::string, std::vector<uint32>> cellIntData;
+    std::map<std::string, std::vector<string>> cellStringData;
+
+    std::vector<uint32> nodesNumbers;
+    std::vector<Vec<3>> nodesPos;
+
+    std::vector<loadBC> loadBcs;
+    std::vector<fixBC> fixBcs;
+
+    std::vector<Mpc*> mpcs; 
+
+    std::map<std::string, FEComponent> feComps;
+
+    void clear();
+    // compress node numbers, we don't compress element numbers because FEStorage will renumber
+    // elements anyway
+    void compressNumbers();
+    // delete nodes which don't belong to any element
+    // void deleteOrphanNodes();
+    // cell is degenerated when some it's nodes have the same numbers
+    std::vector<uint32> getDegeneratedCells();
+    std::vector<uint32> getCellsByAttribute(std::string atr_name, uint32 atr_val);
+
+  private:
+
+};
+
+std::string& strim(std::string& str);
+// std::move semantic to use with rvalues
+std::string&& strim(std::string&& str);
+// return first non-space char
+char sfirstNotBlank(const string& str);
+
+std::string& stolower(std::string& str);
+std::string& stoupper(std::string& str);
+
+// split `line` to substrings with `widths` lengths, if `strict` is false then last substrings could
+// be missed
+std::vector<std::string> ssplit(const std::string& line, const std::vector<int>& widths,
+                                  bool strict = true);
+
+// case insensitive string comparison
+bool iequals(const string& a, const string& b);
+
+// get line from the `is`, supported all types of line endings
+std::istream& getLine(std::istream& is, std::string& t);
+
+
+// class to split line into substrings with delimiters.
+class Tokenizer {
+  public:
+    int tokenize(const std::string& line);
+    int tokenInt(size_t ind);
+    double tokenDouble(size_t ind);
+    std::string& getTokenString(size_t ind);
+
+    std::set<char> delimiters;
+    std::vector<std::string> tokens;
+    std::vector<char> dtokens;
+    bool tolower = false;
+};
+
+
+// read Ansys Mechanical APDL *.cdb file. Nodes, Elements, Displacement BC, Force BC and MPC
+// (Constraint equations) is supported
+bool readCdbFile(std::string filename, MeshData& md);
+
+// read netgen neutral volume mesh format *.neu file. 
+bool readNeuFile(std::string filename, MeshData& md);

libgcdfes/libgcdfes.cbp

+<?xml version="1.0" encoding="UTF-8" standalone="yes" ?>
+<CodeBlocks_project_file>
+	<FileVersion major="1" minor="6" />
+	<Project>
+		<Option title="libgcdfes" />
+		<Option pch_mode="2" />
+		<Option compiler="gcc" />
+		<Build>
+			<Target title="Debug">
+				<Option output="bin/Debug/gcdfes" prefix_auto="1" extension_auto="1" />
+				<Option object_output="obj/Debug/" />
+				<Option type="3" />
+				<Option compiler="gcc" />
+				<Option createDefFile="1" />
+				<Option createStaticLib="1" />
+				<Compiler>
+					<Add option="-m64" />
+					<Add option="-g" />
+					<Add option="-fPIC" />
+					<Add option="-DNLA3D_USE_MKL" />
+					<Add option="-DNLA3D_USE_PARALUTION" />
+					<Add directory="../../cuda_nla3d" />
+					<Add directory="../comsdk/dev" />
+					<Add directory="../easylogging" />
+					<Add directory="../math" />
+					<Add directory="../nla3d" />
+					<Add directory="/opt/intel/mkl/include" />
+					<Add directory="/usr/include/eigen3" />
+					<Add directory="../libgcdfes" />
+				</Compiler>
+				<Linker>
+					<Add option="-m64" />
+					<Add option="-lcomsdk" />
+					<Add option="-lnla3d" />
+					<Add option="-lmath" />
+					<Add option="-lparalution" />
+					<Add option="-lmkl_intel_lp64" />
+					<Add option="-lmkl_sequential" />
+					<Add option="-lmkl_core" />
+					<Add option="-lpthread" />
+					<Add option="-lm" />
+					<Add option="-ldl" />
+					<Add option="-lboost_system" />
+					<Add directory="../comsdk/rls" />
+					<Add directory="../math/bin/Debug" />
+					<Add directory="../nla3d/bin/Debug" />
+					<Add directory="../paralution-1.1.0/build/lib" />
+					<Add directory="/opt/intel/mkl/lib/intel64" />
+				</Linker>
+			</Target>
+			<Target title="Release">
+				<Option output="bin/Release/gcdfes" prefix_auto="1" extension_auto="1" />
+				<Option object_output="obj/Release/" />
+				<Option type="3" />
+				<Option compiler="gcc" />
+				<Option createDefFile="1" />
+				<Option createStaticLib="1" />
+				<Compiler>
+					<Add option="-O2" />
+					<Add option="-m64" />
+					<Add option="-fPIC" />
+					<Add option="-DNLA3D_USE_MKL" />
+					<Add option="-DNLA3D_USE_PARALUTION" />
+					<Add directory="../../cuda_nla3d" />
+					<Add directory="../comsdk/dev" />
+					<Add directory="../easylogging" />
+					<Add directory="../math" />
+					<Add directory="../nla3d" />
+					<Add directory="/opt/intel/mkl/include" />
+					<Add directory="/usr/include/eigen3" />
+					<Add directory="../libgcdfes" />
+				</Compiler>
+				<Linker>
+					<Add option="-s" />
+					<Add option="-m64" />
+					<Add option="-lcomsdk" />
+					<Add option="-lnla3d" />
+					<Add option="-lmath" />
+					<Add option="-lparalution" />
+					<Add option="-lmkl_intel_lp64" />
+					<Add option="-lmkl_sequential" />
+					<Add option="-lmkl_core" />
+					<Add option="-lpthread" />
+					<Add option="-lm" />
+					<Add option="-ldl" />
+					<Add option="-lboost_system" />
+					<Add directory="../comsdk/rls" />
+					<Add directory="../math/bin/Release" />
+					<Add directory="../nla3d/bin/Release" />
+					<Add directory="../paralution-1.1.0/build/lib" />
+					<Add directory="/opt/intel/mkl/lib/intel64" />
+				</Linker>
+			</Target>
+		</Build>
+		<Compiler>
+			<Add option="-Wall" />
+			<Add option="-fexceptions" />
+		</Compiler>
+		<Unit filename="FEReaders.cpp" />
+		<Unit filename="gcdfes.cpp" />
+		<Unit filename="gcdfes.h" />
+		<Extensions>
+			<code_completion />
+			<envvars />
+			<lib_finder disable_auto="1" />
+			<debugger />
+		</Extensions>
+	</Project>
+</CodeBlocks_project_file>

libgcdfes/makefile

+#------------------------------------------------------------------------------#
+# This makefile was generated by 'cbp2make' tool rev.147                       #
+#------------------------------------------------------------------------------#
+
+
+WORKDIR = `pwd`
+
+CC = gcc
+CXX = g++
+AR = ar
+LD = g++
+WINDRES = windres
+
+INC = 
+CFLAGS = -Wall -fexceptions
+RESINC = 
+LIBDIR = 
+LIB = 
+LDFLAGS = 
+
+INC_DEBUG = $(INC) -I../../cuda_nla3d -I../comsdk/dev -I../easylogging -I../math -I../nla3d -I/opt/intel/mkl/include -I/usr/include/eigen3 -I../libgcdfes
+CFLAGS_DEBUG = $(CFLAGS) -m64 -g -fPIC -DNLA3D_USE_MKL -DNLA3D_USE_PARALUTION
+RESINC_DEBUG = $(RESINC)
+RCFLAGS_DEBUG = $(RCFLAGS)
+LIBDIR_DEBUG = $(LIBDIR) -L../comsdk/rls -L../math/bin/Debug -L../nla3d/bin/Debug -L../paralution-1.1.0/build/lib -L/opt/intel/mkl/lib/intel64
+LIB_DEBUG = $(LIB)
+LDFLAGS_DEBUG = $(LDFLAGS) -m64 -lcomsdk -lnla3d -lmath -lparalution -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl -lboost_system
+OBJDIR_DEBUG = obj/Debug
+DEP_DEBUG = 
+OUT_DEBUG = bin/Debug/gcdfes.so
+
+INC_RELEASE = $(INC) -I../../cuda_nla3d -I../comsdk/dev -I../easylogging -I../math -I../nla3d -I/opt/intel/mkl/include -I/usr/include/eigen3 -I../libgcdfes
+CFLAGS_RELEASE = $(CFLAGS) -O2 -m64 -fPIC -DNLA3D_USE_MKL -DNLA3D_USE_PARALUTION
+RESINC_RELEASE = $(RESINC)
+RCFLAGS_RELEASE = $(RCFLAGS)
+LIBDIR_RELEASE = $(LIBDIR) -L../comsdk/rls -L../math/bin/Release -L../nla3d/bin/Release -L../paralution-1.1.0/build/lib -L/opt/intel/mkl/lib/intel64
+LIB_RELEASE = $(LIB)
+LDFLAGS_RELEASE = $(LDFLAGS) -s -m64 -lcomsdk -lnla3d -lmath -lparalution -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl -lboost_system
+OBJDIR_RELEASE = obj/Release
+DEP_RELEASE = 
+OUT_RELEASE = bin/Release/gcdfes.so
+
+OBJ_DEBUG = $(OBJDIR_DEBUG)/FEReaders.o $(OBJDIR_DEBUG)/gcdfes.o
+
+OBJ_RELEASE = $(OBJDIR_RELEASE)/FEReaders.o $(OBJDIR_RELEASE)/gcdfes.o
+
+all: debug release
+
+clean: clean_debug clean_release
+
+before_debug: 
+	test -d bin/Debug || mkdir -p bin/Debug
+	test -d $(OBJDIR_DEBUG) || mkdir -p $(OBJDIR_DEBUG)
+
+after_debug: 
+
+debug: before_debug out_debug after_debug
+
+out_debug: before_debug $(OBJ_DEBUG) $(DEP_DEBUG)
+	$(LD) -shared $(LIBDIR_DEBUG) $(OBJ_DEBUG)  -o $(OUT_DEBUG) $(LDFLAGS_DEBUG) $(LIB_DEBUG)
+
+$(OBJDIR_DEBUG)/FEReaders.o: FEReaders.cpp
+	$(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c FEReaders.cpp -o $(OBJDIR_DEBUG)/FEReaders.o
+
+$(OBJDIR_DEBUG)/gcdfes.o: gcdfes.cpp
+	$(CXX) $(CFLAGS_DEBUG) $(INC_DEBUG) -c gcdfes.cpp -o $(OBJDIR_DEBUG)/gcdfes.o
+
+clean_debug: 
+	rm -f $(OBJ_DEBUG) $(OUT_DEBUG)
+	rm -rf bin/Debug
+	rm -rf $(OBJDIR_DEBUG)
+
+before_release: 
+	test -d bin/Release || mkdir -p bin/Release
+	test -d $(OBJDIR_RELEASE) || mkdir -p $(OBJDIR_RELEASE)
+
+after_release: 
+
+release: before_release out_release after_release
+
+out_release: before_release $(OBJ_RELEASE) $(DEP_RELEASE)
+	$(LD) -shared $(LIBDIR_RELEASE) $(OBJ_RELEASE)  -o $(OUT_RELEASE) $(LDFLAGS_RELEASE) $(LIB_RELEASE)
+
+$(OBJDIR_RELEASE)/FEReaders.o: FEReaders.cpp
+	$(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c FEReaders.cpp -o $(OBJDIR_RELEASE)/FEReaders.o
+
+$(OBJDIR_RELEASE)/gcdfes.o: gcdfes.cpp
+	$(CXX) $(CFLAGS_RELEASE) $(INC_RELEASE) -c gcdfes.cpp -o $(OBJDIR_RELEASE)/gcdfes.o
+
+clean_release: 
+	rm -f $(OBJ_RELEASE) $(OUT_RELEASE)
+	rm -rf bin/Release
+	rm -rf $(OBJDIR_RELEASE)
+
+.PHONY: before_debug after_debug clean_debug before_release after_release clean_release
+