kindred/create
1
1
Fork 0
Code Issues 30 Pull Requests Actions Packages Projects 9 Releases 2 Wiki Activity

FEM: VTK post processing, changes:

Browse Source 
- vtk result file reading
- some changes in vtk post processing pipe line
This commit is contained in:
qingfengxia
2017-02-16 07:54:16 +01:00
committed by Bernd Hahnebach
parent 99fa203f32
commit 1c18867d1e
5 changed files with 332 additions and 275 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
src/Mod/Fem/App/AppFemPy.cpp
View File

@@ -95,8 +95,8 @@ public:
"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 (file format detected from file suffix)"
add_varargs_method("readResult",&Module::readResult,
"Read a CFD or Mechanical result (auto detect) from a file (file format detected from file suffix)"
);
add_varargs_method("writeResult",&Module::writeResult,
"write a CFD or FEM result (auto detect) to a file (file format detected from file suffix)"
@@ -247,7 +247,7 @@ private:
}
#ifdef FC_USE_VTK
Py::Object readCfdResult(const Py::Tuple& args)
Py::Object readResult(const Py::Tuple& args)
{
char* fileName = NULL;
char* objName = NULL;
@@ -263,11 +263,11 @@ private:
{
App::Document* pcDoc = App::GetApplication().getActiveDocument();
App::DocumentObject* obj = pcDoc->getObject(resName.c_str());
FemVTKTools::readFluidicResult(EncodedName.c_str(), obj);
FemVTKTools::readResult(EncodedName.c_str(), obj);
}
else
FemVTKTools::readFluidicResult(EncodedName.c_str()); //assuming activeObject can hold Result
FemVTKTools::readResult(EncodedName.c_str()); //assuming activeObject can hold Result
return Py::None();
}
@@ -291,7 +291,7 @@ private:
}
else
FemVTKTools::writeResult(EncodedName.c_str());
return Py::None();
}
#endif

12
src/Mod/Fem/App/FemPostPipeline.cpp
View File

@@ -116,11 +116,12 @@ DocumentObjectExecReturn* FemPostPipeline::execute(void) {
bool FemPostPipeline::canRead(Base::FileInfo File) {
if (File.hasExtension("vtk") ||
// from FemResult only unstructural mesh is supported in femvtktoools.cpp
File.hasExtension("vtp") ||
File.hasExtension("vts") ||
File.hasExtension("vtr") ||
File.hasExtension("vtu") ||
File.hasExtension("vti"))
File.hasExtension("vti") ||
File.hasExtension("vtu"))
return true;
return false;
@@ -244,9 +245,10 @@ bool FemPostPipeline::holdsPostObject(FemPostObject* obj) {
}
void FemPostPipeline::load(FemResultObject* res) {
if(!res->Mesh.getValue() || !res->Mesh.getValue()->isDerivedFrom(Fem::FemMeshObject::getClassTypeId()))
if(!res->Mesh.getValue() || !res->Mesh.getValue()->isDerivedFrom(Fem::FemMeshObject::getClassTypeId())) {
Base::Console().Warning("Mesh of result object is empty or not derived from Fem::FemMeshObject\n");
return;
}
//first copy the mesh over
//########################
@@ -256,7 +258,7 @@ void FemPostPipeline::load(FemResultObject* res) {
//Now copy the point data over
//############################
if(res->getPropertyByName("Velocity")){
if(res->getPropertyByName("Velocity")){ // consider better way to detect result type, res->Type == "CfdResult"
FemVTKTools::exportFluidicResult(res, grid);
}
else{

551
src/Mod/Fem/App/FemVTKTools.cpp
View File

@@ -1,5 +1,6 @@
/***************************************************************************
* Copyright (c) Jürgen Riegel (juergen.riegel@web.de) 2009 *
* Copyright (c) Qingfeng Xia (qingfeng.xia at oxford uni) 2017 *
* *
* This file is part of the FreeCAD CAx development system. *
* *
@@ -27,6 +28,7 @@
# include <cstdlib>
# include <memory>
# include <cmath>
# include <map>
# include <Bnd_Box.hxx>
# include <BRep_Tool.hxx>
@@ -110,12 +112,7 @@ template<class TWriter> void writeVTKFile(const char* filename, vtkSmartPointer<
writer->Write();
}
/*
double scale = 1000;
p[0] = p[0]* scale; // scale back to mm
p[1] = p[1]* scale;
p[1] = p[1]* scale;
*/
void FemVTKTools::importVTKMesh(vtkSmartPointer<vtkDataSet> dataset, FemMesh* mesh, float scale)
{
const vtkIdType nPoints = dataset->GetNumberOfPoints();
@@ -387,7 +384,7 @@ void exportFemMeshCells(vtkSmartPointer<vtkUnstructuredGrid> grid, const SMDS_Vo
}
void FemVTKTools::exportVTKMesh(const FemMesh* mesh, vtkSmartPointer<vtkUnstructuredGrid> grid)
void FemVTKTools::exportVTKMesh(const FemMesh* mesh, vtkSmartPointer<vtkUnstructuredGrid> grid, float scale)
{
SMESH_Mesh* smesh = const_cast<SMESH_Mesh*>(mesh->getSMesh());
@@ -401,7 +398,7 @@ void FemVTKTools::exportVTKMesh(const FemMesh* mesh, vtkSmartPointer<vtkUnstruct
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())};
double coords[3] = {double(node->X()*scale), double(node->Y()*scale), double(node->Z()*scale)};
points->SetPoint(node->GetID()-1, coords);
}
grid->SetPoints(points);
@@ -489,7 +486,7 @@ App::DocumentObject* createObjectByType(const Base::Type type)
}
App::DocumentObject* FemVTKTools::readFluidicResult(const char* filename, App::DocumentObject* res)
App::DocumentObject* FemVTKTools::readResult(const char* filename, App::DocumentObject* res)
{
Base::TimeInfo Start;
Base::Console().Log("Start: read FemResult with FemMesh from VTK file ======================\n");
@@ -502,6 +499,7 @@ App::DocumentObject* FemVTKTools::readFluidicResult(const char* filename, App::D
{
scale = 1000.0; // convert from meter in length of CFD result file
}
vtkSmartPointer<vtkDataSet> ds;
if(f.hasExtension("vtu"))
{
@@ -531,7 +529,7 @@ App::DocumentObject* FemVTKTools::readFluidicResult(const char* filename, App::D
else
{
Base::Console().Log("FemResultObject pointer is NULL, trying to get the active object\n");
if (obj->getTypeId() == Base::Type::fromName("Fem::FemResultObjectPython"))
if(obj->getTypeId() == Base::Type::fromName("Fem::FemResultObject"))
result = obj;
else
{
@@ -546,10 +544,23 @@ App::DocumentObject* FemVTKTools::readFluidicResult(const char* filename, App::D
static_cast<PropertyFemMesh*>(mesh->getPropertyByName("FemMesh"))->setValue(*fmesh);
static_cast<App::PropertyLink*>(result->getPropertyByName("Mesh"))->setValue(mesh);
// PropertyLink is the property type to store DocumentObject pointer
importFluidicResult(dataset, result);
vtkSmartPointer<vtkPointData> pd = dataset->GetPointData();
vtkSmartPointer<vtkDataArray> displ = pd->GetArray("Displacement"); // name in vtk file, not the property name
vtkSmartPointer<vtkDataArray> vel = pd->GetArray("U"); // name in vtk file, not the property name
if(vel)
{
importFluidicResult(dataset, result);
}
else if (displ)
{
importMechanicalResult(dataset, result);
}
else
{
Base::Console().Error("FemResult type, fluidic (array name of `U`) or mechanical (array name of `Displacement`) can not be detected\n");
}
pcDoc->recompute();
Base::Console().Log(" %f: Done \n", Base::TimeInfo::diffTimeF(Start, Base::TimeInfo()));
return result;
@@ -565,29 +576,38 @@ void FemVTKTools::writeResult(const char* filename, const App::DocumentObject* r
Base::Console().Message("No active document is found thus do nothing and return\n");
return;
}
res = pcDoc->getActiveObject(); //type checking
res = pcDoc->getActiveObject(); //type checking is done by caller
}
if(!res) {
Base::Console().Error("Result object pointer is invalid and it is not active oject");
return;
}
auto hGrp = App::GetApplication().GetParameterGroupByPath("User parameter:BaseApp/Preferences/Units");
int unitSchema = hGrp->GetInt("UserSchema",0);
float scale = 1.0;
if(unitSchema == 0) // standard mm
{
scale = 0.001; // convert from mm in FreeCAD to SI length in result file
}
Base::TimeInfo Start;
Base::Console().Log("Start: write FemResult or CfdResult to VTK unstructuredGrid dataset =======\n");
Base::Console().Message("Start: write FemResult or CfdResult to VTK unstructuredGrid dataset =======\n");
Base::FileInfo f(filename);
vtkSmartPointer<vtkUnstructuredGrid> grid = vtkSmartPointer<vtkUnstructuredGrid>::New();
App::DocumentObject* mesh = static_cast<App::PropertyLink*>(res->getPropertyByName("Mesh"))->getValue();
const FemMesh& fmesh = static_cast<PropertyFemMesh*>(mesh->getPropertyByName("FemMesh"))->getValue();
FemVTKTools::exportVTKMesh(&fmesh, grid);
FemVTKTools::exportVTKMesh(&fmesh, grid, scale);
if(res->getPropertyByName("Velocity")){
if(res->getPropertyByName("Velocity")){ // consider better way to detect result type, res->Type == "CfdResult"
FemVTKTools::exportFluidicResult(res, grid);
}
else if(res->getPropertyByName("DisplacementVectors")){
FemVTKTools::exportMechanicalResult(res, grid);
}
else{
Base::Console().Error("Result type can not be detected from unique property name like Velocity or DisplacementVectors\n");
return;
}
@@ -602,61 +622,19 @@ void FemVTKTools::writeResult(const char* filename, const App::DocumentObject* r
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()));
Base::Console().Message(" %f: result writing is done \n",Base::TimeInfo::diffTimeF(Start, Base::TimeInfo()));
}
void FemVTKTools::importFluidicResult(vtkSmartPointer<vtkDataSet> dataset, App::DocumentObject* res) {
// velocity and pressure are essential, Temperature is optional, so are turbulence related variables
std::map<const char*, const char*> vars; // varable name defined in openfoam -> property 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<double> stats(3*max_var_index, 0.0);
std::map<const char*, int> varids; // must agree with definition in _TaskPanelCfdResult.py
varids["Ux"] = 0;
varids["Uy"] = 1;
varids["Uz"] = 2;
varids["Umag"] = 3;
varids["Pressure"] = 4;
varids["Temperature"] = 5;
varids["TurbulenceEnergy"] = 6;
varids["TurbulenceViscosity"] = 7;
varids["TurbulenceDissipationRate"] = 8;
//varids["TurbulenceThermalDiffusivity"] = 9;
//varids["TurbulenceSpecificDissipation"] = 10;
double ts = 0.0; // t=0.0 for static simulation
static_cast<App::PropertyFloat*>(res->getPropertyByName("Time"))->setValue(ts);
vtkSmartPointer<vtkPointData> pd = dataset->GetPointData();
const vtkIdType nPoints = dataset->GetNumberOfPoints();
if(pd->GetNumberOfArrays() == 0) {
Base::Console().Error("No point data array is found in vtk data set, do nothin\n");
// if pointData is empty, data may be in cellDate, cellData -> pointData interpolation is possible in VTK
return;
}
std::vector<long> nodeIds(nPoints);
vtkSmartPointer<vtkDataArray> vel = pd->GetArray(vars["Velocity"]);
if(nPoints && vel && vel->GetNumberOfComponents() == 3) {
std::vector<Base::Vector3d> vec(nPoints);
double vmin=1.0e100, vmean=0.0, vmax=0.0;
// it is an internal utility func to avoid code duplication
void _calcStat(const std::vector<Base::Vector3d>& vel, std::vector<double>& stats) {
vtkIdType nPoints = vel.size();
double vmin=1.0e100, vmean=0.0, vmax=-1.0e100;
//stat of Vx, Vy, Vz is not necessary
double vmins[3] = {0.0, 0.0, 0.0};
double vmins[3] = {1.0e100, 1.0e100, 1.0e100}; // set up a very big positive float then reduce it
double vmeans[3] = {0.0, 0.0, 0.0};
double vmaxs[3] = {0.0, 0.0, 0.0};
for(vtkIdType i=0; i<nPoints; ++i) {
double *p = vel->GetTuple(i); // both vtkFloatArray and vtkDoubleArray return double* for GetTuple(i)
double vmaxs[3] = {-1.0e100, -1.0e100, -1.0e100};
for(std::vector<Base::Vector3d>::const_iterator it=vel.begin(); it!=vel.end(); ++it) {
double p[] = {it->x, it->y, it->z};
double vmag = std::sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
for(int ii=0; ii<3; ii++) {
vmeans[ii] += p[ii];
@@ -666,9 +644,6 @@ void FemVTKTools::importFluidicResult(vtkSmartPointer<vtkDataSet> dataset, App::
vmean += vmag;
if(vmag > vmax) vmax = vmag;
if(vmag < vmin) vmin = vmag;
vec[i] = (Base::Vector3d(p[0], p[1], p[2]));
nodeIds[i] = i;
}
for(int ii=0; ii<3; ii++) {
@@ -676,41 +651,99 @@ void FemVTKTools::importFluidicResult(vtkSmartPointer<vtkDataSet> dataset, App::
stats[ii*3 + 2] = vmaxs[ii];
stats[ii*3 + 1] = vmeans[ii]/nPoints;
}
int index = varids["Umag"];
int index = 3; // velocity mag or displacement mag
stats[index*3] = vmin;
stats[index*3 + 2] = vmax;
stats[index*3 + 1] = vmean/nPoints;
}
App::PropertyVectorList* velocity = static_cast<App::PropertyVectorList*>(res->getPropertyByName("Velocity"));
if(velocity) {
//PropertyVectorList will not show up in PropertyEditor
velocity->setValues(vec);
static_cast<App::PropertyIntegerList*>(res->getPropertyByName("NodeNumbers"))->setValues(nodeIds);
Base::Console().Message("Velocity field has been loaded \n");
}
else
Base::Console().Error("Velocity property is not found in Cfd Result object \n");
}
else {
Base::Console().Error("Velocity field is not found in Cfd Result vtk file \n");
void _importResult(const vtkSmartPointer<vtkDataSet> dataset, App::DocumentObject* res,
const std::map<std::string, std::string>& vectors, const std::map<std::string, std::string> scalers,
const std::map<std::string, int> varids, const std::string& essential_property){
const int max_var_index = 11;
// all code below can be shared!
std::vector<double> stats(3*max_var_index, 0.0);
double ts = 0.0; // t=0.0 for static simulation
static_cast<App::PropertyFloat*>(res->getPropertyByName("Time"))->setValue(ts);
vtkSmartPointer<vtkPointData> pd = dataset->GetPointData();
const vtkIdType nPoints = dataset->GetNumberOfPoints();
if(pd->GetNumberOfArrays() == 0) {
Base::Console().Error("No point data array is found in vtk data set, do nothing\n");
// if pointData is empty, data may be in cellDate, cellData -> pointData interpolation is possible in VTK
return;
}
for(auto const& kv: vars){
if (std::string(kv.first) == std::string("Velocity"))
continue;
vtkDataArray* vec = vtkDataArray::SafeDownCast(pd->GetArray(kv.second));
auto hGrp = App::GetApplication().GetParameterGroupByPath("User parameter:BaseApp/Preferences/Units");
int unitSchema = hGrp->GetInt("UserSchema",0);
float scale = 1.0;
if(unitSchema == 0) // standard mm
{
scale = 1000; // convert from SI length in result file to mm in FreeCAD
}
const char* essential_var = vectors.at(essential_property).c_str();
vtkSmartPointer<vtkDataArray> essential_array = pd->GetArray(essential_var); // a vector must exist
if(nPoints && essential_array) {
int dim = 3; // Fixme: currently 3D only
for(auto const& kv: vectors){
vtkDataArray* vector_field = vtkDataArray::SafeDownCast(pd->GetArray(kv.second.c_str()));
if(!vector_field)
vector_field = vtkDataArray::SafeDownCast(pd->GetArray(kv.first.c_str())); // name from FreeCAD export
if(vector_field && vector_field->GetNumberOfComponents() == dim) {
App::PropertyVectorList* vector_list = static_cast<App::PropertyVectorList*>(res->getPropertyByName(kv.first.c_str()));
if(vector_list) {
std::vector<Base::Vector3d> vec(nPoints);
if(kv.first == std::string(essential_property)) { // is there any other var need to scale?
for(vtkIdType i=0; i<nPoints; ++i) {
double *p = vector_field->GetTuple(i); // both vtkFloatArray and vtkDoubleArray return double* for GetTuple(i)
vec[i] = (Base::Vector3d(p[0]*scale, p[1]*scale, p[2]*scale));
}
}
else{
for(vtkIdType i=0; i<nPoints; ++i) {
double *p = vector_field->GetTuple(i); // both vtkFloatArray and vtkDoubleArray return double* for GetTuple(i)
vec[i] = (Base::Vector3d(p[0], p[1], p[2]));
}
}
if (kv.first == std::string(essential_property)) // for displacement or velocity calc min and max of each components
_calcStat(vec, stats);
//PropertyVectorList will not show up in PropertyEditor
vector_list->setValues(vec);
Base::Console().Message("PropertyVectorList %s has been loaded \n", kv.first.c_str());
}
else {
Base::Console().Error("static_cast<App::PropertyVectorList*>((res->getPropertyByName(\"%s\")) failed \n", kv.first.c_str());
continue;
}
}
std::vector<long> nodeIds(nPoints);
for(vtkIdType i=0; i<nPoints; ++i) {
nodeIds[i] = i+1;
}
static_cast<App::PropertyIntegerList*>(res->getPropertyByName("NodeNumbers"))->setValues(nodeIds);
}
}
else{
Base::Console().Error("essential_property %s corresponding essential array %s in VTK file is not found", essential_property.c_str(), essential_var);
}
for(auto const& kv: scalers){
vtkDataArray* vec = vtkDataArray::SafeDownCast(pd->GetArray(kv.second.c_str())); // name from OpenFOAM/Fem solver export
if(!vec)
vec = vtkDataArray::SafeDownCast(pd->GetArray(kv.first.c_str())); // name from FreeCAD export
if(nPoints && vec && vec->GetNumberOfComponents() == 1) {
App::PropertyFloatList* field = static_cast<App::PropertyFloatList*>(res->getPropertyByName(kv.first));
App::PropertyFloatList* field = static_cast<App::PropertyFloatList*>(res->getPropertyByName(kv.first.c_str()));
if (!field) {
Base::Console().Error("static_cast<App::PropertyFloatList*>((res->getPropertyByName(\"%s\")) failed \n", kv.first);
Base::Console().Error("static_cast<App::PropertyFloatList*>((res->getPropertyByName(\"%s\")) failed \n", kv.first.c_str());
continue;
}
double vmin=1.0e100, vmean=0.0, vmax=0.0;
double vmin=1.0e100, vmean=0.0, vmax=-1.0e100;
std::vector<double> values(nPoints, 0.0);
for(vtkIdType i = 0; i < vec->GetNumberOfTuples(); i++)
{
for(vtkIdType i = 0; i < vec->GetNumberOfTuples(); i++) {
double v = *(vec->GetTuple(i));
values[i] = v;
vmean += v;
@@ -719,199 +752,209 @@ void FemVTKTools::importFluidicResult(vtkSmartPointer<vtkDataSet> dataset, App::
}
field->setValues(values);
int index = varids[kv.first];
stats[index*3] = vmin;
stats[index*3 + 2] = vmax;
stats[index*3 + 1] = vmean/nPoints;
if(varids.find(kv.first) != varids.end()) {
const int index = varids.at(kv.first);
stats[index*3] = vmin;
stats[index*3 + 1] = vmean/nPoints;
stats[index*3 + 2] = vmax;
}
Base::Console().Message("field \"%s\" has been loaded \n", kv.first);
Base::Console().Message("field \"%s\" has been loaded \n", kv.first.c_str());
}
}
static_cast<App::PropertyFloatList*>(res->getPropertyByName("Stats"))->setValues(stats);
}
/*
void FemVTKTools::importMechanicalResult(const vtkDataSet* grid, App::DocumentObject* res) {
// to be implemented later by FemWorkbench developer
}
* */
void _exportResult(const App::DocumentObject* result, vtkSmartPointer<vtkDataSet> grid,
const std::map<std::string, std::string>& vectors, const std::map<std::string, std::string> scalers,
const std::string& essential_property){
void FemVTKTools::exportFluidicResult(const App::DocumentObject* res, vtkSmartPointer<vtkDataSet> grid) {
if(!res->getPropertyByName("Velocity")){
Base::Console().Message("Warning: essential field like `velocity` is not found in CfdResult\n");
return;
const Fem::FemResultObject* res = static_cast<const Fem::FemResultObject*>(result);
auto hGrp = App::GetApplication().GetParameterGroupByPath("User parameter:BaseApp/Preferences/Units");
int unitSchema = hGrp->GetInt("UserSchema",0);
float scale = 1.0;
if(unitSchema == 0) // standard mm
{
scale = 0.001; // convert from mm in FreeCAD to SI length in result file
}
App::PropertyVectorList* velocity = static_cast<App::PropertyVectorList*>(res->getPropertyByName("Velocity"));
const std::vector<Base::Vector3d>& vel = velocity->getValues();
if(!vel.empty()) {
vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
data->SetNumberOfComponents(3);
data->SetName("Velocity");
for(std::vector<Base::Vector3d>::const_iterator it=vel.begin(); it!=vel.end(); ++it) {
double tuple[] = {it->x, it->y, it->z};
data->InsertNextTuple(tuple);
const vtkIdType nPoints = grid->GetNumberOfPoints();
for (auto const& kv: vectors) {
const int dim = 3; //Fixme, detect dim
App::PropertyVectorList* field = nullptr;
if (res->getPropertyByName(kv.first.c_str()))
field = static_cast<App::PropertyVectorList*>(res->getPropertyByName(kv.first.c_str()));
else
Base::Console().Error("PropertyVectorList %s not found \n", kv.first.c_str());
if(field && field->getValues().size()>1) { // FreeCAD property list
const std::vector<Base::Vector3d>& vel = field->getValues();
vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
if(nPoints != field->getSize())
Base::Console().Error("PropertyVectorList->getSize() = %d, not equal to mesh point number \n", field->getSize());
data->SetNumberOfComponents(dim);
data->SetNumberOfTuples(vel.size());
data->SetName(kv.second.c_str()); // kv.first may be a better name, without space
vtkIdType i=0;
if(kv.first == essential_property) {
for(std::vector<Base::Vector3d>::const_iterator it=vel.begin(); it!=vel.end(); ++it) {
Base::Vector3d v = vel.at(i);
double tuple[] = {v.x*scale, v.y*scale, v.z*scale};
//double tuple[] = {it->x*scale, it->y*scale, it->z*scale};
data->SetTuple(i, tuple);
++i;
}
}
else{
for(std::vector<Base::Vector3d>::const_iterator it=vel.begin(); it!=vel.end(); ++it) {
double tuple[] = {it->x, it->y, it->z};
data->SetTuple(i, tuple);
++i;
}
}
grid->GetPointData()->AddArray(data);
Base::Console().Message("Info: PropertyVectorList %s exported as vtk array name '%s'\n", kv.first.c_str(), kv.second.c_str());
}
else
Base::Console().Error("field = static_cast<App::PropertyVectorList*> failed or empty for field: %s", kv.first.c_str());
}
grid->GetPointData()->AddArray(data);
}
else{
Base::Console().Message("Warning: essential fields pressure and velocity is empty in CfdResult\n");
}
// Temperature is optional, so are other turbulence related variables
std::vector<const char*> vars; // varable names are 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){
for (auto const& kv: scalers) {
App::PropertyFloatList* field = nullptr;
if (res->getPropertyByName(var))
field = static_cast<App::PropertyFloatList*>(res->getPropertyByName(var));
if(!field && !field->getValues().empty()) {
if (res->getPropertyByName(kv.first.c_str()))
field = static_cast<App::PropertyFloatList*>(res->getPropertyByName(kv.first.c_str()));
if(field && field->getValues().size()>1) {
const std::vector<double>& vec = field->getValues();
vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
data->SetNumberOfValues(vec.size());
data->SetName(var);
data->SetName(kv.second.c_str());
for(size_t i=0; i<vec.size(); ++i)
data->SetValue(i, vec[i]);
grid->GetPointData()->AddArray(data);
Base::Console().Message("Info: PropertyFloatList %s exported as vtk array name '%s'\n", kv.first.c_str(), kv.second.c_str());
}
}
}
void FemVTKTools::importFluidicResult(vtkSmartPointer<vtkDataSet> dataset, App::DocumentObject* res) {
// velocity and pressure are essential, Temperature is optional, so are turbulence related variables
std::map<std::string, std::string> cfd_vectors; // vector field defined in openfoam -> property defined in CfdResult.py
cfd_vectors["Velocity"] = "U";
std::map<std::string, std::string> cfd_scalers; // varable name defined in openfoam -> property defined in CfdResult.py
cfd_scalers["Pressure"] = "p";
cfd_scalers["Temperature"] = "T";
cfd_scalers["TurbulenceEnergy"] = "k";
cfd_scalers["TurbulenceViscosity"] = "nut";
cfd_scalers["TurbulenceDissipationRate"] = "epsilon";
cfd_scalers["TurbulenceSpecificDissipation"] = "omega";
cfd_scalers["TurbulenceThermalDiffusivity"] = "alphat";
std::map<std::string, int> cfd_varids; // must agree with definition in Stat calc Cfd/_TaskPanelCfdResult.py
cfd_varids["Ux"] = 0;
cfd_varids["Uy"] = 1;
cfd_varids["Uz"] = 2;
cfd_varids["Umag"] = 3;
cfd_varids["Pressure"] = 4;
cfd_varids["Temperature"] = 5;
cfd_varids["TurbulenceEnergy"] = 6;
cfd_varids["TurbulenceViscosity"] = 7;
cfd_varids["TurbulenceDissipationRate"] = 8;
//cfd_varids["TurbulenceSpecificDissipation"] = 9;
//cfd_varids["TurbulenceThermalDiffusivity"] = 10;
std::string essential_property = std::string("Velocity");
_importResult(dataset, res, cfd_vectors, cfd_scalers, cfd_varids, essential_property);
}
void FemVTKTools::exportFluidicResult(const App::DocumentObject* res, vtkSmartPointer<vtkDataSet> grid) {
// velocity and pressure are essential, Temperature is optional, so are turbulence related variables
static std::map<std::string, std::string> cfd_vectors; // vector field defined in openfoam -> property defined in CfdResult.py
cfd_vectors["Velocity"] = "U";
static std::map<std::string, std::string> cfd_scalers; // varable name defined in openfoam -> property defined in CfdResult.py
cfd_scalers["Pressure"] = "p";
cfd_scalers["Temperature"] = "T";
cfd_scalers["TurbulenceEnergy"] = "k";
cfd_scalers["TurbulenceViscosity"] = "nut";
cfd_scalers["TurbulenceDissipationRate"] = "epsilon";
cfd_scalers["TurbulenceSpecificDissipation"] = "omega";
cfd_scalers["TurbulenceThermalDiffusivity"] = "alphat";
std::string essential_property = std::string("Velocity");
if(!res->getPropertyByName("Velocity")){
Base::Console().Error("essential field like `velocity` is not found in CfdResult\n");
return;
}
_exportResult(res, grid, cfd_vectors, cfd_scalers, essential_property);
}
void FemVTKTools::exportMechanicalResult(const App::DocumentObject* obj, vtkSmartPointer<vtkDataSet> grid) {
// code redundance can be avoided by property inspection, consider refactoring
const FemResultObject* res = static_cast<const FemResultObject*>(obj);
if(!res->StressValues.getValues().empty()) {
const std::vector<double>& vec = res->StressValues.getValues();
if (vec.size()>1) {
vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
data->SetNumberOfValues(vec.size());
data->SetName("Von Mises stress");
void FemVTKTools::importMechanicalResult(vtkSmartPointer<vtkDataSet> dataset, App::DocumentObject* res) {
// field names are defined in this cpp, exportMechanicalResult()
// DisplaceVectors are essential, Temperature and other is optional
std::map<std::string, std::string> vectors; // property defined in MechanicalResult.py -> variable name in vtk
vectors["DisplacementVectors"] = "Displacement";
vectors["StrainVectors"] = "Strain vectors";
vectors["StressVectors"] = "Stress vectors";
std::map<std::string, std::string> scalers; // App::FloatListProperty name -> vtk name
scalers["UserDefined"] = "User Defined Results";
scalers["Temperature"] = "Temperature";
scalers["PrincipalMax"] = "Maximum Principal stress";
scalers["PrincipalMed"] = "Median Principal stress";
scalers["PrincipalMin"] = "Minimum Principal stress";
scalers["MaxShear"] = "Max shear stress (Tresca)";
scalers["StressValues"] = "Von Mises stress";
//scalers["DisplacementLengths"] = ""; // not yet exported in exportMechanicalResult()
for(size_t i=0; i<vec.size(); ++i)
data->SetValue(i, vec[i]);
std::map<std::string, int> varids;
// id sequence must agree with definition in get_result_stats() of Fem/_TaskPanelShowResult.py
varids["U1"] = 0; // U1, displacement x axis
varids["U2"] = 1;
varids["U3"] = 2;
varids["Uabs"] = 3;
varids["StressValues"] = 4; // Sabs
varids["PrincipalMax"] = 5; // MaxPrin
varids["PrincipalMed"] = 6; // MidPrin
varids["PrincipalMin"] = 7; // MinPrin
varids["MaxShear"] = 8; //
grid->GetPointData()->AddArray(data);
}}
std::string essential_property = std::string("DisplacementVectors");
if(!res->MaxShear.getValues().empty()) {
const std::vector<double>& vec = res->MaxShear.getValues();
if (vec.size()>1) {
vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
data->SetNumberOfValues(vec.size());
data->SetName("Max shear stress (Tresca)");
_importResult(dataset, res, vectors, scalers, varids, essential_property);
for(size_t i=0; i<vec.size(); ++i)
data->SetValue(i, vec[i]);
grid->GetPointData()->AddArray(data);
}}
if(!res->PrincipalMax.getValues().empty()) {
const std::vector<double>& vec = res->PrincipalMax.getValues();
if (vec.size()>1) {
vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
data->SetNumberOfValues(vec.size());
data->SetName("Maximum Principal stress");
for(size_t i=0; i<vec.size(); ++i)
data->SetValue(i, vec[i]);
grid->GetPointData()->AddArray(data);
}}
if(!res->PrincipalMax.getValues().empty()) {
const std::vector<double>& vec = res->PrincipalMin.getValues();
if (vec.size()>1) {
vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
data->SetNumberOfValues(vec.size());
data->SetName("Minimum Principal stress");
for(size_t i=0; i<vec.size(); ++i)
data->SetValue(i, vec[i]);
grid->GetPointData()->AddArray(data);
}}
if (!res->Temperature.getValues().empty()) {
const std::vector<double>& vec = res->Temperature.getValues();
if (vec.size()>1) {
vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
data->SetNumberOfValues(vec.size());
data->SetName("Temperature");
for(size_t i=0; i<vec.size(); ++i)
data->SetValue(i, vec[i]);
grid->GetPointData()->AddArray(data);
}}
if (!res->UserDefined.getValues().empty()) {
const std::vector<double>& vec = res->UserDefined.getValues();
if (vec.size()>1) {
vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
data->SetNumberOfValues(vec.size());
data->SetName("User Defined Results");
for(size_t i=0; i<vec.size(); ++i)
data->SetValue(i, vec[i]);
grid->GetPointData()->AddArray(data);
}}
}
if(!res->DisplacementVectors.getValues().empty()) {
const std::vector<Base::Vector3d>& vec = res->DisplacementVectors.getValues();
if (vec.size()>1) {
vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
data->SetNumberOfComponents(3);
data->SetName("Displacement");
void FemVTKTools::exportMechanicalResult(const App::DocumentObject* res, vtkSmartPointer<vtkDataSet> grid) {
if(!res->getPropertyByName("DisplacementVectors")){
Base::Console().Error("essential field like `DisplacementVectors` is not found in this Result object\n");
return;
}
std::map<std::string, std::string> vectors; // property defined in MechanicalResult.py -> variable name in vtk
vectors["DisplacementVectors"] = "Displacement";
vectors["StrainVectors"] = "Strain vectors";
vectors["StressVectors"] = "Stress vectors";
std::map<std::string, std::string> scalers; // App::FloatListProperty name -> vtk name
scalers["UserDefined"] = "User Defined Results";
scalers["Temperature"] = "Temperature";
scalers["PrincipalMax"] = "Maximum Principal stress";
scalers["PrincipalMed"] = "Median Principal stress";
scalers["PrincipalMin"] = "Minimum Principal stress";
scalers["MaxShear"] = "Max shear stress (Tresca)";
scalers["StressValues"] = "Von Mises stress";
//scalers["DisplacementLengths"] = ""; // not yet exported in exportMechanicalResult()
for(std::vector<Base::Vector3d>::const_iterator it=vec.begin(); it!=vec.end(); ++it) {
double tuple[] = {it->x, it->y, it->z};
data->InsertNextTuple(tuple);
}
grid->GetPointData()->AddArray(data);
}}
if(!res->StressVectors.getValues().empty()) {
const std::vector<Base::Vector3d>& vec = res->StressVectors.getValues();
if (vec.size()>1) {
vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
data->SetNumberOfComponents(3);
data->SetName("Stress Vectors");
for(std::vector<Base::Vector3d>::const_iterator it=vec.begin(); it!=vec.end(); ++it) {
double tuple[] = {it->x, it->y , it->z};
data->InsertNextTuple(tuple);
}
grid->GetPointData()->AddArray(data);
}}
if(!res->StrainVectors.getValues().empty()) {
const std::vector<Base::Vector3d>& vec = res->StrainVectors.getValues();
if (vec.size()>1) {
vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
data->SetNumberOfComponents(3);
data->SetName("Strain Vectors");
for(std::vector<Base::Vector3d>::const_iterator it=vec.begin(); it!=vec.end(); ++it) {
double tuple[] = {it->x, it->y, it->z};
data->InsertNextTuple(tuple);
}
grid->GetPointData()->AddArray(data);
}}
std::string essential_property = std::string("DisplacementVectors");
_exportResult(res, grid, vectors, scalers, essential_property);
}

8
src/Mod/Fem/App/FemVTKTools.h
View File

@@ -54,7 +54,7 @@ namespace Fem
/*!
FemMesh export to vtkUnstructuredGrid instance
*/
static void exportVTKMesh(const FemMesh* mesh, vtkSmartPointer<vtkUnstructuredGrid> grid);
static void exportVTKMesh(const FemMesh* mesh, vtkSmartPointer<vtkUnstructuredGrid> grid, float scale = 1.0);
/*!
FemMesh write to vtkUnstructuredGrid data file
*/
@@ -64,7 +64,7 @@ namespace Fem
* FemResult import from vtkUnstructuredGrid object
*/
static void importFluidicResult(vtkSmartPointer<vtkDataSet> dataset, App::DocumentObject* res);
// importMechanicalResult can be defined if necessary in the future
static void importMechanicalResult(vtkSmartPointer<vtkDataSet> dataset, App::DocumentObject* res);
/*!
* FemResult export to vtkUnstructuredGrid object
@@ -73,9 +73,9 @@ namespace Fem
static void exportMechanicalResult(const App::DocumentObject* res, vtkSmartPointer<vtkDataSet> grid);
/*!
* FemResult (active or created if res= NULL) read from vtkUnstructuredGrid dataset file
* FemResult (activeObject or created if res= NULL) read from vtkUnstructuredGrid dataset file
*/
static App::DocumentObject* readFluidicResult(const char* Filename, App::DocumentObject* res = NULL);
static App::DocumentObject* readResult(const char* Filename, App::DocumentObject* res = NULL);
/*!
* write FemResult (activeObject if res= NULL) to vtkUnstructuredGrid dataset file

22
src/Mod/Fem/Gui/Command.cpp
View File

@@ -1525,13 +1525,25 @@ CmdFemPostPipelineFromResult::CmdFemPostPipelineFromResult()
void CmdFemPostPipelineFromResult::activated(int)
{
/*
Gui::SelectionFilter ResultFilter("SELECT Fem::FemResultObject COUNT 1");
if (ResultFilter.match()) {
Base::Console().Message("Debug: `SELECT Fem::FemResultObject COUNT 1` has matched obj");
Fem::FemResultObject* result = static_cast<Fem::FemResultObject*>(ResultFilter.Result[0][0].getObject());
//static_cast failed here
Base::Console().Message("Debug: FemResultObject pointer = %p", result );
*/
App::Document* pcDoc = App::GetApplication().getActiveDocument();
if(!pcDoc)
{
Base::Console().Message("No active document is found thus do nothing and return\n");
return;
}
Fem::FemResultObject* result= static_cast<Fem::FemResultObject*>(pcDoc->getActiveObject());
if(result)
{
std::string FeatName = getUniqueObjectName("Pipeline");
openCommand("Create pipeline from result");
doCommand(Doc,"App.activeDocument().addObject('Fem::FemPostPipeline','%s')",FeatName.c_str());
@@ -1543,8 +1555,8 @@ void CmdFemPostPipelineFromResult::activated(int)
}
else {
QMessageBox::warning(Gui::getMainWindow(),
qApp->translate("CmdFemPostCreateClipFilter", "Wrong selection"),
qApp->translate("CmdFemPostCreateClipFilter", "Select a result, please."));
qApp->translate("CmdFemPostPipelineFromResult", "Wrong selection type"),
qApp->translate("CmdFemPostPipelineFromResult", "Select a result object, please."));
}
}