Fem: Add frd format converter to VTK

src/Mod/Fem/App/AppFemPy.cpp

@@ -66,6 +66,7 @@ public:
                            &Module::read,
                            "Read a mesh from a file and returns a Mesh object.");
 #ifdef FC_USE_VTK
+        add_varargs_method("frdToVTK", &Module::frdToVTK, "Convert a .frd result file to VTK file");
         add_varargs_method("readResult",
                            &Module::readResult,
                            "Read a CFD or Mechanical result (auto detect) from a file (file format "
@@ -248,6 +249,21 @@ private:
     }
 
 #ifdef FC_USE_VTK
+    Py::Object frdToVTK(const Py::Tuple& args)
+    {
+        char* filename = nullptr;
+
+        if (!PyArg_ParseTuple(args.ptr(), "et", "utf-8", &filename)) {
+            throw Py::Exception();
+        }
+        std::string encodedName = std::string(filename);
+        PyMem_Free(filename);
+
+        FemVTKTools::frdToVTK(encodedName.c_str());
+
+        return Py::None();
+    }
+
     Py::Object readResult(const Py::Tuple& args)
     {
         char* fileName = nullptr;

src/Mod/Fem/App/FemVTKTools.cpp

@@ -25,6 +25,7 @@
 
 #ifndef _PreComp_
 #include <Python.h>
+#include <charconv>
 #include <cmath>
 #include <cstdlib>
 #include <map>
@@ -41,6 +42,7 @@
 #include <vtkHexahedron.h>
 #include <vtkIdList.h>
 #include <vtkLine.h>
+#include <vtkMultiBlockDataSet.h>
 #include <vtkPointData.h>
 #include <vtkPyramid.h>
 #include <vtkQuad.h>
@@ -56,6 +58,7 @@
 #include <vtkUnsignedCharArray.h>
 #include <vtkUnstructuredGrid.h>
 #include <vtkWedge.h>
+#include <vtkXMLMultiBlockDataWriter.h>
 #include <vtkXMLPUnstructuredGridReader.h>
 #include <vtkXMLUnstructuredGridReader.h>
 #include <vtkXMLUnstructuredGridWriter.h>
@@ -1029,4 +1032,583 @@ void FemVTKTools::exportFreeCADResult(const App::DocumentObject* result,
     Base::Console().Log("End: Create VTK result data from FreeCAD result data.\n");
 }
 
+
+namespace FRDReader
+{
+
+// number of nodes per CalculiX element type: {type, nodes}
+std::map<int, unsigned int> mapCcxTypeNodes = {
+    {11, 2},
+    {12, 3},
+    {7, 3},
+    {8, 6},
+    {9, 4},
+    {10, 8},
+    {3, 4},
+    {6, 10},
+    {1, 8},
+    {4, 20},
+    {2, 6},
+    {5, 15},
+};
+
+// map CalculiX nodes order to Vtk order
+std::map<int, std::vector<int>> mapCcxToVtk = {
+    {VTK_LINE, {0, 1}},
+    {VTK_QUADRATIC_EDGE, {0, 1, 2}},
+    {VTK_TRIANGLE, {0, 1, 2}},
+    {VTK_QUADRATIC_TRIANGLE, {0, 1, 2, 3, 4, 5}},
+    {VTK_QUAD, {0, 1, 2, 3}},
+    {VTK_QUADRATIC_QUAD, {0, 1, 2, 3, 4, 5, 6, 7}},
+    {VTK_TETRA, {0, 1, 2, 3}},
+    {VTK_QUADRATIC_TETRA, {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}},
+    {VTK_HEXAHEDRON, {0, 1, 2, 3, 4, 5, 6, 7}},
+    {VTK_QUADRATIC_HEXAHEDRON,
+     {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15}},
+    {VTK_WEDGE, {0, 1, 2, 3, 4, 5}},
+    {VTK_QUADRATIC_WEDGE, {0, 1, 2, 3, 4, 5, 6, 7, 8, 12, 13, 14, 9, 10, 11}}};
+
+// give position of first non-blank character of string_view
+size_t getFirstNotBlankPos(const std::string_view& view)
+{
+    size_t pos = view.find_first_not_of(" ");
+    if (pos == std::string_view::npos) {
+        pos = 0;
+    }
+
+    return pos;
+}
+
+// get n-digits value from string_view
+// not used until libc++ std::from_chars supports double values
+// template<typename T>
+// void valueFromLine(const std::string_view::const_iterator& it, int digits, T& value)
+//{
+//    std::string_view sub(it, digits);
+//    auto pos = getFirstNotBlankPos(sub);
+//    std::from_chars(sub.data() + pos, sub.data() + digits, value, 10);
+//}
+template<typename T>
+void valueFromLine(const std::string_view::iterator& it, int digits, T& value)
+{
+    std::string_view sub(&*it, digits);
+    value = std::strtol(sub.data(), nullptr, 10);
+}
+template<>
+void valueFromLine<double>(const std::string_view::iterator& it, int digits, double& value)
+{
+    std::string_view sub(&*it, digits);
+    value = std::strtof(sub.data(), nullptr);
+}
+
+// add cell from sorted nodes
+template<typename T>
+void addCell(vtkSmartPointer<vtkCellArray>& cellArray, const std::vector<int>& topoElem)
+{
+    vtkSmartPointer<T> cell = vtkSmartPointer<T>::New();
+    cell->GetPointIds()->SetNumberOfIds(topoElem.size());
+    int type = cell->GetCellType();
+    for (size_t i = 0; i < topoElem.size(); ++i) {
+        cell->GetPointIds()->SetId(i, topoElem[mapCcxToVtk[type][i]]);
+    }
+    cellArray->InsertNextCell(cell);
+}
+
+// fill cell array
+void fillCell(vtkSmartPointer<vtkCellArray>& cellArray,
+              std::vector<int>& topoElem,
+              std::vector<int>& vtkType,
+              int elemType)
+{
+    switch (elemType) {
+        case 1:
+            addCell<vtkHexahedron>(cellArray, topoElem);
+            vtkType.emplace_back(VTK_HEXAHEDRON);
+            break;
+        case 2:
+            addCell<vtkWedge>(cellArray, topoElem);
+            vtkType.emplace_back(VTK_WEDGE);
+            break;
+        case 3:
+            addCell<vtkTetra>(cellArray, topoElem);
+            vtkType.emplace_back(VTK_TETRA);
+            break;
+        case 4:
+            addCell<vtkQuadraticHexahedron>(cellArray, topoElem);
+            vtkType.emplace_back(VTK_QUADRATIC_HEXAHEDRON);
+            break;
+        case 5:
+            addCell<vtkQuadraticWedge>(cellArray, topoElem);
+            vtkType.emplace_back(VTK_QUADRATIC_WEDGE);
+            break;
+        case 6:
+            addCell<vtkQuadraticTetra>(cellArray, topoElem);
+            vtkType.emplace_back(VTK_QUADRATIC_TETRA);
+            break;
+        case 7:
+            addCell<vtkTriangle>(cellArray, topoElem);
+            vtkType.emplace_back(VTK_TRIANGLE);
+            break;
+        case 8:
+            addCell<vtkQuadraticTriangle>(cellArray, topoElem);
+            vtkType.emplace_back(VTK_QUADRATIC_TRIANGLE);
+            break;
+        case 9:
+            addCell<vtkQuad>(cellArray, topoElem);
+            vtkType.emplace_back(VTK_QUAD);
+            break;
+        case 10:
+            addCell<vtkQuadraticQuad>(cellArray, topoElem);
+            vtkType.emplace_back(VTK_QUADRATIC_QUAD);
+            break;
+        case 11:
+            addCell<vtkLine>(cellArray, topoElem);
+            vtkType.emplace_back(VTK_LINE);
+            break;
+        case 12:
+            addCell<vtkQuadraticEdge>(cellArray, topoElem);
+            vtkType.emplace_back(VTK_QUADRATIC_EDGE);
+            break;
+    }
+}
+
+struct FRDResultInfo
+{
+    double value;
+    long numNodes;
+    int analysisType;
+    int step;
+    int indicator;
+};
+
+// get number of digits from format indicator
+int getDigits(int indicator)
+{
+    int digits = 0;
+    if (indicator == 0) {
+        // short
+        digits = 5;
+    }
+    else if (indicator == 1) {
+        // long
+        digits = 10;
+    }
+
+    return digits;
+}
+
+// get position of scalar entities in line result vector
+std::vector<size_t> identifyScalarEntities(const std::vector<std::vector<int>> entities)
+{
+    std::vector<size_t> pos;
+    for (auto it = entities.begin(); it != entities.end(); ++it) {
+        // check type == 1 or component < 1
+        if ((*it)[0] == 1 || (*it)[1] < 1) {
+            pos.emplace_back(it - entities.begin());
+        }
+    }
+    return pos;
+}
+
+// read nodes and fill vtkPoints object
+std::map<int, int>
+readNodes(std::ifstream& ifstr, const std::string& lines, vtkSmartPointer<vtkPoints>& points)
+{
+    std::string keyCode = "    2C";
+    std::string keyCodeCoord = " -1";
+    long numNodes;
+    int indicator;
+    int node;
+    long nodeID = 0;
+
+    // frd file might have nodes that are not numbered starting from zero.
+    // Use the map to identify them
+    std::map<int, int> mapNodes;
+
+    std::string_view view {lines};
+    std::string_view sub = view.substr(keyCode.length() + 18);
+
+    valueFromLine(sub.begin(), 12, numNodes);
+
+    sub = sub.substr(12 + 37);
+    valueFromLine(sub.begin(), 1, indicator);
+    int digits = getDigits(indicator);
+
+    points->SetNumberOfPoints(numNodes);
+
+    std::string line;
+    while (nodeID < numNodes && std::getline(ifstr, line)) {
+        std::vector<double> coords;
+        std::string_view view {line};
+        if (view.rfind(keyCodeCoord, 0) == 0) {
+            std::string_view v(line.data() + keyCodeCoord.length(), digits);
+            valueFromLine(v.begin(), digits, node);
+
+            std::string_view vi = view.substr(keyCodeCoord.length() + digits);
+            double value;
+            for (auto it = vi.begin(); it != vi.end(); it += 12) {
+                valueFromLine(it, 12, value);
+                coords.emplace_back(value);
+            }
+        }
+
+        points->SetPoint(nodeID, coords.data());
+        mapNodes[node] = nodeID++;
+    }
+
+    return mapNodes;
+}
+
+// fill elements and fill cell array
+std::vector<int> readElements(std::ifstream& ifstr,
+                              const std::string& lines,
+                              const std::map<int, int>& mapNodes,
+                              vtkSmartPointer<vtkCellArray>& cellArray)
+{
+    std::string line;
+    std::string keyCode = "    3C";
+    std::string keyCodeType = " -1";
+    std::string keyCodeNodes = " -2";
+    long numElem;
+    int indicator;
+    int elem;
+    long elemID = 0;
+    // element info: {type, group, material}
+    std::vector<int> info(3);
+    std::map<int, int> mapElem;
+    std::vector<int> topoElem;
+    std::vector<int> vtkType;
+
+    std::string_view view {lines};
+
+    std::string_view sub = view.substr(keyCode.length() + 18);
+    valueFromLine(sub.begin(), 12, numElem);
+
+    sub = sub.substr(12 + 37);
+    valueFromLine(sub.begin(), 1, indicator);
+    int digits = getDigits(indicator);
+    while (elemID < numElem && std::getline(ifstr, line)) {
+        std::string_view view {line};
+        if (view.rfind(keyCodeType, 0) == 0) {
+            std::string_view v(line.data() + keyCodeType.length());
+            valueFromLine(v.begin(), digits, elem);
+            v = v.substr(digits);
+            std::string_view::iterator it1;
+            std::vector<int>::iterator it2;
+            for (it1 = v.begin(), it2 = info.begin(); it1 != v.end() && it2 != info.end();
+                 it1 += 5, ++it2) {
+                valueFromLine(it1, 5, *it2);
+            }
+        }
+        if (view.rfind(keyCodeNodes, 0) == 0) {
+            std::string_view vi = view.substr(keyCodeNodes.length());
+            int node;
+            for (auto it = vi.begin(); it != vi.end(); it += digits) {
+                valueFromLine(it, digits, node);
+                topoElem.emplace_back(mapNodes.at(node));
+            }
+
+            // add cell to cellArray
+            if (topoElem.size() == mapCcxTypeNodes[info[0]]) {
+                fillCell(cellArray, topoElem, vtkType, info[0]);
+                topoElem.clear();
+                mapElem[elem] = elemID++;
+            }
+        }
+    }
+    return vtkType;
+}
+
+// read parameter header (not used)
+void readParameter(std::ifstream& ifstr, const std::string& line)
+{
+    // do nothing
+    (void)ifstr;
+    (void)line;
+}
+
+// read first header from nodal result block
+void readResultInfo(std::ifstream& ifstr, const std::string& lines, FRDResultInfo& info)
+{
+    (void)ifstr;
+
+    std::string keyCode = "  100C";
+
+    std::string_view view {lines};
+    std::string_view sub = view.substr(keyCode.length() + 6);
+    valueFromLine(sub.begin(), 12, info.value);
+
+    sub = sub.substr(12);
+    valueFromLine(sub.begin(), 12, info.numNodes);
+
+    sub = sub.substr(12 + 20);
+    valueFromLine(sub.begin(), 2, info.analysisType);
+
+    sub = sub.substr(2);
+    valueFromLine(sub.begin(), 5, info.step);
+
+    sub = sub.substr(5 + 10);
+    valueFromLine(sub.begin(), 2, info.indicator);
+}
+
+// read result from nodal result block and add result array to grid
+void readResults(std::ifstream& ifstr,
+                 const std::string& lines,
+                 const std::map<int, int>& mapNodes,
+                 const FRDResultInfo& info,
+                 vtkSmartPointer<vtkUnstructuredGrid>& grid)
+{
+    int digits = getDigits(info.indicator);
+    (void)lines;
+
+    // get dataset info, start with " -4"
+    std::string line;
+    std::string keyDataSet = " -4";
+    unsigned int numComps;
+    std::getline(ifstr, line);
+    std::string_view view = line;
+    std::string_view sub = view.substr(keyDataSet.length() + 2);
+    std::string dataSetName {sub.substr(0, 8)};
+    // remove trailing spaces
+    dataSetName.erase(dataSetName.find_last_not_of(" ") + 1);
+    sub = sub.substr(8);
+    valueFromLine(sub.begin(), 5, numComps);
+
+    // get entity info
+    std::string keyEntity = " -5";
+    std::vector<std::string> entityNames;
+    // type: 1: scalar; 2: vector; 4: matrix; 12: vector (3 amp - 3 phase); 14: tensor (6 amp - 6
+    // phase) {type, row, col, exist}
+    std::vector<std::vector<int>> entityTypes;
+    unsigned int countComp = 0;
+    while (countComp < numComps && std::getline(ifstr, line)) {
+        std::string_view view {line};
+        if (view.rfind(keyEntity, 0) == 0) {
+            sub = view.substr(keyEntity.length() + 2);
+            std::string en {sub.substr(0, 8)};
+            // remove trailing spaces
+            en.erase(en.find_last_not_of(" ") + 1);
+            std::vector<int> et = {0, 0, 0, 0};
+            // fill entityType, ignore MENU: "    1"
+            sub = sub.substr(8 + 5, 4 * 5);
+            std::string_view::iterator it1;
+            std::vector<int>::iterator it2;
+            for (it1 = sub.begin(), it2 = et.begin(); it1 != sub.end() && it2 != et.end();
+                 (it1 += 5), ++it2) {
+                valueFromLine(it1, digits, *it2);
+            }
+
+            if (et[3] == 0) {
+                // ignore predefined entity
+                entityNames.emplace_back(en);
+                entityTypes.emplace_back(et);
+            }
+            ++countComp;
+        }
+    }
+
+    // used components
+    numComps = entityNames.size();
+
+    // enter in node values block
+    std::string code1 = " -1";
+    std::string code2 = " -2";
+    int node;
+    double value;
+    std::vector<double> vecValues;
+    std::vector<double> scaValues;
+    std::vector<int> nodes;
+    int countNodes = 0;
+    int countScaPos;
+    // result block could have both vector/matrix and scalar components
+    // save each scalars entity in his own array
+    auto scalarPos = identifyScalarEntities(entityTypes);
+    // array for vector entities (if needed)
+    vtkSmartPointer<vtkDoubleArray> vecArray = vtkSmartPointer<vtkDoubleArray>::New();
+    // arrays for scalar entities (if needed)
+    std::vector<vtkSmartPointer<vtkDoubleArray>> scaArrays;
+    for (size_t i = 0; i < scalarPos.size(); ++i) {
+        scaArrays.emplace_back(vtkSmartPointer<vtkDoubleArray>::New());
+    }
+
+    vecArray->SetNumberOfComponents(numComps - scalarPos.size());
+    vecArray->SetNumberOfTuples(mapNodes.size());
+    vecArray->SetName(dataSetName.c_str());
+    // set all values to zero
+    for (int i = 0; i < vecArray->GetNumberOfComponents(); ++i) {
+        vecArray->FillComponent(i, 0.0);
+    }
+    //    vecArray->Fill(0.0);
+    for (size_t i = 0; i < scaArrays.size(); ++i) {
+        scaArrays[i]->SetNumberOfComponents(1);
+        scaArrays[i]->SetNumberOfTuples(mapNodes.size());
+        std::string name = entityNames[scalarPos[i]];
+        scaArrays[i]->SetName(name.c_str());
+        //        scaArrays[i]->Fill(0.0);
+        for (int j = 0; j < scaArrays[i]->GetNumberOfComponents(); ++j) {
+            scaArrays[i]->FillComponent(j, 0.0);
+        }
+    }
+
+    while (countNodes < info.numNodes && std::getline(ifstr, line)) {
+        std::string_view view {line};
+        if (view.rfind(code1, 0) == 0) {
+            sub = view.substr(code1.length());
+            valueFromLine(sub.begin(), digits, node);
+            // clear values vector for each node result block
+            vecValues.clear();
+            scaValues.clear();
+            countScaPos = 0;
+            try {
+                // result nodes could not exist in .frd file due to element expansion
+                // so mapNodes.at() could throw an exception
+                nodes.emplace_back(mapNodes.at(node));
+                sub = sub.substr(digits);
+                for (auto it = sub.begin(); it != sub.end(); it += 12, ++countScaPos) {
+                    valueFromLine(it, 12, value);
+                    // search if value is scalar or vector/matrix component
+                    auto pos = std::find(scalarPos.begin(), scalarPos.end(), countScaPos);
+                    if (pos == scalarPos.end()) {
+                        vecValues.emplace_back(value);
+                    }
+                    else {
+                        scaValues.emplace_back(value);
+                    }
+                }
+            }
+            catch (const std::out_of_range& ex) {
+                Base::Console().Warning("Invalid node: %d\n", node);
+            }
+            ++countNodes;
+        }
+        else if (view.rfind(code2, 0) == 0) {
+            sub = view.substr(code2.length() + digits);
+            for (auto it = sub.begin(); it != sub.end(); it += 12) {
+                valueFromLine(it, 12, value);
+                // search if value is scalar or vector/matrix component
+                auto pos = std::find(scalarPos.begin(), scalarPos.end(), countScaPos);
+                if (pos == scalarPos.end()) {
+                    vecValues.emplace_back(value);
+                }
+                else {
+                    scaValues.emplace_back(value);
+                }
+            }
+        }
+        if ((vecValues.size() + scaValues.size()) == numComps) {
+            if (!vecValues.empty()) {
+                vecArray->SetTuple(mapNodes.at(node), vecValues.data());
+            }
+            if (!scaValues.empty()) {
+                std::vector<vtkSmartPointer<vtkDoubleArray>>::iterator it1;
+                std::vector<double>::iterator it2;
+                for (it1 = scaArrays.begin(), it2 = scaValues.begin();
+                     it1 != scaArrays.end() && it2 != scaValues.end();
+                     ++it1, ++it2) {
+                    (*it1)->SetTuple1(mapNodes.at(node), *it2);
+                }
+            }
+        }
+    }
+
+    // add vecArray only if not all scalars
+    if (numComps != scalarPos.size()) {
+        grid->GetPointData()->AddArray(vecArray);
+    }
+    for (auto& s : scaArrays) {
+        grid->GetPointData()->AddArray(s);
+    }
+}
+
+vtkSmartPointer<vtkMultiBlockDataSet> readFRD(std::ifstream& ifstr)
+{
+    auto points = vtkSmartPointer<vtkPoints>::New();
+    auto cells = vtkSmartPointer<vtkCellArray>::New();
+    auto multiBlock = vtkSmartPointer<vtkMultiBlockDataSet>::New();
+    vtkSmartPointer<vtkUnstructuredGrid> grid;
+    std::map<std::pair<int, int>, vtkSmartPointer<vtkUnstructuredGrid>> grids;
+    std::string line;
+    std::map<int, int> mapNodes;
+    std::vector<int> cellTypes;
+
+    while (std::getline(ifstr, line)) {
+        std::string keyCode = "    2C";
+        std::string_view view = line;
+
+        if (view.rfind(keyCode, 0) == 0) {
+            // read nodes block
+            mapNodes = readNodes(ifstr, line, points);
+        }
+        keyCode = "    3C";
+        if (view.rfind(keyCode, 0) == 0) {
+            // read elements block
+            cellTypes = readElements(ifstr, line, mapNodes, cells);
+        }
+        keyCode = "    1P";
+        if (view.rfind(keyCode, 0) == 0) {
+            // read parameter
+            readParameter(ifstr, line);
+        }
+        keyCode = "  100C";
+        if (view.rfind(keyCode, 0) == 0) {
+            // read result info block
+            FRDResultInfo info;
+            readResultInfo(ifstr, line, info);
+            auto it = grids.find(std::pair<int, int>(info.step, info.analysisType));
+            if (it == grids.end()) {
+                grid = vtkSmartPointer<vtkUnstructuredGrid>::New();
+                grid->SetPoints(points);
+                grid->SetCells(cellTypes.data(), cells);
+                grids[std::pair<int, int>(info.step, info.analysisType)] = grid;
+            }
+            else {
+                grid = (*it).second;
+            }
+            // read result entries and node results
+            readResults(ifstr, line, mapNodes, info, grid);
+        }
+    }
+    int i = 0;
+
+    multiBlock->SetNumberOfBlocks(grids.size());
+    for (const auto& g : grids) {
+        multiBlock->SetBlock(i, g.second);
+        ++i;
+    }
+
+    // save points and elements even without results
+    if (grids.empty()) {
+        grid = vtkSmartPointer<vtkUnstructuredGrid>::New();
+        grid->SetPoints(points);
+        grid->SetCells(cellTypes.data(), cells);
+        multiBlock->SetBlock(0, grid);
+    }
+
+    return multiBlock;
+}
+
+}  // namespace FRDReader
+
+void FemVTKTools::frdToVTK(const char* filename)
+{
+    Base::FileInfo fi(filename);
+
+    if (!fi.isReadable()) {
+        throw Base::FileException("File to load not existing or not readable", filename);
+    }
+
+    std::ifstream ifstr(filename, std::ios::in | std::ios::binary);
+
+    vtkSmartPointer<vtkMultiBlockDataSet> multiBlock = FRDReader::readFRD(ifstr);
+
+    std::string dir = fi.dirPath();
+
+    auto writer = vtkSmartPointer<vtkXMLMultiBlockDataWriter>::New();
+
+    std::string blockFile = dir + "/" + fi.fileNamePure() + "." + writer->GetDefaultFileExtension();
+    writer->SetFileName(blockFile.c_str());
+    writer->SetInputData(multiBlock);
+    writer->Update();
+}
+
 }  // namespace Fem

src/Mod/Fem/App/FemVTKTools.h

@@ -71,6 +71,8 @@ public:
 
     // write FemResult (activeObject if res= NULL) to vtkUnstructuredGrid dataset file
     static void writeResult(const char* filename, const App::DocumentObject* res = nullptr);
+
+    static void frdToVTK(const char* filename);
 };
 }  // namespace Fem