diff --git a/src/Mod/Fem/App/AppFemPy.cpp b/src/Mod/Fem/App/AppFemPy.cpp index 1764428cf469..ac8441248d75 100644 --- a/src/Mod/Fem/App/AppFemPy.cpp +++ b/src/Mod/Fem/App/AppFemPy.cpp @@ -67,6 +67,7 @@ #include "FemMeshPy.h" #ifdef FC_USE_VTK #include "FemPostPipeline.h" +#include "FemVTKTools.h" #endif #include @@ -93,6 +94,11 @@ class Module : public Py::ExtensionModule add_varargs_method("read",&Module::read, "Read a mesh from a file and returns a Mesh object." ); +#ifdef FC_USE_VTK + add_varargs_method("readCfdResult",&Module::readCfdResult, + "Read a CFD result from a file and returns a Result object." + ); +#endif add_varargs_method("show",&Module::show, "show(shape) -- Add the shape to the active document or create one if no document exists." ); @@ -236,6 +242,30 @@ class Module : public Py::ExtensionModule mesh->read(EncodedName.c_str()); return Py::asObject(new FemMeshPy(mesh.release())); } + +#ifdef FC_USE_VTK + Py::Object readCfdResult(const Py::Tuple& args) + { + PyObject *pcObj; + char* Name; // PythonFeatureT is of type `App::DocumentObjectPy` + if (!PyArg_ParseTuple(args.ptr(), "etO","utf-8", &Name, &(App::DocumentObjectPy::Type), &pcObj)) + throw Py::Exception(); + std::string EncodedName = std::string(Name); + PyMem_Free(Name); + // this function needs the second parameter: App::DocumentObjectPy, since it is created in python + if (pcObj) + { + //App::DocumentObjectPy objpy(pcObj); + App::DocumentObject* obj = static_cast(pcObj)->getDocumentObjectPtr(); + FemVTKTools::readFluidicResult(EncodedName.c_str(), obj); + } + else + FemVTKTools::readFluidicResult(EncodedName.c_str()); + + return Py::None(); + } +#endif + Py::Object show(const Py::Tuple& args) { PyObject *pcObj; diff --git a/src/Mod/Fem/App/CMakeLists.txt b/src/Mod/Fem/App/CMakeLists.txt index 25722b5d6500..7afb6fd1a39b 100755 --- a/src/Mod/Fem/App/CMakeLists.txt +++ b/src/Mod/Fem/App/CMakeLists.txt @@ -240,6 +240,8 @@ if(BUILD_FEM_VTK) FemPostFilter.cpp FemPostFunction.h FemPostFunction.cpp + FemVTKTools.h + FemVTKTools.cpp ) SOURCE_GROUP("PostObjects" FILES ${FemPost_SRCS}) endif(BUILD_FEM_VTK) diff --git a/src/Mod/Fem/App/FemMesh.cpp b/src/Mod/Fem/App/FemMesh.cpp index 716e89d7233e..5ea79820801e 100644 --- a/src/Mod/Fem/App/FemMesh.cpp +++ b/src/Mod/Fem/App/FemMesh.cpp @@ -49,6 +49,9 @@ #include #include "FemMesh.h" +#ifdef FC_USE_VTK + #include "FemVTKTools.h" +#endif #include #include @@ -898,6 +901,12 @@ void FemMesh::read(const char *FileName) // read Nastran-file readNastran(File.filePath()); } +#ifdef FC_USE_VTK + else if (File.hasExtension("vtk") || File.hasExtension("vtu")) { + // read *.vtk legacy format or *.vtu XML unstructure Mesh + FemVTKTools::readVTKMesh(File.filePath().c_str(), this); + } +#endif else{ throw Base::Exception("Unknown extension"); } @@ -1185,6 +1194,12 @@ void FemMesh::write(const char *FileName) const // write ABAQUS Output writeABAQUS(File.filePath()); } +#ifdef FC_USE_VTK + else if (File.hasExtension("vtk") || File.hasExtension("vtu") ) { + // write unstructure mesh to VTK format *.vtk and *.vtu + FemVTKTools::writeVTKMesh(File.filePath().c_str(), this); + } +#endif else{ throw Base::Exception("Unknown extension"); } diff --git a/src/Mod/Fem/App/FemVTKTools.cpp b/src/Mod/Fem/App/FemVTKTools.cpp new file mode 100644 index 000000000000..6d1f57907ee1 --- /dev/null +++ b/src/Mod/Fem/App/FemVTKTools.cpp @@ -0,0 +1,745 @@ +/*************************************************************************** + * Copyright (c) Jürgen Riegel (juergen.riegel@web.de) 2009 * + * * + * This file is part of the FreeCAD CAx development system. * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Library General Public * + * License as published by the Free Software Foundation; either * + * version 2 of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU Library General Public License for more details. * + * * + * You should have received a copy of the GNU Library General Public * + * License along with this library; see the file COPYING.LIB. If not, * + * write to the Free Software Foundation, Inc., 59 Temple Place, * + * Suite 330, Boston, MA 02111-1307, USA * + * * + ***************************************************************************/ + + +#include "PreCompiled.h" + +#ifndef _PreComp_ +# include +# include +# include + +# include +# include +# include +# include +# include +# include +# include +#endif + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +# include +# include +# include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "FemVTKTools.h" +#include "FemMeshProperty.h" +#include "FemAnalysis.h" + +namespace Fem +{ + +template vtkDataSet* readVTKFile(const char*fileName) +{ + vtkSmartPointer reader = + vtkSmartPointer::New(); + reader->SetFileName(fileName); + reader->Update(); + reader->GetOutput()->Register(reader); + return vtkDataSet::SafeDownCast(reader->GetOutput()); +} + +template void writeVTKFile(const char* filename, vtkSmartPointer dataset) +{ + vtkSmartPointer writer = + vtkSmartPointer::New(); + writer->SetFileName(filename); + writer->SetInputData(dataset); + writer->Write(); + //writer->Update(); +} + +void FemVTKTools::importVTKMesh(vtkSmartPointer dataset, FemMesh* mesh) +{ + const vtkIdType nPoints = dataset->GetNumberOfPoints(); + const vtkIdType nCells = dataset->GetNumberOfCells(); + Base::Console().Log("%d nodes/points and %d cells/elements found!\n", nPoints, nCells); + + //vtkSmartPointer cells = dataset->GetCells(); // works only for vtkUnstructuredGrid + vtkSmartPointer idlist= vtkSmartPointer::New(); + + //Now fill the SMESH datastructure + SMESH_Mesh* smesh = const_cast(mesh->getSMesh()); + SMESHDS_Mesh* meshds = smesh->GetMeshDS(); + meshds->ClearMesh(); + + for(vtkIdType i=0; iGetPoint(i); + meshds->AddNodeWithID(p[0], p[1], p[2], i+1); + } + + for(vtkIdType iCell=0; iCellReset(); + idlist = dataset->GetCell(iCell)->GetPointIds(); + vtkIdType *ids = idlist->GetPointer(0); + // 3D cells first + switch(dataset->GetCellType(iCell)) + { + case VTK_TETRA: + meshds->AddVolumeWithID(ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1, iCell+1); + break; + case VTK_HEXAHEDRON: + meshds->AddVolumeWithID(ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1, ids[4]+1, ids[5]+1, ids[6]+1, ids[7]+1, iCell+1); + break; + case VTK_QUADRATIC_TETRA: + meshds->AddVolumeWithID(ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1, ids[4]+1, ids[5]+1, ids[6]+1, ids[7]+1, ids[8]+1, ids[9]+1, iCell+1); + break; + case VTK_QUADRATIC_HEXAHEDRON: + meshds->AddVolumeWithID(ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1, ids[4]+1, ids[5]+1, ids[6]+1, ids[7]+1, ids[8]+1, ids[9]+1,\ + ids[10]+1, ids[11]+1, ids[12]+1, ids[13]+1, ids[14]+1, ids[15]+1, ids[16]+1, ids[17]+1, ids[18]+1, ids[19]+1,\ + iCell+1); + break; + case VTK_WEDGE: + meshds->AddVolumeWithID(ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1, ids[4]+1, ids[5]+1, iCell+1); + break; + case VTK_PYRAMID: + meshds->AddVolumeWithID(ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1, ids[4]+1, iCell+1); + break; + // 2D elements + case VTK_TRIANGLE: + meshds->AddFaceWithID(ids[0]+1, ids[1]+1, ids[2]+1, iCell+1); + break; + case VTK_QUADRATIC_TRIANGLE: + meshds->AddFaceWithID(ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1, ids[4]+1, ids[5]+1, iCell+1); + break; + case VTK_QUAD: + meshds->AddFaceWithID(ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1, iCell+1); + break; + case VTK_QUADRATIC_QUAD: + meshds->AddFaceWithID(ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1, ids[4]+1, ids[5]+1, ids[6]+1, ids[7]+1, iCell+1); + break; + default: + { + Base::Console().Error("Only common 2D and 3D Cells are supported in VTK mesh import\n"); + break; + } + } + } +} + +FemMesh* FemVTKTools::readVTKMesh(const char* filename, FemMesh* mesh) +{ + Base::TimeInfo Start; + Base::Console().Log("Start: read FemMesh from VTK unstructuredGrid ======================\n"); + Base::FileInfo f(filename); + + if(f.hasExtension("vtu")) + { + vtkSmartPointer dataset = readVTKFile(filename); + importVTKMesh(dataset, mesh); + } + else if(f.hasExtension("vtk")) + { + vtkSmartPointer dataset = readVTKFile(filename); + importVTKMesh(dataset, mesh); + } + else + { + Base::Console().Error("file name extension is not supported\n"); + return NULL; + } + + Base::Console().Log(" %f: Done \n",Base::TimeInfo::diffTimeF(Start,Base::TimeInfo())); + return mesh; +} + +void exportFemMeshFaces(vtkSmartPointer grid, const SMDS_FaceIteratorPtr& aFaceIter) +{ + + vtkSmartPointer triangleArray = vtkSmartPointer::New(); + vtkSmartPointer quadTriangleArray = vtkSmartPointer::New(); + vtkSmartPointer quadArray = vtkSmartPointer::New(); + vtkSmartPointer quadQuadArray = vtkSmartPointer::New(); + + for (;aFaceIter->more();) + { + const SMDS_MeshFace* aFace = aFaceIter->next(); + + //triangle + if(aFace->NbNodes() == 3) + { + vtkSmartPointer tria = vtkSmartPointer::New(); + tria->GetPointIds()->SetId(0, aFace->GetNode(0)->GetID()-1); + tria->GetPointIds()->SetId(1, aFace->GetNode(1)->GetID()-1); + tria->GetPointIds()->SetId(2, aFace->GetNode(2)->GetID()-1); + + triangleArray->InsertNextCell(tria); + } + //quad + else if(aFace->NbNodes() == 4) + { + vtkSmartPointer quad = vtkSmartPointer::New(); + quad->GetPointIds()->SetId(0, aFace->GetNode(0)->GetID()-1); + quad->GetPointIds()->SetId(1, aFace->GetNode(1)->GetID()-1); + quad->GetPointIds()->SetId(2, aFace->GetNode(2)->GetID()-1); + quad->GetPointIds()->SetId(3, aFace->GetNode(3)->GetID()-1); + + quadArray->InsertNextCell(quad); + } + //quadratic triangle + else if (aFace->NbNodes() == 6) + { + vtkSmartPointer tria = vtkSmartPointer::New(); + tria->GetPointIds()->SetId(0, aFace->GetNode(0)->GetID()-1); + tria->GetPointIds()->SetId(1, aFace->GetNode(1)->GetID()-1); + tria->GetPointIds()->SetId(2, aFace->GetNode(2)->GetID()-1); + tria->GetPointIds()->SetId(3, aFace->GetNode(3)->GetID()-1); + tria->GetPointIds()->SetId(4, aFace->GetNode(4)->GetID()-1); + tria->GetPointIds()->SetId(5, aFace->GetNode(5)->GetID()-1); + quadTriangleArray->InsertNextCell(tria); + } + //quadratic quad + else if(aFace->NbNodes() == 8) + { + vtkSmartPointer quad = vtkSmartPointer::New(); + quad->GetPointIds()->SetId(0, aFace->GetNode(0)->GetID()-1); + quad->GetPointIds()->SetId(1, aFace->GetNode(1)->GetID()-1); + quad->GetPointIds()->SetId(2, aFace->GetNode(2)->GetID()-1); + quad->GetPointIds()->SetId(3, aFace->GetNode(3)->GetID()-1); + quad->GetPointIds()->SetId(4, aFace->GetNode(4)->GetID()-1); + quad->GetPointIds()->SetId(5, aFace->GetNode(5)->GetID()-1); + quad->GetPointIds()->SetId(6, aFace->GetNode(6)->GetID()-1); + quad->GetPointIds()->SetId(7, aFace->GetNode(7)->GetID()-1); + + quadQuadArray->InsertNextCell(quad); + } + } + if(triangleArray->GetNumberOfCells()>0) + grid->SetCells(VTK_TRIANGLE, triangleArray); + + if(quadArray->GetNumberOfCells()>0) + grid->SetCells(VTK_QUAD, quadArray); + + if(quadTriangleArray->GetNumberOfCells()>0) + grid->SetCells(VTK_QUADRATIC_TRIANGLE, quadTriangleArray); + + if(quadQuadArray->GetNumberOfCells()>0) + grid->SetCells(VTK_QUADRATIC_QUAD, quadQuadArray); + +} + +void exportFemMeshCells(vtkSmartPointer grid, const SMDS_VolumeIteratorPtr& aVolIter) +{ + // add common CFD cells like hex, wedge, prism, in addition to tetra + vtkSmartPointer tetraArray = vtkSmartPointer::New(); + vtkSmartPointer hexaArray = vtkSmartPointer::New(); + vtkSmartPointer wedgeArray = vtkSmartPointer::New(); + vtkSmartPointer pyramidArray = vtkSmartPointer::New(); + // quadratic elemnts with 13 and 15 nodes are not added yet + vtkSmartPointer quadTetraArray = vtkSmartPointer::New(); + vtkSmartPointer quadHexaArray = vtkSmartPointer::New(); + + for (;aVolIter->more();) + { + const SMDS_MeshVolume* aVol = aVolIter->next(); + + //tetrahedra + if(aVol->NbNodes() == 4) { + vtkSmartPointer tetra = vtkSmartPointer::New(); + tetra->GetPointIds()->SetId(0, aVol->GetNode(0)->GetID()-1); + tetra->GetPointIds()->SetId(1, aVol->GetNode(1)->GetID()-1); + tetra->GetPointIds()->SetId(2, aVol->GetNode(2)->GetID()-1); + tetra->GetPointIds()->SetId(3, aVol->GetNode(3)->GetID()-1); + + tetraArray->InsertNextCell(tetra); + } + // common cell types for CFD + if(aVol->NbNodes() == 5) { + vtkSmartPointer cell= vtkSmartPointer::New(); + cell->GetPointIds()->SetId(0, aVol->GetNode(0)->GetID()-1); + cell->GetPointIds()->SetId(1, aVol->GetNode(1)->GetID()-1); + cell->GetPointIds()->SetId(2, aVol->GetNode(2)->GetID()-1); + cell->GetPointIds()->SetId(3, aVol->GetNode(3)->GetID()-1); + cell->GetPointIds()->SetId(4, aVol->GetNode(4)->GetID()-1); + + pyramidArray->InsertNextCell(cell); + } + if(aVol->NbNodes() == 6) { + vtkSmartPointer cell = vtkSmartPointer::New(); + cell->GetPointIds()->SetId(0, aVol->GetNode(0)->GetID()-1); + cell->GetPointIds()->SetId(1, aVol->GetNode(1)->GetID()-1); + cell->GetPointIds()->SetId(2, aVol->GetNode(2)->GetID()-1); + cell->GetPointIds()->SetId(3, aVol->GetNode(3)->GetID()-1); + cell->GetPointIds()->SetId(4, aVol->GetNode(4)->GetID()-1); + cell->GetPointIds()->SetId(5, aVol->GetNode(5)->GetID()-1); + + wedgeArray->InsertNextCell(cell); + } + if(aVol->NbNodes() == 8) { + vtkSmartPointer cell= vtkSmartPointer::New(); + cell->GetPointIds()->SetId(0, aVol->GetNode(0)->GetID()-1); + cell->GetPointIds()->SetId(1, aVol->GetNode(1)->GetID()-1); + cell->GetPointIds()->SetId(2, aVol->GetNode(2)->GetID()-1); + cell->GetPointIds()->SetId(3, aVol->GetNode(3)->GetID()-1); + cell->GetPointIds()->SetId(4, aVol->GetNode(4)->GetID()-1); + cell->GetPointIds()->SetId(5, aVol->GetNode(5)->GetID()-1); + cell->GetPointIds()->SetId(6, aVol->GetNode(6)->GetID()-1); + cell->GetPointIds()->SetId(7, aVol->GetNode(7)->GetID()-1); + + hexaArray->InsertNextCell(cell); + } + //quadratic tetrahedra + else if( aVol->NbNodes() == 10) { + + vtkSmartPointer tetra = vtkSmartPointer::New(); + for(int i=0; i<10; i++){ + tetra->GetPointIds()->SetId(i, aVol->GetNode(i)->GetID()-1); + } + + quadTetraArray->InsertNextCell(tetra); + } + if(aVol->NbNodes() == 20) { // not tested, no sure about vertex sequence + vtkSmartPointer cell= vtkSmartPointer::New(); + for(int i=0; i<20; i++){ + cell->GetPointIds()->SetId(i, aVol->GetNode(i)->GetID()-1); + } + hexaArray->InsertNextCell(cell); + } + } + + if(tetraArray->GetNumberOfCells()>0) + grid->SetCells(VTK_TETRA, tetraArray); + + if(pyramidArray->GetNumberOfCells()>0) + grid->SetCells(VTK_PYRAMID, pyramidArray); + + if(wedgeArray->GetNumberOfCells()>0) + grid->SetCells(VTK_WEDGE, wedgeArray); + + if(hexaArray->GetNumberOfCells()>0) + grid->SetCells(VTK_HEXAHEDRON, hexaArray); + + if(quadTetraArray->GetNumberOfCells()>0) + grid->SetCells(VTK_QUADRATIC_TETRA, quadTetraArray); + + if(quadHexaArray->GetNumberOfCells()>0) + grid->SetCells(VTK_QUADRATIC_HEXAHEDRON, quadHexaArray); + +} + +void FemVTKTools::exportVTKMesh(const FemMesh* mesh, vtkSmartPointer grid) +{ + + SMESH_Mesh* smesh = const_cast(mesh->getSMesh()); + SMESHDS_Mesh* meshDS = smesh->GetMeshDS(); + const SMDS_MeshInfo& info = meshDS->GetMeshInfo(); + + //start with the nodes + vtkSmartPointer points = vtkSmartPointer::New(); + SMDS_NodeIteratorPtr aNodeIter = meshDS->nodesIterator(); + + points->SetNumberOfPoints(info.NbNodes()); + for(; aNodeIter->more(); ) { + const SMDS_MeshNode* node = aNodeIter->next(); // why float, not double? + double coords[3] = {double(node->X()), double(node->Y()), double(node->Z())}; + points->SetPoint(node->GetID()-1, coords); + } + grid->SetPoints(points); + //start with 2d elements + SMDS_FaceIteratorPtr aFaceIter = meshDS->facesIterator(); + exportFemMeshFaces(grid, aFaceIter); + //3D volume elements + SMDS_VolumeIteratorPtr aVolIter = meshDS->volumesIterator(); + exportFemMeshCells(grid, aVolIter); + +} + +void FemVTKTools::writeVTKMesh(const char* filename, const FemMesh* mesh) +{ + + Base::TimeInfo Start; + Base::Console().Log("Start: write FemMesh from VTK unstructuredGrid ======================\n"); + Base::FileInfo f(filename); + + vtkSmartPointer grid = vtkSmartPointer::New(); + exportVTKMesh(mesh, grid); + //vtkSmartPointer dataset = vtkDataSet::SafeDownCast(grid); + if(f.hasExtension("vtu")){ + writeVTKFile(filename, grid); + } + else if(f.hasExtension("vtk")){ + writeVTKFile(filename, grid); + } + else{ + Base::Console().Error("file name extension is not supported to write VTK\n"); + } + + Base::Console().Log(" %f: Done \n",Base::TimeInfo::diffTimeF(Start, Base::TimeInfo())); +} + + +template App::DocumentObject* getActiveFemObject(bool creating = false) +{ + App::Document* pcDoc = App::GetApplication().getActiveDocument(); + if(!pcDoc) + { + Base::Console().Message("No active document is found\n"); + pcDoc = App::GetApplication().newDocument(); + } + + App::DocumentObject* obj = pcDoc->getActiveObject(); + FemObjectT tobj; // will be added to document? will it destruct out of scope? + //Base::Type meshId = Base::Type::fromName("Fem::FemMeshObject"); + if(obj->getTypeId() == tobj.getTypeId()) + { + return obj; + } + else if (creating) + { + if(obj->getTypeId() == FemAnalysis::getClassTypeId()) + { + std::vector fem = (static_cast(obj))->Member.getValues(); + for (std::vector::iterator it = fem.begin(); it != fem.end(); ++it) { + if ((*it)->getTypeId().isDerivedFrom(tobj.getClassTypeId())) + return static_cast(*it); // return the first of that type + } + App::DocumentObject* newobj = pcDoc->addObject(tobj.getTypeId().getName()); + fem.push_back(newobj); // FemAnalysis is not a DocumentGroup derived class but DocumentObject + (static_cast(obj))->Member.setValues(fem); + return newobj; + } + else + { + return pcDoc->addObject(tobj.getTypeId().getName()); // create in the acitive document + } + } + return NULL; +} + + +App::DocumentObject* FemVTKTools::readFluidicResult(const char* filename, App::DocumentObject* res) +{ + Base::TimeInfo Start; + Base::Console().Log("Start: read FemResult with FemMesh from VTK file ======================\n"); + Base::FileInfo f(filename); + + vtkSmartPointer ds; + //vtkDataSet* ds; // + if(f.hasExtension("vtu")) + { + ds = readVTKFile(filename); + } + else if(f.hasExtension("vtk")) + { + ds = readVTKFile(filename); + } + else + { + Base::Console().Error("file name extension is not supported\n"); + } + + vtkSmartPointer dataset = ds; + App::DocumentObject* result = NULL; + if(!res) + result = static_cast(res); + else + //result = getActiveFemObject(); + Base::Console().Error("FemResultObject pointer is invalid\n"); + + App::DocumentObject* mesh = getActiveFemObject(true); + FemMesh* fmesh = new FemMesh(); // PropertyFemMesh instance is responsible to relase FemMesh ?? + static_cast(mesh->getPropertyByName("FemMesh"))->setValue(*fmesh); + importVTKMesh(dataset, fmesh); + static_cast(result->getPropertyByName("Mesh"))->setValue(mesh); + //PropertyLink is the property type to store DocumentObject pointer + + importFluidicResult(dataset, result); + App::Document *pcDoc = App::GetApplication().getActiveDocument(); + pcDoc->recompute(); + + Base::Console().Log(" %f: Done \n", Base::TimeInfo::diffTimeF(Start, Base::TimeInfo())); + + return result; +} + +void FemVTKTools::importFluidicResult(vtkSmartPointer dataset, App::DocumentObject* res) { + + // Temperature is optional, same for other turbulence related + std::map vars; // varable name openfoam -> defined in CfdResult.py + vars["Pressure"] = "p"; + vars["Velocity"] = "U"; + vars["Temperature"] = "T"; + vars["TurbulenceThermalDiffusivity"] = "alphat"; + vars["TurbulenceViscosity"] = "nut"; + vars["TurbulenceEnergy"] = "k"; + vars["TurbulenceDissipationRate"] = "epsilon"; + vars["TurbulenceSpecificDissipation"] = "omega"; + + const int max_var_index = 11; + std::vector stats(3*max_var_index); + for(int i=0; i<3*max_var_index; i++) + stats[i] = 0.0; + + std::map varids; + varids["Ux"] = 0; + varids["Uy"] = 1; + varids["Uz"] = 2; + varids["Umag"] = 3; + varids["Pressure"] = 4; + varids["Temperature"] = 5; + varids["TurbulenceThermalDiffusivity"] = 6; + varids["TurbulenceViscosity"] = 7; + varids["TurbulenceEnergy"] = 8; + varids["TurbulenceDissipationRate"] = 9; + varids["TurbulenceSpecificDissipation"] = 10; + + double ts = 0.0; // t=0.0 for static simulation + static_cast(res->getPropertyByName("Time"))->setValue(ts); + + vtkSmartPointer pd = dataset->GetPointData(); + const vtkIdType nPoints = dataset->GetNumberOfPoints(); + + double vmin=1.0e100, vmean=0.0, vmax=0.0; + vtkSmartPointer vel = pd->GetArray(vars["Velocity"]); + //vtkSmartPointer vel = vtkDoubleArray::SafeDownCast(pd->GetArray("Velocity")); + if(nPoints && vel && vel->GetNumberOfComponents() == 3) { + std::vector vec(nPoints); + for(vtkIdType i=0; iGetTuple(i); // GetTuple3(): only for vtkDataArray + double vmag = std::sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); + vmean += vmag; + if(vmag > vmax) vmax = vmag; + if(vmag < vmin) vmin = vmag; + vec.push_back(Base::Vector3d(p[0], p[1], p[2])); + } + int index = varids["Vmag"]; + stats[index*3] = vmin; + stats[index*3 + 2] = vmax; + stats[index*3 + 1] = vmean/nPoints; + App::PropertyVectorList* velocity = static_cast(res->getPropertyByName("Velocity")); + velocity->setValues(vec); // is that possible to avoid copy, using std::move() ? + } + else{ + Base::Console().Message("Error: pressure or velocity fields is not found in Cfd Result vtk file \n"); + } + + for(auto const& kv: vars){ + vtkSmartPointer vec = vtkDoubleArray::SafeDownCast(pd->GetArray(kv.second)); + if(!vec) { + App::PropertyFloatList* field = static_cast(res->getPropertyByName(kv.first)); + int index = varids[kv.first]; + stats[index*3] = vec->GetDataTypeValueMin(); + stats[index*3 + 2] = vec->GetDataTypeValueMin(); + stats[index*3 + 1] = (stats[index*3 + 2] + stats[index*3] ) * 0.5; // not mean value, but the middle range + for(vtkIdType i=0; iGetTuple1(); // both PropertyFloatList and vtkArray support [] operator, + field->set1Value(i, vec->GetTuple1(i)); + } + } + } + static_cast(res->getPropertyByName("Stats"))->setValues(stats); +} + +/* +void FemVTKTools::importMechanicalResult(const vtkDataSet* grid, App::DocumentObject* res) { + // to be implemented later by FemWorkbench developer +} + * */ + +void FemVTKTools::exportFluidicResult(const App::DocumentObject* res, vtkSmartPointer grid) { + //FemResultObject* res = static_cast(obj); + // Property defined in Python is dynamicProperty, which can not be accessed by ->PropertyName in c++ + if(!res->getPropertyByName("Velocity")){ + Base::Console().Message("Warning: essential field like `velocity` is not found in CfdResult\n"); + return; + } + App::PropertyVectorList* velocity = static_cast(res->getPropertyByName("Velocity")); + const std::vector& vel = velocity->getValues(); + if(!vel.empty()) { + vtkSmartPointer data = vtkSmartPointer::New(); + data->SetNumberOfComponents(3); + data->SetName("Velocity"); + + for(std::vector::const_iterator it=vel.begin(); it!=vel.end(); ++it) { + double tuple[] = {it->x, it->y, it->z}; + data->InsertNextTuple(tuple); + } + + grid->GetPointData()->AddArray(data); + } + else{ + Base::Console().Message("Warning: essential fields pressure and velocity is empty in CfdResult\n"); + } + // Temperature is optional, same for other turbulence related + std::vector vars; // varable name is defined in CfdResult.py + vars.push_back("Pressure"); + vars.push_back("Temperature"); + vars.push_back("TurbulenceThermalDiffusivity"); + vars.push_back("TurbulenceViscosity"); + vars.push_back("TurbulenceEnergy"); + vars.push_back("TurbulenceDissipationRate"); + vars.push_back("TurbulenceSpecificDissipation"); + for(auto const& var: vars){ + App::PropertyFloatList* field; + if (res->getPropertyByName(var)) + field = static_cast(res->getPropertyByName(var)); + if(!field && !field->getValues().empty()) { + const std::vector& vec = field->getValues(); + vtkSmartPointer data = vtkSmartPointer::New(); + data->SetNumberOfValues(vec.size()); + data->SetName(var); + + for(size_t i=0; iSetValue(i, vec[i]); + + grid->GetPointData()->AddArray(data); + } + } +} + + +void FemVTKTools::exportMechanicalResult(const App::DocumentObject* obj, vtkSmartPointer grid) { + // code redundance can be avoided by property inspection, consider refactoring + const FemResultObject* res = static_cast(obj); + if(!res->StressValues.getValues().empty()) { + const std::vector& vec = res->StressValues.getValues(); + vtkSmartPointer data = vtkSmartPointer::New(); + data->SetNumberOfValues(vec.size()); + data->SetName("Von Mises stress"); + + for(size_t i=0; iSetValue(i, vec[i]); + + grid->GetPointData()->AddArray(data); + } + + if(!res->StressValues.getValues().empty()) { + const std::vector& vec = res->MaxShear.getValues(); + vtkSmartPointer data = vtkSmartPointer::New(); + data->SetNumberOfValues(vec.size()); + data->SetName("Max shear stress (Tresca)"); + + for(size_t i=0; iSetValue(i, vec[i]); + + grid->GetPointData()->AddArray(data); + } + + if(!res->StressValues.getValues().empty()) { + const std::vector& vec = res->PrincipalMax.getValues(); + vtkSmartPointer data = vtkSmartPointer::New(); + data->SetNumberOfValues(vec.size()); + data->SetName("Maximum Principal stress"); + + for(size_t i=0; iSetValue(i, vec[i]); + + grid->GetPointData()->AddArray(data); + } + + if(!res->StressValues.getValues().empty()) { + const std::vector& vec = res->PrincipalMin.getValues(); + vtkSmartPointer data = vtkSmartPointer::New(); + data->SetNumberOfValues(vec.size()); + data->SetName("Minimum Principal stress"); + + for(size_t i=0; iSetValue(i, vec[i]); + + grid->GetPointData()->AddArray(data); + } + + if(!res->StressValues.getValues().empty()) { + const std::vector& vec = res->Temperature.getValues(); + vtkSmartPointer data = vtkSmartPointer::New(); + data->SetNumberOfValues(vec.size()); + data->SetName("Temperature"); + + for(size_t i=0; iSetValue(i, vec[i]); + + grid->GetPointData()->AddArray(data); + } + + if(!res->StressValues.getValues().empty()) { + const std::vector& vec = res->UserDefined.getValues(); + vtkSmartPointer data = vtkSmartPointer::New(); + data->SetNumberOfValues(vec.size()); + data->SetName("User Defined Results"); + + for(size_t i=0; iSetValue(i, vec[i]); + + grid->GetPointData()->AddArray(data); + } + + + if(!res->StressValues.getValues().empty()) { + const std::vector& vec = res->DisplacementVectors.getValues(); + vtkSmartPointer data = vtkSmartPointer::New(); + data->SetNumberOfComponents(3); + data->SetName("Displacement"); + + for(std::vector::const_iterator it=vec.begin(); it!=vec.end(); ++it) { + double tuple[] = {it->x, it->y, it->z}; + data->InsertNextTuple(tuple); + } + + grid->GetPointData()->AddArray(data); + } +} + +} // namespace diff --git a/src/Mod/Fem/App/FemVTKTools.h b/src/Mod/Fem/App/FemVTKTools.h new file mode 100644 index 000000000000..33124eaacca1 --- /dev/null +++ b/src/Mod/Fem/App/FemVTKTools.h @@ -0,0 +1,75 @@ +/*************************************************************************** + * Copyright (c) 2015 Werner Mayer * + * * + * This file is part of the FreeCAD CAx development system. * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Library General Public * + * License as published by the Free Software Foundation; either * + * version 2 of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU Library General Public License for more details. * + * * + * You should have received a copy of the GNU Library General Public * + * License along with this library; see the file COPYING.LIB. If not, * + * write to the Free Software Foundation, Inc., 59 Temple Place, * + * Suite 330, Boston, MA 02111-1307, USA * + * * + ***************************************************************************/ + + +#ifndef FEM_VTK_TOOLS_H +#define FEM_VTK_TOOLS_H + +#include "FemMesh.h" +#include "FemMeshObject.h" +#include +#include "FemResultObject.h" + +#include +#include +#include + +#include + +namespace Fem +{ + /*! + utitly class to import/export read/write vtk mesh and result + */ + class AppFemExport FemVTKTools + { + public: + /*! + FemMesh import from vtkUnstructuredGrid instance + */ + static void importVTKMesh(vtkSmartPointer grid, FemMesh* mesh); + /*! + FemMesh read from vtkUnstructuredGrid data file + */ + static FemMesh* readVTKMesh(const char* filename, FemMesh* mesh); + /*! + FemMesh export to vtkUnstructuredGrid instance + */ + static void exportVTKMesh(const FemMesh* mesh, vtkSmartPointer grid); + /*! + FemMesh write to vtkUnstructuredGrid data file + */ + static void writeVTKMesh(const char* Filename, const FemMesh* mesh); + + /*! + * FemResult export and import from vtkUnstructuredGrid, + */ + static void importFluidicResult(vtkSmartPointer dataset, App::DocumentObject* res); + static void exportFluidicResult(const App::DocumentObject* res, vtkSmartPointer grid); + static void exportMechanicalResult(const App::DocumentObject* res, vtkSmartPointer grid); + + static App::DocumentObject* readFluidicResult(const char* Filename, App::DocumentObject* res=NULL); + }; +} + +#endif //FEM_VTK_TOOLS_H + \ No newline at end of file diff --git a/src/Mod/Fem/Init.py b/src/Mod/Fem/Init.py index be027c9707f0..0dcfbc573b94 100644 --- a/src/Mod/Fem/Init.py +++ b/src/Mod/Fem/Init.py @@ -33,6 +33,8 @@ FreeCAD.addImportType("FEM formats (*.unv *.med *.dat *.bdf)", "Fem") if("BUILD_FEM_VTK" in FreeCAD.__cmake__): FreeCAD.addImportType("FEM results (*.vtk *.vtp *.vts *.vtr *.vtu *.vti)", "Fem") + FreeCAD.addImportType("FEM CFD Unstructure Mesh (*.vtk *.vtu)", "Fem") + FreeCAD.addExportType("FEM CFD Unstructure Mesh (*.vtk *.vtu)", "Fem") FreeCAD.addExportType("FEM formats (*.unv *.med *.dat *.inp)", "Fem") FreeCAD.addImportType("CalculiX result (*.frd)", "ccxFrdReader")