From e84eb920dcd396d4daafc5a9501dfd1a93d557c1 Mon Sep 17 00:00:00 2001 From: marioalexis Date: Wed, 29 Jan 2025 19:18:29 -0300 Subject: [PATCH] Fem: Add frd format converter to VTK --- src/Mod/Fem/App/AppFemPy.cpp | 16 + src/Mod/Fem/App/FemVTKTools.cpp | 582 ++++++++++++++++++++++++++++++++ src/Mod/Fem/App/FemVTKTools.h | 2 + 3 files changed, 600 insertions(+) diff --git a/src/Mod/Fem/App/AppFemPy.cpp b/src/Mod/Fem/App/AppFemPy.cpp index b758250e02..bbbe4520d4 100644 --- a/src/Mod/Fem/App/AppFemPy.cpp +++ b/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; diff --git a/src/Mod/Fem/App/FemVTKTools.cpp b/src/Mod/Fem/App/FemVTKTools.cpp index ff9550584c..522b461e0c 100644 --- a/src/Mod/Fem/App/FemVTKTools.cpp +++ b/src/Mod/Fem/App/FemVTKTools.cpp @@ -25,6 +25,7 @@ #ifndef _PreComp_ #include +#include #include #include #include @@ -41,6 +42,7 @@ #include #include #include +#include #include #include #include @@ -56,6 +58,7 @@ #include #include #include +#include #include #include #include @@ -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 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> 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 +// 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 +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(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 +void addCell(vtkSmartPointer& cellArray, const std::vector& topoElem) +{ + vtkSmartPointer cell = vtkSmartPointer::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& cellArray, + std::vector& topoElem, + std::vector& vtkType, + int elemType) +{ + switch (elemType) { + case 1: + addCell(cellArray, topoElem); + vtkType.emplace_back(VTK_HEXAHEDRON); + break; + case 2: + addCell(cellArray, topoElem); + vtkType.emplace_back(VTK_WEDGE); + break; + case 3: + addCell(cellArray, topoElem); + vtkType.emplace_back(VTK_TETRA); + break; + case 4: + addCell(cellArray, topoElem); + vtkType.emplace_back(VTK_QUADRATIC_HEXAHEDRON); + break; + case 5: + addCell(cellArray, topoElem); + vtkType.emplace_back(VTK_QUADRATIC_WEDGE); + break; + case 6: + addCell(cellArray, topoElem); + vtkType.emplace_back(VTK_QUADRATIC_TETRA); + break; + case 7: + addCell(cellArray, topoElem); + vtkType.emplace_back(VTK_TRIANGLE); + break; + case 8: + addCell(cellArray, topoElem); + vtkType.emplace_back(VTK_QUADRATIC_TRIANGLE); + break; + case 9: + addCell(cellArray, topoElem); + vtkType.emplace_back(VTK_QUAD); + break; + case 10: + addCell(cellArray, topoElem); + vtkType.emplace_back(VTK_QUADRATIC_QUAD); + break; + case 11: + addCell(cellArray, topoElem); + vtkType.emplace_back(VTK_LINE); + break; + case 12: + addCell(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 identifyScalarEntities(const std::vector> entities) +{ + std::vector 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 +readNodes(std::ifstream& ifstr, const std::string& lines, vtkSmartPointer& 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 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 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 readElements(std::ifstream& ifstr, + const std::string& lines, + const std::map& mapNodes, + vtkSmartPointer& 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 info(3); + std::map mapElem; + std::vector topoElem; + std::vector 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::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& mapNodes, + const FRDResultInfo& info, + vtkSmartPointer& 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 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> 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 et = {0, 0, 0, 0}; + // fill entityType, ignore MENU: " 1" + sub = sub.substr(8 + 5, 4 * 5); + std::string_view::iterator it1; + std::vector::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 vecValues; + std::vector scaValues; + std::vector 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 vecArray = vtkSmartPointer::New(); + // arrays for scalar entities (if needed) + std::vector> scaArrays; + for (size_t i = 0; i < scalarPos.size(); ++i) { + scaArrays.emplace_back(vtkSmartPointer::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>::iterator it1; + std::vector::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 readFRD(std::ifstream& ifstr) +{ + auto points = vtkSmartPointer::New(); + auto cells = vtkSmartPointer::New(); + auto multiBlock = vtkSmartPointer::New(); + vtkSmartPointer grid; + std::map, vtkSmartPointer> grids; + std::string line; + std::map mapNodes; + std::vector 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(info.step, info.analysisType)); + if (it == grids.end()) { + grid = vtkSmartPointer::New(); + grid->SetPoints(points); + grid->SetCells(cellTypes.data(), cells); + grids[std::pair(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::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 multiBlock = FRDReader::readFRD(ifstr); + + std::string dir = fi.dirPath(); + + auto writer = vtkSmartPointer::New(); + + std::string blockFile = dir + "/" + fi.fileNamePure() + "." + writer->GetDefaultFileExtension(); + writer->SetFileName(blockFile.c_str()); + writer->SetInputData(multiBlock); + writer->Update(); +} + } // namespace Fem diff --git a/src/Mod/Fem/App/FemVTKTools.h b/src/Mod/Fem/App/FemVTKTools.h index 34a3552e76..00db4d0763 100644 --- a/src/Mod/Fem/App/FemVTKTools.h +++ b/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