From e84eb920dcd396d4daafc5a9501dfd1a93d557c1 Mon Sep 17 00:00:00 2001 From: marioalexis Date: Wed, 29 Jan 2025 19:18:29 -0300 Subject: [PATCH 1/2] 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 From 5042b944630b3bb64a7e1ba4c84cc1bf11868f2c Mon Sep 17 00:00:00 2001 From: marioalexis Date: Mon, 17 Feb 2025 08:26:03 -0300 Subject: [PATCH 2/2] Fem: Add metadata info to blocks --- src/Mod/Fem/App/FemVTKTools.cpp | 233 ++++++++++++++++++++++++-------- 1 file changed, 179 insertions(+), 54 deletions(-) diff --git a/src/Mod/Fem/App/FemVTKTools.cpp b/src/Mod/Fem/App/FemVTKTools.cpp index 522b461e0c..382ff07455 100644 --- a/src/Mod/Fem/App/FemVTKTools.cpp +++ b/src/Mod/Fem/App/FemVTKTools.cpp @@ -39,6 +39,7 @@ #include #include #include +#include #include #include #include @@ -53,6 +54,7 @@ #include #include #include +#include #include #include #include @@ -1036,20 +1038,60 @@ void FemVTKTools::exportFreeCADResult(const App::DocumentObject* result, namespace FRDReader { +enum class ElementType +{ + Edge = 11, + QuadEdge = 12, + Triangle = 7, + QuadTriangle = 8, + Quadrangle = 9, + QuadQuadrangle = 10, + Tetra = 3, + QuadTetra = 6, + Hexa = 1, + QuadHexa = 4, + Penta = 2, + QuadPenta = 5 +}; + +enum class AnalysisType +{ + Static = 0, + TimeStep = 1, + Frequency = 2, + LoadStep = 3, + UserNamed = 4 +}; + +std::map mapAnalysisTypeToStr = {{AnalysisType::Static, "Static"}, + {AnalysisType::TimeStep, "TimeStep"}, + {AnalysisType::Frequency, "Frequency"}, + {AnalysisType::LoadStep, "LoadStep"}, + {AnalysisType::UserNamed, "User"}}; + +// value format indicator +enum class Indicator +{ + Short = 0, + Long = 1, + // BinaryFloat = 2, not used + // BinaryDouble = 3 not used +}; + // 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}, +std::map mapCcxTypeNodes = { + {ElementType::Edge, 2}, + {ElementType::QuadEdge, 3}, + {ElementType::Triangle, 3}, + {ElementType::QuadTriangle, 6}, + {ElementType::Quadrangle, 4}, + {ElementType::QuadQuadrangle, 8}, + {ElementType::Tetra, 4}, + {ElementType::QuadTetra, 10}, + {ElementType::Hexa, 8}, + {ElementType::QuadHexa, 20}, + {ElementType::Penta, 6}, + {ElementType::QuadPenta, 15}, }; // map CalculiX nodes order to Vtk order @@ -1118,54 +1160,54 @@ void addCell(vtkSmartPointer& cellArray, const std::vector& t void fillCell(vtkSmartPointer& cellArray, std::vector& topoElem, std::vector& vtkType, - int elemType) + ElementType elemType) { switch (elemType) { - case 1: + case ElementType::Hexa: addCell(cellArray, topoElem); vtkType.emplace_back(VTK_HEXAHEDRON); break; - case 2: + case ElementType::Penta: addCell(cellArray, topoElem); vtkType.emplace_back(VTK_WEDGE); break; - case 3: + case ElementType::Tetra: addCell(cellArray, topoElem); vtkType.emplace_back(VTK_TETRA); break; - case 4: + case ElementType::QuadHexa: addCell(cellArray, topoElem); vtkType.emplace_back(VTK_QUADRATIC_HEXAHEDRON); break; - case 5: + case ElementType::QuadPenta: addCell(cellArray, topoElem); vtkType.emplace_back(VTK_QUADRATIC_WEDGE); break; - case 6: + case ElementType::QuadTetra: addCell(cellArray, topoElem); vtkType.emplace_back(VTK_QUADRATIC_TETRA); break; - case 7: + case ElementType::Triangle: addCell(cellArray, topoElem); vtkType.emplace_back(VTK_TRIANGLE); break; - case 8: + case ElementType::QuadTriangle: addCell(cellArray, topoElem); vtkType.emplace_back(VTK_QUADRATIC_TRIANGLE); break; - case 9: + case ElementType::Quadrangle: addCell(cellArray, topoElem); vtkType.emplace_back(VTK_QUAD); break; - case 10: + case ElementType::QuadQuadrangle: addCell(cellArray, topoElem); vtkType.emplace_back(VTK_QUADRATIC_QUAD); break; - case 11: + case ElementType::Edge: addCell(cellArray, topoElem); vtkType.emplace_back(VTK_LINE); break; - case 12: + case ElementType::QuadEdge: addCell(cellArray, topoElem); vtkType.emplace_back(VTK_QUADRATIC_EDGE); break; @@ -1176,22 +1218,42 @@ struct FRDResultInfo { double value; long numNodes; - int analysisType; + AnalysisType analysisType; int step; - int indicator; + Indicator indicator; + + bool operator==(const FRDResultInfo& other) const + { + return (this->step == other.step) && (this->analysisType == other.analysisType); + } + bool operator<(const FRDResultInfo& other) const + { + if (this->step < other.step) { + return true; + } + else if (this->step > other.step) { + return false; + } + else if (static_cast(this->analysisType) < static_cast(other.analysisType)) { + return true; + } + else { + return false; + } + } }; // get number of digits from format indicator -int getDigits(int indicator) +int getDigits(Indicator indicator) { int digits = 0; - if (indicator == 0) { - // short - digits = 5; - } - else if (indicator == 1) { - // long - digits = 10; + switch (indicator) { + case Indicator::Short: + digits = 5; + break; + case Indicator::Long: + digits = 10; + break; } return digits; @@ -1232,7 +1294,7 @@ readNodes(std::ifstream& ifstr, const std::string& lines, vtkSmartPointer(indicator)); points->SetNumberOfPoints(numNodes); @@ -1286,7 +1348,7 @@ std::vector readElements(std::ifstream& ifstr, sub = sub.substr(12 + 37); valueFromLine(sub.begin(), 1, indicator); - int digits = getDigits(indicator); + int digits = getDigits(static_cast(indicator)); while (elemID < numElem && std::getline(ifstr, line)) { std::string_view view {line}; if (view.rfind(keyCodeType, 0) == 0) { @@ -1309,8 +1371,8 @@ std::vector readElements(std::ifstream& ifstr, } // add cell to cellArray - if (topoElem.size() == mapCcxTypeNodes[info[0]]) { - fillCell(cellArray, topoElem, vtkType, info[0]); + if (topoElem.size() == mapCcxTypeNodes[static_cast(info[0])]) { + fillCell(cellArray, topoElem, vtkType, static_cast(info[0])); topoElem.clear(); mapElem[elem] = elemID++; } @@ -1342,13 +1404,17 @@ void readResultInfo(std::ifstream& ifstr, const std::string& lines, FRDResultInf valueFromLine(sub.begin(), 12, info.numNodes); sub = sub.substr(12 + 20); - valueFromLine(sub.begin(), 2, info.analysisType); + int anType; + valueFromLine(sub.begin(), 2, anType); + info.analysisType = static_cast(anType); sub = sub.substr(2); valueFromLine(sub.begin(), 5, info.step); sub = sub.substr(5 + 10); - valueFromLine(sub.begin(), 2, info.indicator); + int ind; + valueFromLine(sub.begin(), 2, ind); + info.indicator = static_cast(ind); } // read result from nodal result block and add result array to grid @@ -1519,6 +1585,25 @@ void readResults(std::ifstream& ifstr, grid->GetPointData()->AddArray(s); } } +vtkSmartPointer createTimeInfo(const std::string& type) +{ + auto timeInfo = vtkSmartPointer::New(); + timeInfo->SetName("TimeInfo"); + timeInfo->InsertNextValue(type); + // set unit to empty string + timeInfo->InsertNextValue(""); + + return timeInfo; +} + +vtkSmartPointer createTimeValue(const double& value) +{ + auto stepValue = vtkSmartPointer::New(); + stepValue->SetName("TimeValue"); + stepValue->InsertNextValue(value); + + return stepValue; +} vtkSmartPointer readFRD(std::ifstream& ifstr) { @@ -1526,7 +1611,9 @@ vtkSmartPointer readFRD(std::ifstream& ifstr) auto cells = vtkSmartPointer::New(); auto multiBlock = vtkSmartPointer::New(); vtkSmartPointer grid; - std::map, vtkSmartPointer> grids; + vtkSmartPointer block; + std::map> grids; + std::map> blocks; std::string line; std::map mapNodes; std::vector cellTypes; @@ -1554,12 +1641,34 @@ vtkSmartPointer readFRD(std::ifstream& ifstr) // read result info block FRDResultInfo info; readResultInfo(ifstr, line, info); - auto it = grids.find(std::pair(info.step, info.analysisType)); + auto it = grids.find(info); if (it == grids.end()) { + // create TimeInfo metadata + auto timeInfo = createTimeInfo(mapAnalysisTypeToStr[info.analysisType]); + // search analysis type block and create it if necessary + auto it2 = blocks.find(info.analysisType); + if (it2 == blocks.end()) { + block = vtkSmartPointer::New(); + block->GetFieldData()->AddArray(timeInfo); + blocks[info.analysisType] = block; + } + else { + block = it2->second; + } + // create unstructured grid grid = vtkSmartPointer::New(); grid->SetPoints(points); grid->SetCells(cellTypes.data(), cells); - grids[std::pair(info.step, info.analysisType)] = grid; + + // create TimeValue metadata + auto stepValue = createTimeValue(info.value); + + grid->GetFieldData()->AddArray(stepValue); + grid->GetFieldData()->AddArray(timeInfo); + + grids[info] = grid; + unsigned int nb = block->GetNumberOfBlocks(); + block->SetBlock(nb, grid); } else { grid = (*it).second; @@ -1570,18 +1679,25 @@ vtkSmartPointer readFRD(std::ifstream& ifstr) } int i = 0; - multiBlock->SetNumberOfBlocks(grids.size()); - for (const auto& g : grids) { - multiBlock->SetBlock(i, g.second); + for (const auto& b : blocks) { + multiBlock->SetBlock(i, b.second); ++i; } // save points and elements even without results if (grids.empty()) { + block = vtkSmartPointer::New(); grid = vtkSmartPointer::New(); grid->SetPoints(points); grid->SetCells(cellTypes.data(), cells); - multiBlock->SetBlock(0, grid); + auto timeInfo = createTimeInfo(""); + auto stepValue = createTimeValue(0); + grid->GetFieldData()->AddArray(stepValue); + grid->GetFieldData()->AddArray(timeInfo); + + block->SetBlock(0, grid); + block->GetFieldData()->AddArray(timeInfo); + multiBlock->SetBlock(0, block); } return multiBlock; @@ -1603,12 +1719,21 @@ void FemVTKTools::frdToVTK(const char* filename) std::string dir = fi.dirPath(); - auto writer = vtkSmartPointer::New(); + for (unsigned int i = 0; i < multiBlock->GetNumberOfBlocks(); ++i) { + vtkDataObject* block = multiBlock->GetBlock(i); + // get TimeInfo + vtkSmartPointer info = + vtkStringArray::SafeDownCast(block->GetFieldData()->GetAbstractArray(0)); + std::string type = info->GetValue(0).c_str(); - std::string blockFile = dir + "/" + fi.fileNamePure() + "." + writer->GetDefaultFileExtension(); - writer->SetFileName(blockFile.c_str()); - writer->SetInputData(multiBlock); - writer->Update(); + auto writer = vtkSmartPointer::New(); + + std::string blockFile = + dir + "/" + fi.fileNamePure() + type + "." + writer->GetDefaultFileExtension(); + writer->SetFileName(blockFile.c_str()); + writer->SetInputData(block); + writer->Update(); + } } } // namespace Fem