From: Anthony Geay Date: Mon, 17 Sep 2018 13:34:15 +0000 (+0200) Subject: [EDF17832] : Move VoroGauss GaussToCell StaticMesh DevelopedSurface from EDF to parav... X-Git-Tag: V9_2_0a1~1 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=7356ff2dc331f1109527b7776eb9eb40eac880f1;p=modules%2Fparavis.git [EDF17832] : Move VoroGauss GaussToCell StaticMesh DevelopedSurface from EDF to paravis. Cotech104 --- diff --git a/src/Plugins/CMakeLists.txt b/src/Plugins/CMakeLists.txt index c60e0645..5cf700c4 100755 --- a/src/Plugins/CMakeLists.txt +++ b/src/Plugins/CMakeLists.txt @@ -30,6 +30,10 @@ SET(_subdirs DifferenceTimesteps ArrayRenamer JSONReader + DevelopedSurface + StaticMesh + GaussToCell + VoroGauss ) IF(NOT SALOME_LIGHT_ONLY) diff --git a/src/Plugins/DevelopedSurface/CMakeLists.txt b/src/Plugins/DevelopedSurface/CMakeLists.txt new file mode 100644 index 00000000..5db84d5e --- /dev/null +++ b/src/Plugins/DevelopedSurface/CMakeLists.txt @@ -0,0 +1,70 @@ +# Copyright (C) 2017 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +PROJECT(DevelopedSurface) +CMAKE_MINIMUM_REQUIRED(VERSION 2.8) + +CMAKE_POLICY(SET CMP0003 NEW) +IF(${CMAKE_VERSION} VERSION_GREATER "3.0.0") + CMAKE_POLICY(SET CMP0022 OLD) + CMAKE_POLICY(SET CMP0023 OLD) +ENDIF() + + +ENABLE_TESTING() + +# Common CMake macros +# =================== +SET(CONFIGURATION_ROOT_DIR $ENV{CONFIGURATION_ROOT_DIR} CACHE PATH "Path to the Salome CMake configuration files") +IF(EXISTS ${CONFIGURATION_ROOT_DIR}) + LIST(APPEND CMAKE_MODULE_PATH "${CONFIGURATION_ROOT_DIR}/cmake") + INCLUDE(SalomeMacros) +ELSE() + MESSAGE(FATAL_ERROR "We absolutely need the Salome CMake configuration files, please define CONFIGURATION_ROOT_DIR !") +ENDIF() +FIND_PACKAGE(SalomePythonInterp REQUIRED) +FIND_PACKAGE(SalomePythonLibs REQUIRED) + +FIND_PACKAGE(ParaView REQUIRED) +IF(NOT ParaView_FOUND) + MESSAGE(FATAL_ERROR "Please locate ParaView." ) +ENDIF(NOT ParaView_FOUND) +INCLUDE(${PARAVIEW_USE_FILE}) + +SET(MEDCOUPLING_ROOT_DIR $ENV{MEDCOUPLING_ROOT_DIR} CACHE PATH "MEDCOUPLING_ROOT_DIR") +LIST(APPEND CMAKE_MODULE_PATH "${MEDCOUPLING_ROOT_DIR}/cmake_files") + +## +FIND_PACKAGE(SalomeMEDCoupling REQUIRED) + +OPTION(BUILD_SHARED_LIBS "Build with shared libraries." ${VTK_BUILD_SHARED_LIBS}) + +SET(VTK_INSTALL_RUNTIME_DIR lib) +SET(VTK_INSTALL_LIBRARY_DIR lib) +SET(VTK_INSTALL_ARCHIVE_DIR lib) + +PV_PROCESS_MODULES() +INCLUDE_DIRECTORIES( ${MEDCOUPLING_INCLUDE_DIRS} ) +INCLUDE_DIRECTORIES( "${PROJECT_SOURCE_DIR}/../MEDWriter/IO" ) +#INCLUDE_DIRECTORIES( "${PARAVIS_ROOT_DIR}/include/salome" ) +ADD_SUBDIRECTORY(ParaViewPlugin) + + +ADD_SUBDIRECTORY(Test) diff --git a/src/Plugins/DevelopedSurface/IO/vtkDevelopedSurface.cxx b/src/Plugins/DevelopedSurface/IO/vtkDevelopedSurface.cxx new file mode 100644 index 00000000..78872ed0 --- /dev/null +++ b/src/Plugins/DevelopedSurface/IO/vtkDevelopedSurface.cxx @@ -0,0 +1,520 @@ +// Copyright (C) 2017 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 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 +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay (EDF R&D) + +#include "vtkDevelopedSurface.h" +#include "VTKToMEDMem.hxx" + +#include "vtkAdjacentVertexIterator.h" +#include "vtkIntArray.h" +#include "vtkCellData.h" +#include "vtkPointData.h" +#include "vtkCylinder.h" +#include "vtkNew.h" +#include "vtkCutter.h" +#include "vtkTransform.h" + +#include "vtkStreamingDemandDrivenPipeline.h" +#include "vtkUnstructuredGrid.h" +#include "vtkMultiBlockDataSet.h" + +#include "vtkInformationStringKey.h" +#include "vtkAlgorithmOutput.h" +#include "vtkObjectFactory.h" +#include "vtkMutableDirectedGraph.h" +#include "vtkMultiBlockDataSet.h" +#include "vtkDataSet.h" +#include "vtkInformationVector.h" +#include "vtkInformation.h" +#include "vtkDataArraySelection.h" +#include "vtkTimeStamp.h" +#include "vtkInEdgeIterator.h" +#include "vtkInformationDataObjectKey.h" +#include "vtkExecutive.h" +#include "vtkVariantArray.h" +#include "vtkStringArray.h" +#include "vtkDoubleArray.h" +#include "vtkFloatArray.h" +#include "vtkCharArray.h" +#include "vtkUnsignedCharArray.h" +#include "vtkDataSetAttributes.h" +#include "vtkDemandDrivenPipeline.h" +#include "vtkDataObjectTreeIterator.h" +#include "vtkWarpScalar.h" + +#include "MEDCouplingMemArray.hxx" + +#include "VTKMEDTraits.hxx" + +#include +#include +#include +#include + +vtkStandardNewMacro(vtkDevelopedSurface); + +/////////////////// + +template +struct VTKTraits +{ +}; + +template<> +struct VTKTraits +{ + typedef vtkDoubleArray ArrayType; +}; + +template<> +struct VTKTraits +{ + typedef vtkFloatArray ArrayType; +}; + +void ExtractInfo(vtkInformationVector *inputVector, vtkDataSet *& usgIn) +{ + vtkInformation *inputInfo(inputVector->GetInformationObject(0)); + vtkDataSet *input(0); + vtkDataSet *input0(vtkDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT()))); + vtkMultiBlockDataSet *input1(vtkMultiBlockDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT()))); + if(input0) + input=input0; + else + { + if(!input1) + throw MZCException("Input dataSet must be a DataSet or single elt multi block dataset expected !"); + if(input1->GetNumberOfBlocks()!=1) + throw MZCException("Input dataSet is a multiblock dataset with not exactly one block ! Use MergeBlocks or ExtractBlocks filter before calling this filter !"); + vtkDataObject *input2(input1->GetBlock(0)); + if(!input2) + throw MZCException("Input dataSet is a multiblock dataset with exactly one block but this single element is NULL !"); + vtkDataSet *input2c(vtkDataSet::SafeDownCast(input2)); + if(!input2c) + throw MZCException("Input dataSet is a multiblock dataset with exactly one block but this single element is not a dataset ! Use MergeBlocks or ExtractBlocks filter before calling this filter !"); + input=input2c; + } + if(!input) + throw MZCException("Input data set is NULL !"); + vtkPointData *att(input->GetPointData()); + if(!att) + throw MZCException("Input dataset has no point data attribute ! Impossible to deduce a developed surface on it !"); + usgIn=input; +} + +class vtkDevelopedSurface::vtkInternals +{ +public: + vtkNew Cutter; +}; + +//////////////////// + +vtkDevelopedSurface::vtkDevelopedSurface():_cyl(nullptr),Internal(new vtkInternals),InvertStatus(false),OffsetInRad(0.) +{ + //this->RegisterFilter(this->Internal->Cutter.GetPointer()); +} + +vtkDevelopedSurface::~vtkDevelopedSurface() +{ + delete this->Internal; +} + +void vtkDevelopedSurface::SetInvertWay(bool invertStatus) +{ + this->InvertStatus=invertStatus; + this->Modified(); +} + +void vtkDevelopedSurface::SetThetaOffset(double offsetInDegrees) +{ + double tmp(std::min(offsetInDegrees,180.)); + tmp=std::max(tmp,-180.); + this->OffsetInRad=tmp/180.*M_PI; + this->Modified(); +} + +int vtkDevelopedSurface::RequestInformation(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) +{ + //std::cerr << "########################################## vtkDevelopedSurface::RequestInformation ##########################################" << std::endl; + try + { + vtkDataSet *usgIn(0); + ExtractInfo(inputVector[0],usgIn); + } + catch(MZCException& e) + { + std::ostringstream oss; + oss << "Exception has been thrown in vtkDevelopedSurface::RequestInformation : " << e.what() << std::endl; + if(this->HasObserver("ErrorEvent") ) + this->InvokeEvent("ErrorEvent",const_cast(oss.str().c_str())); + else + vtkOutputWindowDisplayErrorText(const_cast(oss.str().c_str())); + vtkObject::BreakOnError(); + return 0; + } + return 1; +} + +std::vector UnWrapByDuplicatingNodes(vtkCellArray *ca, vtkIdType& offset, const MEDCoupling::DataArrayDouble *thetas) +{ + std::vector ret; + vtkIdType nbCells(ca->GetNumberOfCells()); + vtkIdType *conn(ca->GetPointer()); + const double *tptr(thetas->begin()); + for(vtkIdType i=0;i::max()),mi(std::numeric_limits::max()); + for(vtkIdType j=0;jM_PI) + { + for(vtkIdType j=0;j +void DealArray(vtkDataSetAttributes *pd, int pos, typename MEDFileVTKTraits::VtkType *arr, std::vector& nodeSel) +{ + int nbc(arr->GetNumberOfComponents()); + std::size_t nbt(nodeSel.size()); + vtkSmartPointer< typename MEDFileVTKTraits::VtkType > newArr; + newArr.TakeReference(MEDFileVTKTraits::VtkType::New()); + newArr->SetNumberOfComponents(nbc); + newArr->SetNumberOfTuples(nbt); + T *ptr(newArr->GetPointer(0)); + const T *inPtr(arr->GetPointer(0)); + for(std::size_t i=0;iSetName(arr->GetName()); + arr->DeepCopy(newArr); +} + +void ToDouble(vtkDataArray *coords, vtkSmartPointer& coordsOut) +{ + vtkDoubleArray *coords2(vtkDoubleArray::SafeDownCast(coords)); + vtkFloatArray *coords3(vtkFloatArray::SafeDownCast(coords)); + if(!coords2 && !coords3) + throw MZCException("Input coordinates are neither float64 or float32 !"); + // + if(coords2) + { + coordsOut.TakeReference(coords2); + coords2->Register(0); + } + else + { + coordsOut.TakeReference(vtkDoubleArray::New()); + coordsOut->SetNumberOfComponents(3); + vtkIdType nbTuples(coords3->GetNumberOfTuples()); + coordsOut->SetNumberOfTuples(nbTuples); + std::copy(coords3->GetPointer(0),coords3->GetPointer(0)+3*nbTuples,coordsOut->GetPointer(0)); + } +} + +void dealWith(vtkPolyData *outdata, const double center[3], const double axis[3], double radius, double eps, bool invertThetaInc, double offsetInRad) +{ + vtkDataArray *coords(outdata->GetPoints()->GetData()); + if(coords->GetNumberOfComponents()!=3) + throw MZCException("Input coordinates are expected to have 3 components !"); + // + vtkIdType nbNodes(coords->GetNumberOfTuples()); + if(nbNodes==0) + throw MZCException("No points -> impossible to develop anything !"); + // + vtkSmartPointer zeCoords; + ToDouble(coords,zeCoords); + // + double axis_cross_Z[3]={axis[1],-axis[0],0.}; + double n_axis(sqrt(axis_cross_Z[0]*axis_cross_Z[0]+axis_cross_Z[1]*axis_cross_Z[1])); + if(n_axis>eps) + { + double ang(asin(n_axis)); + if(axis[2]<0.) + ang=M_PI-ang; + MEDCoupling::DataArrayDouble::Rotate3DAlg(center,axis_cross_Z,ang,nbNodes,zeCoords->GetPointer(0),zeCoords->GetPointer(0)); + } + // + MEDCoupling::MCAuto c_cyl; + { + MEDCoupling::MCAuto cc(MEDCoupling::DataArrayDouble::New()); cc->alloc(nbNodes,3); + double *ccPtr(cc->getPointer()); + const double *zeCoordsPtr(zeCoords->GetPointer(0)); + for(vtkIdType i=0;i()); + } + c_cyl=cc->fromCartToCyl(); + } + MEDCoupling::MCAuto mfd(MEDCoupling::MEDFileData::New()); + WriteMEDFileFromVTKDataSet(mfd,outdata,{},0.,0); + bool a; + { + MEDCoupling::MEDFileMeshes *ms(mfd->getMeshes()); + if(ms->getNumberOfMeshes()!=1) + throw MZCException("Unexpected number of meshes !"); + MEDCoupling::MEDFileMesh *mm(ms->getMeshAtPos(0)); + MEDCoupling::MEDFileUMesh *mmu(dynamic_cast(mm)); + if(!mmu) + throw MZCException("Expecting unstructured one !"); + MEDCoupling::MCAuto m0(mmu->getMeshAtLevel(0)); + { + int v(0); + MEDCoupling::MCAuto c0s(m0->getCellIdsLyingOnNodes(&v,&v+1,false)); + if(c0s->empty()) + throw MZCException("Orphan node 0 !"); + std::vector nodes0; + m0->getNodeIdsOfCell(c0s->getIJ(0,0),nodes0); + MEDCoupling::MCAuto tmp0(c_cyl->selectByTupleIdSafe(nodes0.data(),nodes0.data()+nodes0.size())); + tmp0=tmp0->keepSelectedComponents({1}); + double tmp(tmp0->getMaxAbsValueInArray()); + a=tmp>0.; + } + } + // + constexpr double EPS_FOR_RADIUS=1e-2; + MEDCoupling::MCAuto rs(c_cyl->keepSelectedComponents({0})); + if(!rs->isUniform(radius,radius*EPS_FOR_RADIUS)) + { + double mi(rs->getMinValueInArray()),ma(rs->getMaxValueInArray()); + std::ostringstream oss; oss << "Looks not really a cylinder within given precision ! Range is [" << mi << "," << ma << "] expecting " << radius << " within precision of " << radius*EPS_FOR_RADIUS << " !"; + throw MZCException(oss.str()); + } + double tetha0(c_cyl->getIJ(0,1)); + { + double *ccylptr(c_cyl->getPointer()+1); + double mi02(std::numeric_limits::max()); + for(vtkIdType i=0;i0.) + { + ccylptr=c_cyl->getPointer()+1; + for(vtkIdType i=0;i=2*M_PI) + { + *ccylptr+=-2*M_PI; + } + } + } + } + { + MEDCoupling::MCAuto c_cyl_2(c_cyl->keepSelectedComponents({1})); + c_cyl_2->abs(); + MEDCoupling::MCAuto poses(c_cyl_2->findIdsInRange(0.,eps)); + c_cyl->setPartOfValuesSimple3(0.,poses->begin(),poses->end(),1,2,1); + } + // + if(a ^ (!invertThetaInc)) + { + MEDCoupling::MCAuto tmp(c_cyl->keepSelectedComponents({1})); + tmp=tmp->negate(); + std::for_each(tmp->getPointer(),tmp->getPointer()+tmp->getNumberOfTuples(),[](double& v) { if(v==-0.) v=0.; }); + c_cyl->setPartOfValues1(tmp,0,nbNodes,1,1,2,1); + } + MEDCoupling::MCAuto c_cyl_post(c_cyl->keepSelectedComponents({1})); + { + double *c_cyl_post_ptr(c_cyl_post->getPointer()); + for(vtkIdType i=0;iGetPolys()); + vtkIdType offset(nbNodes); + std::vector dupNodes(UnWrapByDuplicatingNodes(cb,offset,c_cyl_post)); + // + MEDCoupling::MCAuto c_cyl_post2(c_cyl_post->selectByTupleId(dupNodes.data(),dupNodes.data()+dupNodes.size())); + c_cyl_post2->applyLin(1.,2*M_PI); + c_cyl_post=MEDCoupling::DataArrayDouble::Aggregate(c_cyl_post,c_cyl_post2); + MEDCoupling::MCAuto z0(c_cyl->keepSelectedComponents({2})); + MEDCoupling::MCAuto z1(z0->selectByTupleId(dupNodes.data(),dupNodes.data()+dupNodes.size())); + z0=MEDCoupling::DataArrayDouble::Aggregate(z0,z1); + // + std::size_t outNbNodes(z0->getNumberOfTuples()); + vtkSmartPointer zeCoords2; + zeCoords2.TakeReference(vtkDoubleArray::New()); + zeCoords2->SetNumberOfComponents(3); + zeCoords2->SetNumberOfTuples(outNbNodes); + { + const double *tptr(c_cyl_post->begin()),*zptr(z0->begin()); + double *outPtr(zeCoords2->GetPointer(0)); + for(std::size_t i=0;iGetPoints()->SetData(zeCoords2); + // now post process nodes + std::vector nodeSel(nbNodes+dupNodes.size()); + { + int cnt(0); + std::for_each(nodeSel.begin(),nodeSel.begin()+nbNodes,[&cnt](int& v){ v=cnt++; }); + std::copy(dupNodes.begin(),dupNodes.end(),nodeSel.begin()+nbNodes); + } + vtkDataSetAttributes *pd(outdata->GetPointData()); + int nba(pd->GetNumberOfArrays()); + for(int i=0;iGetArray(i)); + { + vtkIntArray *arr0(vtkIntArray::SafeDownCast(arr)); + if(arr0) + { + DealArray(pd,i,arr0,nodeSel); + continue; + } + } + { + vtkFloatArray *arr0(vtkFloatArray::SafeDownCast(arr)); + if(arr0) + { + DealArray(pd,i,arr0,nodeSel); + continue; + } + } + { + vtkDoubleArray *arr0(vtkDoubleArray::SafeDownCast(arr)); + if(arr0) + { + DealArray(pd,i,arr0,nodeSel); + continue; + } + } + } +} + +int vtkDevelopedSurface::RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) +{ + //std::cerr << "########################################## vtkDevelopedSurface::RequestData ##########################################" << std::endl; + try + { + if(!_cyl) + throw MZCException("No cylinder object as cut function !"); + double center[3],axis[3],radius; + vtkAbstractTransform* trf(_cyl->GetTransform()); + { + _cyl->GetCenter(center); + _cyl->GetAxis(axis[0],axis[1],axis[2]); + radius=_cyl->GetRadius(); + } + if(trf) + { + double axis3[3]={center[0]+0.,center[1]+1.,center[2]+0.},axis4[3]; + trf->TransformPoint(axis3,axis4); + std::transform(axis4,axis4+3,center,axis,[](double a, double b) { return b-a; }); + axis[1]=-axis[1]; + if(std::isnan(axis[0]) && std::isnan(axis[1]) && std::isnan(axis[2])) + { axis[0]=0.; axis[1]=-1.; axis[2]=0.; } + } + //std::cerr << trf << " jjj " << axis[0] << " " << axis[1] << " " << axis[2] << " : " << center[0] << " " << center[1] << " " << center[2] << " " " " << " -> " << radius << std::endl; + vtkDataSet *usgIn(0); + ExtractInfo(inputVector[0],usgIn); + vtkSmartPointer outData; + { + vtkNew Cutter; + Cutter->SetInputData(usgIn); + Cutter->SetCutFunction(_cyl); + Cutter->Update(); + vtkDataSet *zeComputedOutput(Cutter->GetOutput()); + vtkPolyData *zeComputedOutput2(vtkPolyData::SafeDownCast(zeComputedOutput)); + if(!zeComputedOutput2) + throw MZCException("Unexpected output of cutter !"); + outData.TakeReference(zeComputedOutput2); + zeComputedOutput2->Register(0); + } + if(outData->GetNumberOfCells()==0) + return 1;// no cells -> nothing to do + // + dealWith(outData,center,axis,radius,1e-7,this->InvertStatus,this->OffsetInRad); + //finish + vtkInformation *outInfo(outputVector->GetInformationObject(0)); + vtkPolyData *output(vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()))); + output->ShallowCopy(outData); + } + catch(MZCException& e) + { + std::ostringstream oss; + oss << "Exception has been thrown in vtkDevelopedSurface::RequestInformation : " << e.what() << std::endl; + if(this->HasObserver("ErrorEvent") ) + this->InvokeEvent("ErrorEvent",const_cast(oss.str().c_str())); + else + vtkOutputWindowDisplayErrorText(const_cast(oss.str().c_str())); + vtkObject::BreakOnError(); + return 0; + } + return 1; +} + +int vtkDevelopedSurface::FillOutputPortInformation( int vtkNotUsed(port), vtkInformation* info) +{ + info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData"); + return 1; +} + + +void vtkDevelopedSurface::PrintSelf(ostream& os, vtkIndent indent) +{ + this->Superclass::PrintSelf(os, indent); +} + +void vtkDevelopedSurface::SetCutFunction(vtkImplicitFunction* func) +{ + vtkCylinder *cyl(vtkCylinder::SafeDownCast(func)); + if(cyl) + { + _cyl=cyl; + this->Modified(); + } +} + +vtkMTimeType vtkDevelopedSurface::GetMTime() +{ + vtkMTimeType maxMTime = this->Superclass::GetMTime(); // My MTime + if(_cyl) + { + maxMTime=std::max(maxMTime,_cyl->GetMTime()); + } + return maxMTime; +} diff --git a/src/Plugins/DevelopedSurface/IO/vtkDevelopedSurface.h b/src/Plugins/DevelopedSurface/IO/vtkDevelopedSurface.h new file mode 100644 index 00000000..dfd4126a --- /dev/null +++ b/src/Plugins/DevelopedSurface/IO/vtkDevelopedSurface.h @@ -0,0 +1,63 @@ +// Copyright (C) 2017 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 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 +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay (EDF R&D) + +#ifndef vtkDevelopedSurface_h__ +#define vtkDevelopedSurface_h__ + +#include "vtkDataSetAlgorithm.h" + +class vtkMutableDirectedGraph; +class vtkImplicitFunction; +class vtkCylinder; + +class VTK_EXPORT vtkDevelopedSurface : public vtkDataSetAlgorithm +{ +public: + static vtkDevelopedSurface* New(); + vtkTypeMacro(vtkDevelopedSurface, vtkDataSetAlgorithm) + void PrintSelf(ostream& os, vtkIndent indent); + void SetCutFunction(vtkImplicitFunction* func); + vtkMTimeType GetMTime(); + void SetInvertWay(bool invertStatus); + void SetThetaOffset(double offsetInDegrees); +protected: + vtkDevelopedSurface(); + ~vtkDevelopedSurface(); + int FillOutputPortInformation( int vtkNotUsed(port), vtkInformation* info); + int RequestInformation(vtkInformation *request, + vtkInformationVector **inputVector, vtkInformationVector *outputVector); + + int RequestData(vtkInformation *request, vtkInformationVector **inputVector, + vtkInformationVector *outputVector); +private: + vtkDevelopedSurface(const vtkDevelopedSurface&); + void operator=(const vtkDevelopedSurface&); // Not implemented. + private: + //BTX + vtkCylinder *_cyl; + //ETX + + class vtkInternals; + vtkInternals *Internal; + bool InvertStatus; + double OffsetInRad; +}; + +#endif diff --git a/src/Plugins/DevelopedSurface/ParaViewPlugin/CMakeLists.txt b/src/Plugins/DevelopedSurface/ParaViewPlugin/CMakeLists.txt new file mode 100644 index 00000000..f84b17b1 --- /dev/null +++ b/src/Plugins/DevelopedSurface/ParaViewPlugin/CMakeLists.txt @@ -0,0 +1,26 @@ +# Copyright (C) 2017 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +INCLUDE_DIRECTORIES( ${CMAKE_CURRENT_SOURCE_DIR}/../IO ) +ADD_PARAVIEW_PLUGIN(DevelopedSurfacePlugin "4.0" + SERVER_MANAGER_SOURCES ${SM_SRCS} ${PROJECT_SOURCE_DIR}/IO/vtkDevelopedSurface.cxx + SERVER_MANAGER_XML Resources/DevelopedSurfaceServer.xml) +TARGET_LINK_LIBRARIES(DevelopedSurfacePlugin VTKToMEDMem) +INSTALL(TARGETS DevelopedSurfacePlugin RUNTIME DESTINATION lib/paraview LIBRARY DESTINATION lib/paraview ARCHIVE DESTINATION lib/paraview) diff --git a/src/Plugins/DevelopedSurface/ParaViewPlugin/Resources/DevelopedSurfaceServer.xml b/src/Plugins/DevelopedSurface/ParaViewPlugin/Resources/DevelopedSurfaceServer.xml new file mode 100644 index 00000000..15570245 --- /dev/null +++ b/src/Plugins/DevelopedSurface/ParaViewPlugin/Resources/DevelopedSurfaceServer.xml @@ -0,0 +1,47 @@ + + + + + + + + + + + + + This property specifies the input to the Level Scalars filter. + + + + + + + + + + This property sets the parameters of cylinder used for slice. + + + + Specify if way used to develop theta is inverted or not relative to the reference one. By default no. The reference way used is those defined by the first cell sharing node 0. + + + + By default, node 0 theta parameter is used as starting point. This property allows to change this reference by applying an offset on it. Offset is expressed in degrees. + + + + + + + diff --git a/src/Plugins/DevelopedSurface/Test/CMakeLists.txt b/src/Plugins/DevelopedSurface/Test/CMakeLists.txt new file mode 100644 index 00000000..0ebbac38 --- /dev/null +++ b/src/Plugins/DevelopedSurface/Test/CMakeLists.txt @@ -0,0 +1,14 @@ +# Copyright (C) 2017 EDF R&D + +SET(TEMP_DIR "${CMAKE_CURRENT_BINARY_DIR}/Testing/Temporary") + +IF(NOT EXISTS ${TEMP_DIR}) + FILE(MAKE_DIRECTORY ${TEMP_DIR}) +ENDIF(NOT EXISTS ${TEMP_DIR}) + +SET(DEV_SURFACE_TESTS test_dev_surface2 test_dev_surface3) + +FOREACH(tfile ${DEV_SURFACE_TESTS}) + ADD_TEST(${tfile} ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/${tfile}.py ) + SET_TESTS_PROPERTIES(${tfile} PROPERTIES LABELS "PVS_ADD_ONS") +ENDFOREACH(tfile ${DEV_SURFACE_TESTS}) diff --git a/src/Plugins/DevelopedSurface/Test/test_dev_surface.py b/src/Plugins/DevelopedSurface/Test/test_dev_surface.py new file mode 100644 index 00000000..69317369 --- /dev/null +++ b/src/Plugins/DevelopedSurface/Test/test_dev_surface.py @@ -0,0 +1,162 @@ +# Copyright (C) 2017 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +#### import the simple module from the paraview +from paraview.simple import * +#### disable automatic camera reset on 'Show' +paraview.simple._DisableFirstRenderCameraReset() + +# create a new 'MED Reader' +multiTSmed = MEDReader(FileName='multiTS.med') +multiTSmed.AllArrays = ['TS0/Mesh/ComSup0/Pressure@@][@@P0'] +multiTSmed.AllTimeSteps = ['0000', '0001', '0002', '0003', '0004', '0005', '0006', '0007', '0008', '0009'] + +# get animation scene +animationScene1 = GetAnimationScene() + +# update animation scene based on data timesteps +animationScene1.UpdateAnimationUsingDataTimeSteps() + +# get active view +renderView1 = GetActiveViewOrCreate('RenderView') +# uncomment following to set a specific view size +# renderView1.ViewSize = [1499, 582] + +# show data in view +multiTSmedDisplay = Show(multiTSmed, renderView1) + +# trace defaults for the display properties. +multiTSmedDisplay.Representation = 'Surface' +multiTSmedDisplay.ColorArrayName = [None, ''] +multiTSmedDisplay.OSPRayScaleArray = 'FamilyIdNode' +multiTSmedDisplay.OSPRayScaleFunction = 'PiecewiseFunction' +multiTSmedDisplay.SelectOrientationVectors = 'FamilyIdNode' +multiTSmedDisplay.ScaleFactor = 0.07399989366531372 +multiTSmedDisplay.SelectScaleArray = 'FamilyIdNode' +multiTSmedDisplay.GlyphType = 'Arrow' +multiTSmedDisplay.GlyphTableIndexArray = 'FamilyIdNode' +multiTSmedDisplay.DataAxesGrid = 'GridAxesRepresentation' +multiTSmedDisplay.PolarAxes = 'PolarAxesRepresentation' +multiTSmedDisplay.ScalarOpacityUnitDistance = 0.017316274962626298 +multiTSmedDisplay.GaussianRadius = 0.03699994683265686 +multiTSmedDisplay.SetScaleArray = ['POINTS', 'FamilyIdNode'] +multiTSmedDisplay.ScaleTransferFunction = 'PiecewiseFunction' +multiTSmedDisplay.OpacityArray = ['POINTS', 'FamilyIdNode'] +multiTSmedDisplay.OpacityTransferFunction = 'PiecewiseFunction' +multiTSmedDisplay.InputVectors = [None, ''] +multiTSmedDisplay.SelectInputVectors = [None, ''] +multiTSmedDisplay.WriteLog = '' + +# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction' +multiTSmedDisplay.ScaleTransferFunction.Points = [0.0, 0.0, 0.5, 0.0, 1.1757813367477812e-38, 1.0, 0.5, 0.0] + +# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction' +multiTSmedDisplay.OpacityTransferFunction.Points = [0.0, 0.0, 0.5, 0.0, 1.1757813367477812e-38, 1.0, 0.5, 0.0] + +# reset view to fit data +renderView1.ResetCamera() + +# update the view to ensure updated data information +renderView1.Update() + +# create a new 'Developed Surface' +developedSurface1 = DevelopedSurface(Input=multiTSmed) +developedSurface1.SliceType = 'Cylinder' + +# init the 'Cylinder' selected for 'SliceType' +developedSurface1.SliceType.Center = [0.0, 0.0, 0.05000000074505806] +developedSurface1.SliceType.Radius = 0.3699994683265686 + +# Properties modified on developedSurface1.SliceType +developedSurface1.SliceType.Center = [0.0, 0.0, 0.05] +developedSurface1.SliceType.Axis = [0.0, 0.0, 1.0] +developedSurface1.SliceType.Radius = 0.07 + +# Properties modified on developedSurface1.SliceType +developedSurface1.SliceType.Center = [0.0, 0.0, 0.05] +developedSurface1.SliceType.Axis = [0.0, 0.0, 1.0] +developedSurface1.SliceType.Radius = 0.07 + +# show data in view +developedSurface1Display = Show(developedSurface1, renderView1) + +# trace defaults for the display properties. +developedSurface1Display.Representation = 'Surface' +developedSurface1Display.ColorArrayName = [None, ''] +developedSurface1Display.OSPRayScaleArray = 'FamilyIdNode' +developedSurface1Display.OSPRayScaleFunction = 'PiecewiseFunction' +developedSurface1Display.SelectOrientationVectors = 'FamilyIdNode' +developedSurface1Display.ScaleFactor = 0.043982297150257116 +developedSurface1Display.SelectScaleArray = 'FamilyIdNode' +developedSurface1Display.GlyphType = 'Arrow' +developedSurface1Display.GlyphTableIndexArray = 'FamilyIdNode' +developedSurface1Display.DataAxesGrid = 'GridAxesRepresentation' +developedSurface1Display.PolarAxes = 'PolarAxesRepresentation' +developedSurface1Display.GaussianRadius = 0.021991148575128558 +developedSurface1Display.SetScaleArray = ['POINTS', 'FamilyIdNode'] +developedSurface1Display.ScaleTransferFunction = 'PiecewiseFunction' +developedSurface1Display.OpacityArray = ['POINTS', 'FamilyIdNode'] +developedSurface1Display.OpacityTransferFunction = 'PiecewiseFunction' +developedSurface1Display.InputVectors = [None, ''] +developedSurface1Display.SelectInputVectors = [None, ''] +developedSurface1Display.WriteLog = '' + +# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction' +developedSurface1Display.ScaleTransferFunction.Points = [0.0, 0.0, 0.5, 0.0, 1.1757813367477812e-38, 1.0, 0.5, 0.0] + +# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction' +developedSurface1Display.OpacityTransferFunction.Points = [0.0, 0.0, 0.5, 0.0, 1.1757813367477812e-38, 1.0, 0.5, 0.0] + +# hide data in view +Hide(multiTSmed, renderView1) + +# update the view to ensure updated data information +renderView1.Update() + +#change interaction mode for render view +renderView1.InteractionMode = '2D' + +# toggle 3D widget visibility (only when running from the GUI) +Hide3DWidgets(proxy=developedSurface1.SliceType) + +# set scalar coloring +ColorBy(developedSurface1Display, ('CELLS', 'Pressure')) + +# rescale color and/or opacity maps used to include current data range +developedSurface1Display.RescaleTransferFunctionToDataRange(True, False) + +# show color bar/color legend +developedSurface1Display.SetScalarBarVisibility(renderView1, True) + +# get color transfer function/color map for 'Pressure' +pressureLUT = GetColorTransferFunction('Pressure') + +#### saving camera placements for all active views + +# current camera placement for renderView1 +renderView1.InteractionMode = '2D' +renderView1.CameraPosition = [0.18935662797765695, 0.01726656182167085, 2.08092363470839] +renderView1.CameraFocalPoint = [0.18935662797765695, 0.01726656182167085, 0.05000000074505806] +renderView1.CameraParallelScale = 0.16748564967020724 + +#### uncomment the following to render all views +# RenderAllViews() +# alternatively, if you want to write images, you can use SaveScreenshot(...). +Render() diff --git a/src/Plugins/DevelopedSurface/Test/test_dev_surface2.py b/src/Plugins/DevelopedSurface/Test/test_dev_surface2.py new file mode 100644 index 00000000..a3c35806 --- /dev/null +++ b/src/Plugins/DevelopedSurface/Test/test_dev_surface2.py @@ -0,0 +1,170 @@ +# Copyright (C) 2017 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +#### import the simple module from the paraview +from paraview.simple import * +from math import pi +TMPFileName="test2.med" + +#### disable automatic camera reset on 'Show' +paraview.simple._DisableFirstRenderCameraReset() + +# create a new 'Mandelbrot' +mandelbrot1 = Mandelbrot() + +# Properties modified on mandelbrot1 +mandelbrot1.WholeExtent = [0, 50, 0, 50, 0, 50] + +# get active view +renderView1 = GetActiveViewOrCreate('RenderView') +# uncomment following to set a specific view size +# renderView1.ViewSize = [1017, 317] + +# show data in view +mandelbrot1Display = Show(mandelbrot1, renderView1) + +# trace defaults for the display properties. +mandelbrot1Display.Representation = 'Outline' +mandelbrot1Display.ColorArrayName = ['POINTS', ''] +mandelbrot1Display.OSPRayScaleArray = 'Iterations' +mandelbrot1Display.OSPRayScaleFunction = 'PiecewiseFunction' +mandelbrot1Display.SelectOrientationVectors = 'Iterations' +mandelbrot1Display.ScaleFactor = 0.25 +mandelbrot1Display.SelectScaleArray = 'Iterations' +mandelbrot1Display.GlyphType = 'Arrow' +mandelbrot1Display.GlyphTableIndexArray = 'Iterations' +mandelbrot1Display.DataAxesGrid = 'GridAxesRepresentation' +mandelbrot1Display.PolarAxes = 'PolarAxesRepresentation' +mandelbrot1Display.ScalarOpacityUnitDistance = 0.08124038404635964 +mandelbrot1Display.Slice = 25 +mandelbrot1Display.GaussianRadius = 0.125 +mandelbrot1Display.SetScaleArray = ['POINTS', 'Iterations'] +mandelbrot1Display.ScaleTransferFunction = 'PiecewiseFunction' +mandelbrot1Display.OpacityArray = ['POINTS', 'Iterations'] +mandelbrot1Display.OpacityTransferFunction = 'PiecewiseFunction' +mandelbrot1Display.InputVectors = [None, ''] +mandelbrot1Display.SelectInputVectors = [None, ''] +mandelbrot1Display.WriteLog = '' + +# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction' +mandelbrot1Display.ScaleTransferFunction.Points = [1.0, 0.0, 0.5, 0.0, 100.0, 1.0, 0.5, 0.0] + +# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction' +mandelbrot1Display.OpacityTransferFunction.Points = [1.0, 0.0, 0.5, 0.0, 100.0, 1.0, 0.5, 0.0] + +# reset view to fit data +renderView1.ResetCamera() + +# update the view to ensure updated data information +renderView1.Update() + +# create a new 'Developed Surface' +developedSurface1 = DevelopedSurface(Input=mandelbrot1) +developedSurface1.SliceType = 'Cylinder' + +# init the 'Cylinder' selected for 'SliceType' +developedSurface1.SliceType.Center = [-0.5, 0.0, 1.0] +developedSurface1.SliceType.Radius = 0.5 #1.25 + +# show data in view +developedSurface1Display = Show(developedSurface1, renderView1) + +# get color transfer function/color map for 'Iterations' +iterationsLUT = GetColorTransferFunction('Iterations') + +# trace defaults for the display properties. +developedSurface1Display.Representation = 'Surface' +developedSurface1Display.ColorArrayName = ['POINTS', 'Iterations'] +developedSurface1Display.LookupTable = iterationsLUT +developedSurface1Display.OSPRayScaleArray = 'Iterations' +developedSurface1Display.OSPRayScaleFunction = 'PiecewiseFunction' +developedSurface1Display.SelectOrientationVectors = 'Iterations' +developedSurface1Display.ScaleFactor = 0.7853981633974483 +developedSurface1Display.SelectScaleArray = 'Iterations' +developedSurface1Display.GlyphType = 'Arrow' +developedSurface1Display.GlyphTableIndexArray = 'Iterations' +developedSurface1Display.DataAxesGrid = 'GridAxesRepresentation' +developedSurface1Display.PolarAxes = 'PolarAxesRepresentation' +developedSurface1Display.GaussianRadius = 0.39269908169872414 +developedSurface1Display.SetScaleArray = ['POINTS', 'Iterations'] +developedSurface1Display.ScaleTransferFunction = 'PiecewiseFunction' +developedSurface1Display.OpacityArray = ['POINTS', 'Iterations'] +developedSurface1Display.OpacityTransferFunction = 'PiecewiseFunction' +developedSurface1Display.InputVectors = [None, ''] +developedSurface1Display.SelectInputVectors = [None, ''] +developedSurface1Display.WriteLog = '' + +# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction' +developedSurface1Display.ScaleTransferFunction.Points = [1.0, 0.0, 0.5, 0.0, 100.0, 1.0, 0.5, 0.0] + +# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction' +developedSurface1Display.OpacityTransferFunction.Points = [1.0, 0.0, 0.5, 0.0, 100.0, 1.0, 0.5, 0.0] + +# hide data in view +Hide(mandelbrot1, renderView1) + +# show color bar/color legend +developedSurface1Display.SetScalarBarVisibility(renderView1, True) + +# update the view to ensure updated data information +renderView1.Update() + +# toggle 3D widget visibility (only when running from the GUI) +Hide3DWidgets(proxy=developedSurface1.SliceType) + +#### saving camera placements for all active views + +# current camera placement for renderView1 +renderView1.CameraPosition = [4.090024784500779, -0.15919161102314858, 7.485304552729019] +renderView1.CameraFocalPoint = [4.090024784500779, -0.15919161102314858, 1.0] +renderView1.CameraParallelScale = 2.03100960115899 + +#### uncomment the following to render all views +# RenderAllViews() +# alternatively, if you want to write images, you can use SaveScreenshot(...). + +mand=servermanager.Fetch(mandelbrot1) +axisId=1 +high_out=mand.GetSpacing()[axisId]*(mand.GetExtent()[2*axisId+1]-mand.GetExtent()[2*axisId+0]) + +vtp=servermanager.Fetch(developedSurface1) +arr=vtp.GetPointData().GetArray(0) +assert(arr.GetName()=="Iterations") +a,b=arr.GetRange() +assert(a>=1 and a<=2) +assert(b==100.) +SaveData(TMPFileName, proxy=developedSurface1) +from MEDLoader import * + +mm=MEDFileMesh.New(TMPFileName) +m0=mm[0] +area=m0.getMeasureField(True).getArray().accumulate()[0] + +zeResu0=area/high_out/developedSurface1.SliceType.Radius +assert(abs(zeResu0-2*pi)<1e-5) + +fs=MEDFileFields(TMPFileName) +f=fs["Iterations"][0].field(mm) +nodeIds=f.getArray().convertToDblArr().findIdsInRange(99.,101.) +cellIds=m0.getCellIdsLyingOnNodes(nodeIds,True) +zeResu1=m0[cellIds].getMeasureField(True).getArray().accumulate()[0] + +assert(abs(zeResu1-1.1427)<1e-2) + diff --git a/src/Plugins/DevelopedSurface/Test/test_dev_surface3.py b/src/Plugins/DevelopedSurface/Test/test_dev_surface3.py new file mode 100644 index 00000000..231de21c --- /dev/null +++ b/src/Plugins/DevelopedSurface/Test/test_dev_surface3.py @@ -0,0 +1,166 @@ +# Copyright (C) 2017 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +#### import the simple module from the paraview +from paraview.simple import * +from math import pi +TMPFileName="test3.med" + +#### disable automatic camera reset on 'Show' +paraview.simple._DisableFirstRenderCameraReset() + +# create a new 'Mandelbrot' +mandelbrot1 = Mandelbrot() + +# Properties modified on mandelbrot1 +mandelbrot1.WholeExtent = [0, 50, 0, 50, 0, 50] + +# get active view +renderView1 = GetActiveViewOrCreate('RenderView') +# uncomment following to set a specific view size +# renderView1.ViewSize = [1017, 317] + +# show data in view +mandelbrot1Display = Show(mandelbrot1, renderView1) + +# trace defaults for the display properties. +mandelbrot1Display.Representation = 'Outline' +mandelbrot1Display.ColorArrayName = ['POINTS', ''] +mandelbrot1Display.OSPRayScaleArray = 'Iterations' +mandelbrot1Display.OSPRayScaleFunction = 'PiecewiseFunction' +mandelbrot1Display.SelectOrientationVectors = 'Iterations' +mandelbrot1Display.ScaleFactor = 0.25 +mandelbrot1Display.SelectScaleArray = 'Iterations' +mandelbrot1Display.GlyphType = 'Arrow' +mandelbrot1Display.GlyphTableIndexArray = 'Iterations' +mandelbrot1Display.DataAxesGrid = 'GridAxesRepresentation' +mandelbrot1Display.PolarAxes = 'PolarAxesRepresentation' +mandelbrot1Display.ScalarOpacityUnitDistance = 0.08124038404635964 +mandelbrot1Display.Slice = 25 +mandelbrot1Display.GaussianRadius = 0.125 +mandelbrot1Display.SetScaleArray = ['POINTS', 'Iterations'] +mandelbrot1Display.ScaleTransferFunction = 'PiecewiseFunction' +mandelbrot1Display.OpacityArray = ['POINTS', 'Iterations'] +mandelbrot1Display.OpacityTransferFunction = 'PiecewiseFunction' +mandelbrot1Display.InputVectors = [None, ''] +mandelbrot1Display.SelectInputVectors = [None, ''] +mandelbrot1Display.WriteLog = '' + +# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction' +mandelbrot1Display.ScaleTransferFunction.Points = [1.0, 0.0, 0.5, 0.0, 100.0, 1.0, 0.5, 0.0] + +# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction' +mandelbrot1Display.OpacityTransferFunction.Points = [1.0, 0.0, 0.5, 0.0, 100.0, 1.0, 0.5, 0.0] + +# reset view to fit data +renderView1.ResetCamera() + +# update the view to ensure updated data information +renderView1.Update() + +# create a new 'Developed Surface' +developedSurface1 = DevelopedSurface(Input=mandelbrot1) +developedSurface1.SliceType = 'Cylinder' + +# init the 'Cylinder' selected for 'SliceType' +developedSurface1.SliceType.Center = [-0.5, 0.0, 1.0] +developedSurface1.SliceType.Radius = 0.5 #1.25 +developedSurface1.SliceType.Axis = [-0.5065630563269753, -0.6288876685363318, -0.5898255422814533] + +# show data in view +developedSurface1Display = Show(developedSurface1, renderView1) + +# get color transfer function/color map for 'Iterations' +iterationsLUT = GetColorTransferFunction('Iterations') + +# trace defaults for the display properties. +developedSurface1Display.Representation = 'Surface' +developedSurface1Display.ColorArrayName = ['POINTS', 'Iterations'] +developedSurface1Display.LookupTable = iterationsLUT +developedSurface1Display.OSPRayScaleArray = 'Iterations' +developedSurface1Display.OSPRayScaleFunction = 'PiecewiseFunction' +developedSurface1Display.SelectOrientationVectors = 'Iterations' +developedSurface1Display.ScaleFactor = 0.7853981633974483 +developedSurface1Display.SelectScaleArray = 'Iterations' +developedSurface1Display.GlyphType = 'Arrow' +developedSurface1Display.GlyphTableIndexArray = 'Iterations' +developedSurface1Display.DataAxesGrid = 'GridAxesRepresentation' +developedSurface1Display.PolarAxes = 'PolarAxesRepresentation' +developedSurface1Display.GaussianRadius = 0.39269908169872414 +developedSurface1Display.SetScaleArray = ['POINTS', 'Iterations'] +developedSurface1Display.ScaleTransferFunction = 'PiecewiseFunction' +developedSurface1Display.OpacityArray = ['POINTS', 'Iterations'] +developedSurface1Display.OpacityTransferFunction = 'PiecewiseFunction' +developedSurface1Display.InputVectors = [None, ''] +developedSurface1Display.SelectInputVectors = [None, ''] +developedSurface1Display.WriteLog = '' + +# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction' +developedSurface1Display.ScaleTransferFunction.Points = [1.0, 0.0, 0.5, 0.0, 100.0, 1.0, 0.5, 0.0] + +# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction' +developedSurface1Display.OpacityTransferFunction.Points = [1.0, 0.0, 0.5, 0.0, 100.0, 1.0, 0.5, 0.0] + +# hide data in view +Hide(mandelbrot1, renderView1) + +# show color bar/color legend +developedSurface1Display.SetScalarBarVisibility(renderView1, True) + +# update the view to ensure updated data information +renderView1.Update() + +# toggle 3D widget visibility (only when running from the GUI) +Hide3DWidgets(proxy=developedSurface1.SliceType) + +#### saving camera placements for all active views + +# current camera placement for renderView1 +renderView1.CameraPosition = [4.090024784500779, -0.15919161102314858, 7.485304552729019] +renderView1.CameraFocalPoint = [4.090024784500779, -0.15919161102314858, 1.0] +renderView1.CameraParallelScale = 2.03100960115899 + +#### uncomment the following to render all views +# RenderAllViews() +# alternatively, if you want to write images, you can use SaveScreenshot(...). + + +vtp=servermanager.Fetch(developedSurface1) +arr=vtp.GetPointData().GetArray(0) +assert(arr.GetName()=="Iterations") +a,b=arr.GetRange() +assert(a>=1 and a<=2) +assert(b==100.) +SaveData(TMPFileName, proxy=developedSurface1) +from MEDLoader import * + +mm=MEDFileMesh.New(TMPFileName) +m0=mm[0] +area=m0.getMeasureField(True).getArray().accumulate()[0] + + +fs=MEDFileFields(TMPFileName) +f=fs["Iterations"][0].field(mm) +nodeIds=f.getArray().convertToDblArr().findIdsInRange(99.,101.) +cellIds=m0.getCellIdsLyingOnNodes(nodeIds,True) +zeResu1=m0[cellIds].getMeasureField(True).getArray().accumulate()[0] + +assert(abs(zeResu1-1.3564)<1e-2) + diff --git a/src/Plugins/DevelopedSurface/pp.py b/src/Plugins/DevelopedSurface/pp.py new file mode 100644 index 00000000..cee04bf5 --- /dev/null +++ b/src/Plugins/DevelopedSurface/pp.py @@ -0,0 +1,107 @@ +# Copyright (C) 2017 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +from MEDLoader import * +from math import sqrt +import numpy as np +import scipy +import scipy.sparse.linalg + +def f0(sample,p,v,r): + d=sample-p + return d.magnitude()[0]-r + +def f(sample,p,v,r): + d=sample-p + l=d-DataArrayDouble.Dot(d,v)[0]*v + return l.magnitude()[0]-r + +def f2(sample,zev,ff): + p=zev[0,:3] ; v=zev[0,3:6] ; r=zev[0,6] + return ff(sample,p,v,r) + +def df(sample,zev,varid,ff): + eps=0.0001 + zev=zev[:] + zev2=zev[:] + zev[0,varid]+=eps ; zev2[0,varid]-=eps + return (f2(sample,zev,ff)-f2(sample,zev2,ff))/(2*eps) + +#def df2(sample,p,v,r,varid): +# zev=DataArrayDouble.Meld([p_s,v_s,DataArrayDouble([r_s])]) +# return df(sample,zev,varid) + +def jacob(sample,p,v,r,ff): + zev=DataArrayDouble.Meld([p_s,v_s,DataArrayDouble([r_s])]) + return DataArrayDouble([df(sample,zev,i,ff) for i in range(7)],1,7) + +def jacob0(sample,p,v,r): + zev=DataArrayDouble.Meld([p_s,v_s,DataArrayDouble([r_s])]) + return DataArrayDouble([df0(sample,zev,i) for i in range(7)],1,7) + +mm=MEDFileMesh.New("example2.med") +c=mm.getCoords() +p_s=DataArrayDouble(c.accumulate(),1,3)/float(len(c)) +v_s=DataArrayDouble([1,0,0],1,3) +o=DataArrayDouble(c.getMinMaxPerComponent(),3,2) +o0,o1=o.explodeComponents() +o=o1-o0 +o.abs() +r_s=o.getMaxValue()[0]/2. +# +r_s=0.215598 +p_s=DataArrayDouble([0.,0.,0.],1,3) +v_s=DataArrayDouble([0.,0.,1.],1,3) +# +r_s=0.2 +p_s=DataArrayDouble([1.,1.,1.],1,3) +v_s=DataArrayDouble([1.,1.,1.],1,3) +#probes=[0,979,1167,2467,2862,3706,3819] +probes=[1000]+[c[:,i].getMaxValue()[1] for i in range(3)]+[c[:,i].getMinValue()[1] for i in range(3)] + +p=p_s ; v=v_s ; r=r_s +for ii in range(1): + mat=DataArrayDouble.Aggregate([jacob(c[probes[0]],p,v,r,f0)]+[jacob(c[probe],p,v,r,f) for probe in probes[1:]])#+[DataArrayDouble([0,0,0,2*v[0,0],2*v[0,1],2*v[0,2],0],1,7)] + y=DataArrayDouble([f0(c[probes[0]],p,v,r)]+[f(c[probe],p,v,r) for probe in probes[1:]]) + #y.pushBackSilent(v.magnitude()[0]-1.) + delta2=y[:] ; delta2.abs() + if delta2.getMaxValueInArray()<1e-5: + print("finished") + break; + mat=scipy.matrix(mat.toNumPyArray()) + print ii,np.linalg.cond(mat) + y=scipy.matrix(y.toNumPyArray()) + y=y.transpose() + delta=np.linalg.solve(mat,-y) + #delta=scipy.sparse.linalg.gmres(mat,-y)[0] # cg + delta=DataArrayDouble(delta) ; delta.rearrange(7) + #p+=delta[0,:3] + #v+=delta[0,3:6] + #v/=v.magnitude()[0] + #r+=delta[0,6] + print ii,delta.getValues()#,yy.transpose().tolist() + + pass + +mat1=mat[[0,1,2,6]] +mat1=mat1.transpose() +mat1=mat1[[0,1,2,6]] +mat1=mat1.transpose() +mm=np.linalg.solve(mat1,y[[0,1,2,6]]) diff --git a/src/Plugins/DevelopedSurface/pp1.py b/src/Plugins/DevelopedSurface/pp1.py new file mode 100644 index 00000000..54393116 --- /dev/null +++ b/src/Plugins/DevelopedSurface/pp1.py @@ -0,0 +1,27 @@ +# Copyright (C) 2017 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +from MEDLoader import * +from math import pi +mm=MEDFileMesh.New("example2.med") +c=mm.getCoords() +MEDCouplingUMesh.Rotate3DAlg([0,0,0],[1,1,0],pi/4,c) +c+=[3,2,0] +mm.write("example2_2.med",2) diff --git a/src/Plugins/DevelopedSurface/pp2.py b/src/Plugins/DevelopedSurface/pp2.py new file mode 100644 index 00000000..6c8d08e2 --- /dev/null +++ b/src/Plugins/DevelopedSurface/pp2.py @@ -0,0 +1,121 @@ +# Copyright (C) 2017 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +from MEDLoader import * +from math import pi,asin,acos,atan2 +import numpy as np + +invertThetaInc=False + +mm=MEDFileMesh.New("example2_2.med") +c=mm.getCoords() + +v=DataArrayDouble([0.70710678118617443,-0.70710678118572623,1],1,3) +v/=v.magnitude()[0] +vRef=DataArrayDouble([0,0,1],1,3) +vRot=DataArrayDouble.CrossProduct(v,vRef) +ang=asin(vRot.magnitude()[0]) +MEDCouplingUMesh.Rotate3DAlg([0,0,0],vRot.getValues(),ang,c) + +c2=c[:,[0,1]] +# +a,b=c2.findCommonTuples(1e-5) +pool=c2[a[b[:-1]]] +probes=[pool[:,0].getMaxValue()[1],pool[:,1].getMaxValue()[1],pool[:,0].getMinValue()[1]] +center,radius,_=pool[probes].asArcOfCircle() +assert(((c2-center).magnitude()).isUniform(radius,1e-5)) +theta=((c2-center)/radius) +theta=DataArrayDouble([atan2(y,x) for x,y in theta]) +zeTheta=theta-theta[0] + +mm.write("tmp.med",2) +########### +c_cyl=DataArrayDouble.fromCartToCyl(c-(center+[0.])) +assert(c_cyl[:,0].isUniform(radius,1e-5)) +tetha0=c_cyl[0,1] +c_cyl[:,1]-=tetha0 +c_cyl_2=c_cyl[:,1] ; c_cyl_2.abs() +c_cyl[c_cyl_2.findIdsInRange(0.,1e-6),1]=0. +######## +m0=mm[0] +c0s=m0.getCellIdsLyingOnNodes([0],False) +assert(len(c0s)!=0) +tmp=c_cyl[m0.getNodeIdsOfCell(c0s[0]),1].getMaxAbsValueInArray() +a=tmp>0. +if a^(not invertThetaInc): + c_cyl[:,1]=-c_cyl[:,1] + pass +# +c_cyl_post=c_cyl[:,1] +m0.convertAllToPoly() +tmp=MEDCoupling1DGTUMesh(m0) ; c=tmp.getNodalConnectivity() ; ci=tmp.getNodalConnectivityIndex() +for elt in c_cyl_post: + if float(elt)<0: + elt[:]=float(elt)+2*pi + +nbCells=m0.getNumberOfCells() +nbNodes=m0.getNumberOfNodes() +newCoo=DataArrayInt(0,1) +cellsWithPb=[] +newConn=[] +for i in xrange(nbCells): + tmp=c_cyl_post[c[ci[i]:ci[i+1]]] + if tmp.getMaxValueInArray()-tmp.getMinValueInArray()>pi: + cellsWithPb.append(i) + # dup of low val + newLocConn=[] + for elt in c[ci[i]:ci[i+1]]: + if float(c_cyl_post[elt])<=pi: + newLocConn.append(nbNodes+len(newCoo)) + newCoo.pushBackSilent(int(elt)) + pass + else: + newLocConn.append(int(elt)) + pass + pass + newConn.append(newLocConn[:]) + pass + + +c_cyl_post=DataArrayDouble.Aggregate([c_cyl_post,c_cyl_post[newCoo]+2*pi]) +z_post=DataArrayDouble.Aggregate([c_cyl[:,2],c_cyl[newCoo,2]]) +newCoords=DataArrayDouble.Meld(radius*c_cyl_post,z_post) + +z=DataArrayInt([len(elt) for elt in newConn]) +z.computeOffsetsFull() +# +cellsWithPbMesh=MEDCoupling1DGTUMesh("",NORM_POLYGON) +cellsWithPbMesh.setNodalConnectivity(DataArrayInt(sum(newConn,[])),z) +cellsWithPbMesh.setCoords(newCoords) +# +ko_part_ids=DataArrayInt(cellsWithPb) +part_ids=ko_part_ids.buildComplement(nbCells) +part=m0[part_ids] +part.setCoords(newCoords) +whole=MEDCouplingUMesh.MergeUMeshesOnSameCoords([part,cellsWithPbMesh.buildUnstructured()]) + +o2n=DataArrayInt.Aggregate([part_ids,ko_part_ids]) +whole=whole[o2n.invertArrayO2N2N2O(len(o2n))] +whole.setName("mesh") +WriteMesh("tmp3.med",whole,True) + +#sk=m0.computeSkin().computeFetchedNodeIds() +#sk_part=part.computeSkin().computeFetchedNodeIds() +#m0[DataArrayInt(cellsWithPb)].computeFetchedNodeIds() diff --git a/src/Plugins/DevelopedSurface/test.py b/src/Plugins/DevelopedSurface/test.py new file mode 100644 index 00000000..e1e3b23f --- /dev/null +++ b/src/Plugins/DevelopedSurface/test.py @@ -0,0 +1,139 @@ +# Copyright (C) 2017 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +#### import the simple module from the paraview +from paraview.simple import * +#### disable automatic camera reset on 'Show' +paraview.simple._DisableFirstRenderCameraReset() + +# create a new 'Wavelet' +wavelet1 = Wavelet() + +# Properties modified on wavelet1 +wavelet1.WholeExtent = [0, 10, 0, 10, 0, 10] + +# get active view +renderView1 = GetActiveViewOrCreate('RenderView') +# uncomment following to set a specific view size +# renderView1.ViewSize = [1168, 582] + +# show data in view +wavelet1Display = Show(wavelet1, renderView1) + +# trace defaults for the display properties. +wavelet1Display.Representation = 'Outline' +wavelet1Display.ColorArrayName = ['POINTS', ''] +wavelet1Display.OSPRayScaleArray = 'RTData' +wavelet1Display.OSPRayScaleFunction = 'PiecewiseFunction' +wavelet1Display.SelectOrientationVectors = 'None' +wavelet1Display.SelectScaleArray = 'RTData' +wavelet1Display.GlyphType = 'Arrow' +wavelet1Display.GlyphTableIndexArray = 'RTData' +wavelet1Display.DataAxesGrid = 'GridAxesRepresentation' +wavelet1Display.PolarAxes = 'PolarAxesRepresentation' +wavelet1Display.ScalarOpacityUnitDistance = 1.7320508075688779 +wavelet1Display.Slice = 5 +wavelet1Display.GaussianRadius = 0.5 +wavelet1Display.SetScaleArray = ['POINTS', 'RTData'] +wavelet1Display.ScaleTransferFunction = 'PiecewiseFunction' +wavelet1Display.OpacityArray = ['POINTS', 'RTData'] +wavelet1Display.OpacityTransferFunction = 'PiecewiseFunction' +wavelet1Display.InputVectors = [None, ''] +wavelet1Display.SelectInputVectors = [None, ''] +wavelet1Display.WriteLog = '' + +# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction' +wavelet1Display.ScaleTransferFunction.Points = [-16.577068328857422, 0.0, 0.5, 0.0, 260.0, 1.0, 0.5, 0.0] + +# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction' +wavelet1Display.OpacityTransferFunction.Points = [-16.577068328857422, 0.0, 0.5, 0.0, 260.0, 1.0, 0.5, 0.0] + +# reset view to fit data +renderView1.ResetCamera() + +# update the view to ensure updated data information +renderView1.Update() + +# create a new 'Slice' +slice1 = Slice(Input=wavelet1) +slice1.SliceType = 'Plane' +slice1.SliceOffsetValues = [0.0] + +# init the 'Plane' selected for 'SliceType' +slice1.SliceType.Origin = [5.0, 5.0, 5.0] + +# toggle 3D widget visibility (only when running from the GUI) +Show3DWidgets(proxy=slice1.SliceType) + +# Properties modified on slice1 +slice1.SliceType = 'Cylinder' + +# show data in view +slice1Display = Show(slice1, renderView1) + +# get color transfer function/color map for 'RTData' +rTDataLUT = GetColorTransferFunction('RTData') + +# trace defaults for the display properties. +slice1Display.Representation = 'Surface' +slice1Display.ColorArrayName = ['POINTS', 'RTData'] +slice1Display.LookupTable = rTDataLUT +slice1Display.OSPRayScaleArray = 'RTData' +slice1Display.OSPRayScaleFunction = 'PiecewiseFunction' +slice1Display.SelectOrientationVectors = 'None' +slice1Display.SelectScaleArray = 'RTData' +slice1Display.GlyphType = 'Arrow' +slice1Display.GlyphTableIndexArray = 'RTData' +slice1Display.DataAxesGrid = 'GridAxesRepresentation' +slice1Display.PolarAxes = 'PolarAxesRepresentation' +slice1Display.GaussianRadius = 0.5 +slice1Display.SetScaleArray = ['POINTS', 'RTData'] +slice1Display.ScaleTransferFunction = 'PiecewiseFunction' +slice1Display.OpacityArray = ['POINTS', 'RTData'] +slice1Display.OpacityTransferFunction = 'PiecewiseFunction' +slice1Display.InputVectors = [None, ''] +slice1Display.SelectInputVectors = [None, ''] +slice1Display.WriteLog = '' + +# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction' +slice1Display.ScaleTransferFunction.Points = [-14.761041641235352, 0.0, 0.5, 0.0, 232.8310546875, 1.0, 0.5, 0.0] + +# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction' +slice1Display.OpacityTransferFunction.Points = [-14.761041641235352, 0.0, 0.5, 0.0, 232.8310546875, 1.0, 0.5, 0.0] + +# show color bar/color legend +slice1Display.SetScalarBarVisibility(renderView1, True) + +# update the view to ensure updated data information +renderView1.Update() + +# create a new 'Developed Surface' +developedSurface1 = DevelopedSurface(Input=slice1) + +#### saving camera placements for all active views + +# current camera placement for renderView1 +renderView1.CameraPosition = [5.0, 5.0, 38.46065214951232] +renderView1.CameraFocalPoint = [5.0, 5.0, 5.0] +renderView1.CameraParallelScale = 8.660254037844387 + +#### uncomment the following to render all views +# RenderAllViews() +# alternatively, if you want to write images, you can use SaveScreenshot(...). diff --git a/src/Plugins/GaussToCell/CMakeLists.txt b/src/Plugins/GaussToCell/CMakeLists.txt new file mode 100644 index 00000000..aa9a784b --- /dev/null +++ b/src/Plugins/GaussToCell/CMakeLists.txt @@ -0,0 +1,70 @@ +# Copyright (C) 2018 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +PROJECT(VoroGauss) +CMAKE_MINIMUM_REQUIRED(VERSION 2.8) + +CMAKE_POLICY(SET CMP0003 NEW) +IF(${CMAKE_VERSION} VERSION_GREATER "3.0.0") + CMAKE_POLICY(SET CMP0022 OLD) + CMAKE_POLICY(SET CMP0023 OLD) +ENDIF() + +SET(MED_READER_VERSION "0.0.0") + +# Common CMake macros +# =================== +SET(CONFIGURATION_ROOT_DIR $ENV{CONFIGURATION_ROOT_DIR} CACHE PATH "Path to the Salome CMake configuration files") +IF(EXISTS ${CONFIGURATION_ROOT_DIR}) + LIST(APPEND CMAKE_MODULE_PATH "${CONFIGURATION_ROOT_DIR}/cmake") + INCLUDE(SalomeMacros) +ELSE() + MESSAGE(FATAL_ERROR "We absolutely need the Salome CMake configuration files, please define CONFIGURATION_ROOT_DIR !") +ENDIF() +FIND_PACKAGE(SalomePythonInterp REQUIRED) +FIND_PACKAGE(SalomePythonLibs REQUIRED) + +FIND_PACKAGE(ParaView REQUIRED) +IF(NOT ParaView_FOUND) + MESSAGE(FATAL_ERROR "Please locate ParaView." ) +ENDIF(NOT ParaView_FOUND) +INCLUDE(${PARAVIEW_USE_FILE}) + +SET(MEDCOUPLING_ROOT_DIR $ENV{MEDCOUPLING_ROOT_DIR} CACHE PATH "MEDCOUPLING_ROOT_DIR") +LIST(APPEND CMAKE_MODULE_PATH "${MEDCOUPLING_ROOT_DIR}/cmake_files") +FIND_PACKAGE(SalomeMEDCoupling REQUIRED) + +OPTION(BUILD_SHARED_LIBS "Build with shared libraries." ${VTK_BUILD_SHARED_LIBS}) + +SET(VTK_INSTALL_RUNTIME_DIR lib) +SET(VTK_INSTALL_LIBRARY_DIR lib) +SET(VTK_INSTALL_ARCHIVE_DIR lib) + +PV_PROCESS_MODULES() + +INCLUDE_DIRECTORIES( + ${MEDCOUPLING_INCLUDE_DIRS} +) + +ADD_PARAVIEW_PLUGIN(GaussToCellPlugin "1.0" + SERVER_MANAGER_SOURCES ${SM_SRCS} vtkGaussToCell.cxx + SERVER_MANAGER_XML GaussToCellServer.xml) +TARGET_LINK_LIBRARIES(GaussToCellPlugin ${MEDCoupling_medcoupling}) +INSTALL(TARGETS GaussToCellPlugin RUNTIME DESTINATION lib/paraview LIBRARY DESTINATION lib/paraview ARCHIVE DESTINATION lib/paraview) diff --git a/src/Plugins/GaussToCell/GaussToCellServer.xml b/src/Plugins/GaussToCell/GaussToCellServer.xml new file mode 100644 index 00000000..162081fb --- /dev/null +++ b/src/Plugins/GaussToCell/GaussToCellServer.xml @@ -0,0 +1,42 @@ + + + + + + + + + + + + + This property specifies the input to the Level Scalars filter. + + + + + Foreach field on Gauss Points : computes a cell field lying of input mesh that is the average of the values associated to the set of Gauss points for each cell. + + + + Foreach field on Gauss Points : computes a cell field lying of input mesh that is the max of the values associated to the set of Gauss points for each cell. + + + + Foreach field on Gauss Points : computes a cell field lying of input mesh that is the min of the values associated to the set of Gauss points for each cell. + + + + + + + diff --git a/src/Plugins/GaussToCell/PG_3D.med b/src/Plugins/GaussToCell/PG_3D.med new file mode 100644 index 00000000..9c7e51e7 Binary files /dev/null and b/src/Plugins/GaussToCell/PG_3D.med differ diff --git a/src/Plugins/GaussToCell/TestCase.py b/src/Plugins/GaussToCell/TestCase.py new file mode 100644 index 00000000..7ebb27f5 --- /dev/null +++ b/src/Plugins/GaussToCell/TestCase.py @@ -0,0 +1,54 @@ +from MEDLoader import * + +fname="VoroGauss1.med" +meshName="mesh" +mm=MEDFileUMesh() +coords=DataArrayDouble([0,0, 1,0, 2,0, 3,0, 4,0, 5,0, 0,1, 1,1, 2,1, 0,2, 1,2, 3,1, 4,1],13,2) +m0=MEDCouplingUMesh(meshName,2) +m0.setCoords(coords) +m0.allocateCells() +m0.insertNextCell(NORM_TRI3,[2,3,8]) +m0.insertNextCell(NORM_TRI3,[3,4,11]) +m0.insertNextCell(NORM_TRI3,[4,5,12]) +m0.insertNextCell(NORM_TRI3,[6,7,9]) +m0.insertNextCell(NORM_TRI3,[7,8,10]) +m0.insertNextCell(NORM_QUAD4,[0,1,7,6]) +m0.insertNextCell(NORM_QUAD4,[1,2,8,7]) +mm[0]=m0 +m1=MEDCouplingUMesh(meshName,1) +m1.setCoords(coords) +m1.allocateCells() +m1.insertNextCell(NORM_SEG2,[0,1]) +m1.insertNextCell(NORM_SEG2,[1,2]) +m1.insertNextCell(NORM_SEG2,[2,3]) +m1.insertNextCell(NORM_SEG2,[3,4]) +m1.insertNextCell(NORM_SEG2,[4,5]) +mm[-1]=m1 +mm.setFamilyFieldArr(0,DataArrayInt([-1,-1,-2,-3,-3,-1,-3])) +mm.setFamilyFieldArr(-1,DataArrayInt([-1,-4,-4,-4,-1])) +for i in [-1,-2,-3,-4]: + mm.setFamilyId("Fam_%d"%i,i) + mm.setFamiliesOnGroup("G%d"%(abs(i)),["Fam_%d"%i]) + pass +mm.write(fname,2) +# +f0=MEDCouplingFieldDouble(ON_GAUSS_PT) +f0.setMesh(m0) +f0.setName("MyFieldPG") ; f0.setMesh(m0) +f0.setGaussLocalizationOnType(NORM_TRI3,[0,0, 1,0, 0,1],[0.1,0.1, 0.8,0.1, 0.1,0.8],[0.3,0.3,0.4]) +f0.setGaussLocalizationOnType(NORM_QUAD4,[-1,-1, 1,-1, 1,1, -1,1],[-0.57735,-0.57735,0.57735,-0.57735,0.57735,0.57735,-0.57735,0.57735],[0.25,0.25,0.25,0.25]) +arr=DataArrayDouble(f0.getNumberOfTuplesExpected()) ; arr.iota() +arr=DataArrayDouble.Meld(arr,arr) +arr.setInfoOnComponents(["comp0","comp1"]) +f0.setArray(arr) +WriteFieldUsingAlreadyWrittenMesh(fname,f0) +# +f1=MEDCouplingFieldDouble(ON_CELLS) +f1.setMesh(m0) +f1.setName("MyFieldCell") ; f1.setMesh(m0) +arr=DataArrayDouble(f1.getNumberOfTuplesExpected()) ; arr.iota() +arr=DataArrayDouble.Meld(arr,arr) +arr.setInfoOnComponents(["comp2","comp3"]) +f1.setArray(arr) +WriteFieldUsingAlreadyWrittenMesh(fname,f1) + diff --git a/src/Plugins/GaussToCell/testMEDReader14.med b/src/Plugins/GaussToCell/testMEDReader14.med new file mode 100644 index 00000000..c7ff867f Binary files /dev/null and b/src/Plugins/GaussToCell/testMEDReader14.med differ diff --git a/src/Plugins/GaussToCell/vtkGaussToCell.cxx b/src/Plugins/GaussToCell/vtkGaussToCell.cxx new file mode 100644 index 00000000..038b0b80 --- /dev/null +++ b/src/Plugins/GaussToCell/vtkGaussToCell.cxx @@ -0,0 +1,388 @@ +// Copyright (C) 2018 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 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 +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay (EDF R&D) + +#include "vtkGaussToCell.h" + +#include "vtkAdjacentVertexIterator.h" +#include "vtkIntArray.h" +#include "vtkCellData.h" +#include "vtkPointData.h" +#include "vtkCellType.h" +#include "vtkCell.h" +#include "vtkCellArray.h" +#include "vtkIdTypeArray.h" + +#include "vtkStreamingDemandDrivenPipeline.h" +#include "vtkUnstructuredGrid.h" +#include "vtkMultiBlockDataSet.h" + +#include "vtkInformationStringKey.h" +#include "vtkAlgorithmOutput.h" +#include "vtkObjectFactory.h" +#include "vtkMutableDirectedGraph.h" +#include "vtkMultiBlockDataSet.h" +#include "vtkDataSet.h" +#include "vtkInformationVector.h" +#include "vtkInformation.h" +#include "vtkDataArraySelection.h" +#include "vtkTimeStamp.h" +#include "vtkInEdgeIterator.h" +#include "vtkInformationDataObjectKey.h" +#include "vtkInformationDataObjectMetaDataKey.h" +#include "vtkInformationDoubleVectorKey.h" +#include "vtkExecutive.h" +#include "vtkVariantArray.h" +#include "vtkStringArray.h" +#include "vtkDoubleArray.h" +#include "vtkFloatArray.h" +#include "vtkCharArray.h" +#include "vtkLongArray.h" +#include "vtkUnsignedCharArray.h" +#include "vtkDataSetAttributes.h" +#include "vtkDemandDrivenPipeline.h" +#include "vtkDataObjectTreeIterator.h" +#include "vtkWarpScalar.h" +#include "vtkQuadratureSchemeDefinition.h" +#include "vtkInformationQuadratureSchemeDefinitionVectorKey.h" +#include "vtkCompositeDataToUnstructuredGridFilter.h" +#include "vtkMultiBlockDataGroupFilter.h" + +#include "MEDCouplingMemArray.hxx" +#include "MEDCouplingUMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" +#include "InterpKernelAutoPtr.hxx" +#include "InterpKernelGaussCoords.hxx" + +#include +#include +#include +#include + +using MEDCoupling::DataArray; +using MEDCoupling::DataArrayInt; +using MEDCoupling::DataArrayDouble; +using MEDCoupling::MEDCouplingMesh; +using MEDCoupling::MEDCouplingUMesh; +using MEDCoupling::DynamicCastSafe; +using MEDCoupling::MEDCouplingFieldDouble; +using MEDCoupling::ON_GAUSS_PT; +using MEDCoupling::MCAuto; + +vtkStandardNewMacro(vtkGaussToCell); + +vtkInformationDoubleVectorKey *GetMEDReaderMetaDataIfAny() +{ + static const char ZE_KEY[]="vtkMEDReader::GAUSS_DATA"; + MEDCoupling::GlobalDict *gd(MEDCoupling::GlobalDict::GetInstance()); + if(!gd->hasKey(ZE_KEY)) + return 0; + std::string ptSt(gd->value(ZE_KEY)); + void *pt(0); + std::istringstream iss(ptSt); iss >> pt; + return reinterpret_cast(pt); +} + +bool IsInformationOK(vtkInformation *info, std::vector& data) +{ + vtkInformationDoubleVectorKey *key(GetMEDReaderMetaDataIfAny()); + if(!key) + return false; + // Check the information contain meta data key + if(!info->Has(key)) + return false; + int lgth(key->Length(info)); + const double *data2(info->Get(key)); + data.insert(data.end(),data2,data2+lgth); + return true; +} + +void ExtractInfo(vtkInformationVector *inputVector, vtkUnstructuredGrid *& usgIn) +{ + vtkInformation *inputInfo(inputVector->GetInformationObject(0)); + vtkDataSet *input(0); + vtkDataSet *input0(vtkDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT()))); + vtkMultiBlockDataSet *input1(vtkMultiBlockDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT()))); + if(input0) + input=input0; + else + { + if(!input1) + throw INTERP_KERNEL::Exception("Input dataSet must be a DataSet or single elt multi block dataset expected !"); + if(input1->GetNumberOfBlocks()!=1) + throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with not exactly one block ! Use MergeBlocks or ExtractBlocks filter before calling this filter !"); + vtkDataObject *input2(input1->GetBlock(0)); + if(!input2) + throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with exactly one block but this single element is NULL !"); + vtkDataSet *input2c(vtkDataSet::SafeDownCast(input2)); + if(!input2c) + throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with exactly one block but this single element is not a dataset ! Use MergeBlocks or ExtractBlocks filter before calling this filter !"); + input=input2c; + } + if(!input) + throw INTERP_KERNEL::Exception("Input data set is NULL !"); + usgIn=vtkUnstructuredGrid::SafeDownCast(input); + if(!usgIn) + throw INTERP_KERNEL::Exception("Input data set is not an unstructured mesh ! This filter works only on unstructured meshes !"); +} + +vtkGaussToCell::vtkGaussToCell():avgStatus(true),maxStatus(false),minStatus(false) +{ + this->SetNumberOfInputPorts(1); + this->SetNumberOfOutputPorts(1); +} + +vtkGaussToCell::~vtkGaussToCell() +{ +} + +void vtkGaussToCell::SetAvgFlag(bool avgStatus) +{ + if(this->avgStatus!=avgStatus) + { + this->avgStatus=avgStatus; + this->Modified(); + } +} + +void vtkGaussToCell::SetMaxFlag(bool maxStatus) +{ + if(this->maxStatus!=maxStatus) + { + this->maxStatus=maxStatus; + this->Modified(); + } +} + +void vtkGaussToCell::SetMinFlag(bool minStatus) +{ + if(this->minStatus!=minStatus) + { + this->minStatus=minStatus; + this->Modified(); + } +} + +int vtkGaussToCell::RequestInformation(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) +{ + //std::cerr << "########################################## vtkGaussToCell::RequestInformation ##########################################" << std::endl; + try + { + vtkUnstructuredGrid *usgIn(0); + ExtractInfo(inputVector[0],usgIn); + } + catch(INTERP_KERNEL::Exception& e) + { + std::ostringstream oss; + oss << "Exception has been thrown in vtkGaussToCell::RequestInformation : " << e.what() << std::endl; + if(this->HasObserver("ErrorEvent") ) + this->InvokeEvent("ErrorEvent",const_cast(oss.str().c_str())); + else + vtkOutputWindowDisplayErrorText(const_cast(oss.str().c_str())); + vtkObject::BreakOnError(); + return 0; + } + return 1; +} + +typedef void (*DataComputer) (const double *inData, const vtkIdType *offData, const std::vector *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData); + +void ComputeAvg(const double *inData, const vtkIdType *offData, const std::vector *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData) +{ + std::fill(outData,outData+outNbCells*zeNbCompo,0.); + for(auto i=0;i *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData) +{ + std::fill(outData,outData+outNbCells*zeNbCompo,-std::numeric_limits::max()); + for(auto i=0;i *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData) +{ + std::fill(outData,outData+outNbCells*zeNbCompo,std::numeric_limits::max()); + for(auto i=0;i *nbgPerCell, vtkIdType outNbCells, vtkDoubleArray *zeOutArray, DataComputer dc) +{ + std::ostringstream oss; + oss << zearray->GetName() << '_' << postName; + int zeNbCompo(zearray->GetNumberOfComponents()); + { + std::string st(oss.str()); + zeOutArray->SetName(st.c_str()); + } + zeOutArray->SetNumberOfComponents(zeNbCompo); + zeOutArray->SetNumberOfTuples(outNbCells); + double *outData(zeOutArray->GetPointer(0)); + const double *inData(zearray->GetPointer(0)); + const auto *offData(offsets->GetPointer(0)); + for(auto i=0;iGetComponentName(i)); + if(comp) + zeOutArray->SetComponentName(i,comp); + } + dc(inData,offData,nbgPerCell,zeNbCompo,outNbCells,outData); +} + +int vtkGaussToCell::RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) +{ + //std::cerr << "########################################## vtkGaussToCell::RequestData ##########################################" << std::endl; + try + { + std::vector GaussAdvData; + bool isOK(IsInformationOK(inputVector[0]->GetInformationObject(0),GaussAdvData)); + if(!isOK) + throw INTERP_KERNEL::Exception("Sorry but no advanced gauss info found ! Expect to be called right after a MEDReader containing Gauss Points !"); + vtkUnstructuredGrid *usgIn(0); + ExtractInfo(inputVector[0],usgIn); + vtkInformation *outInfo(outputVector->GetInformationObject(0)); + vtkUnstructuredGrid *output(vtkUnstructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()))); + output->ShallowCopy(usgIn); + // + std::string zeArrOffset; + int nArrays(usgIn->GetFieldData()->GetNumberOfArrays()); + std::map > offsetKeyMap;//Map storing for each offsets array the corresponding nb of Gauss Points per cell + for(int i=0;iGetFieldData()->GetArray(i)); + if(!array) + continue; + const char* arrayOffsetName(array->GetInformation()->Get(vtkQuadratureSchemeDefinition::QUADRATURE_OFFSET_ARRAY_NAME())); + if(!arrayOffsetName) + continue; + std::string arrOffsetNameCpp(arrayOffsetName); + if(arrOffsetNameCpp.find("ELGA@")==std::string::npos) + continue; + if(zeArrOffset.empty()) + zeArrOffset=arrOffsetNameCpp; + else + if(zeArrOffset!=arrOffsetNameCpp) + { + throw INTERP_KERNEL::Exception("ComputeGaussToCell : error in QUADRATURE_OFFSET_ARRAY_NAME for Gauss fields array !"); + } + vtkDataArray *offTmp(usgIn->GetCellData()->GetArray(zeArrOffset.c_str())); + if(!offTmp) + { + std::ostringstream oss; oss << "ComputeGaussToCell : cell field " << zeArrOffset << " not found !"; + throw INTERP_KERNEL::Exception(oss.str()); + } + vtkIdTypeArray *offsets(vtkIdTypeArray::SafeDownCast(offTmp)); + if(!offsets) + { + std::ostringstream oss; oss << "ComputeGaussToCell : cell field " << zeArrOffset << " exists but not with the right type of data !"; + throw INTERP_KERNEL::Exception(oss.str()); + } + vtkDoubleArray *zearray(vtkDoubleArray::SafeDownCast(array)); + if(!zearray) + continue ; + // + std::map >::iterator nbgPerCellPt(offsetKeyMap.find(offsets)); + const std::vector *nbgPerCell(nullptr); + if(nbgPerCellPt==offsetKeyMap.end()) + { + // fini la parlote + vtkInformation *info(offsets->GetInformation()); + if(!info) + throw INTERP_KERNEL::Exception("info is null ! Internal error ! Looks bad !"); + vtkInformationQuadratureSchemeDefinitionVectorKey *key(vtkQuadratureSchemeDefinition::DICTIONARY()); + if(!key->Has(info)) + throw INTERP_KERNEL::Exception("No quadrature key in info included in offets array ! Internal error ! Looks bad !"); + int dictSize(key->Size(info)); + INTERP_KERNEL::AutoPtr dict(new vtkQuadratureSchemeDefinition *[dictSize]); + key->GetRange(info,dict,0,0,dictSize); + auto nbOfCells(output->GetNumberOfCells()); + std::vector nbg(nbOfCells); + for(auto cellId=0;cellIdGetCellType(cellId)); + vtkQuadratureSchemeDefinition *gaussLoc(dict[ct]); + if(!gaussLoc) + { + std::ostringstream oss; oss << "For cell " << cellId << " no Gauss info attached !"; + throw INTERP_KERNEL::Exception(oss.str()); + } + int np(gaussLoc->GetNumberOfQuadraturePoints()); + nbg[cellId]=np; + } + nbgPerCell=&((*(offsetKeyMap.emplace(offsets,std::move(nbg)).first)).second); + } + else + { + nbgPerCell=&((*nbgPerCellPt).second); + } + auto outNbCells(nbgPerCell->size()); + if(this->avgStatus) + { + vtkNew zeOutArray; + DealWith("avg",zearray,offsets,nbgPerCell,outNbCells,zeOutArray,ComputeAvg); + output->GetCellData()->AddArray(zeOutArray); + } + if(this->maxStatus) + { + vtkNew zeOutArray; + DealWith("max",zearray,offsets,nbgPerCell,outNbCells,zeOutArray,ComputeMax); + output->GetCellData()->AddArray(zeOutArray); + } + if(this->minStatus) + { + vtkNew zeOutArray; + DealWith("min",zearray,offsets,nbgPerCell,outNbCells,zeOutArray,ComputeMin); + output->GetCellData()->AddArray(zeOutArray); + } + } + } + catch(INTERP_KERNEL::Exception& e) + { + std::ostringstream oss; + oss << "Exception has been thrown in vtkGaussToCell::RequestData : " << e.what() << std::endl; + if(this->HasObserver("ErrorEvent") ) + this->InvokeEvent("ErrorEvent",const_cast(oss.str().c_str())); + else + vtkOutputWindowDisplayErrorText(const_cast(oss.str().c_str())); + vtkObject::BreakOnError(); + return 0; + } + return 1; +} + +void vtkGaussToCell::PrintSelf(ostream& os, vtkIndent indent) +{ + this->Superclass::PrintSelf(os, indent); +} diff --git a/src/Plugins/GaussToCell/vtkGaussToCell.h b/src/Plugins/GaussToCell/vtkGaussToCell.h new file mode 100644 index 00000000..8ea188f8 --- /dev/null +++ b/src/Plugins/GaussToCell/vtkGaussToCell.h @@ -0,0 +1,62 @@ +// Copyright (C) 2018 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 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 +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay (EDF R&D) + +#ifndef vtkGaussToCell_h__ +#define vtkGaussToCell_h__ + +#include "vtkUnstructuredGridAlgorithm.h" + +class vtkMutableDirectedGraph; + +class VTK_EXPORT vtkGaussToCell : public vtkUnstructuredGridAlgorithm +{ +public: + static vtkGaussToCell* New(); + vtkTypeMacro(vtkGaussToCell, vtkUnstructuredGridAlgorithm) + void PrintSelf(ostream& os, vtkIndent indent); + + void SetAvgFlag(bool avgStatus); + + void SetMaxFlag(bool maxStatus); + + void SetMinFlag(bool minStatus); + +protected: + vtkGaussToCell(); + ~vtkGaussToCell(); + + int RequestInformation(vtkInformation *request, + vtkInformationVector **inputVector, vtkInformationVector *outputVector); + + int RequestData(vtkInformation *request, vtkInformationVector **inputVector, + vtkInformationVector *outputVector); + +private: + vtkGaussToCell(const vtkGaussToCell&); + void operator=(const vtkGaussToCell&); // Not implemented. + private: + //BTX + //ETX + bool avgStatus; + bool maxStatus; + bool minStatus; +}; + +#endif diff --git a/src/Plugins/StaticMesh/CMakeLists.txt b/src/Plugins/StaticMesh/CMakeLists.txt new file mode 100644 index 00000000..69f10876 --- /dev/null +++ b/src/Plugins/StaticMesh/CMakeLists.txt @@ -0,0 +1,58 @@ +# Copyright (C) 2018 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# + +if (NOT ParaView_BINARY_DIR) + cmake_minimum_required(VERSION 3.3) + find_package(ParaView REQUIRED) + include(${PARAVIEW_USE_FILE}) +endif() + +include_directories( + ${VTK_INCLUDE_DIRS} + ${CMAKE_CURRENT_SOURCE_DIR} + ${CMAKE_CURRENT_BINARY_DIR} + ) + +set(Module_SRCS + vtkStaticDataSetSurfaceFilter.cxx + vtkStaticEnSight6BinaryReader.cxx + vtkStaticEnSight6Reader.cxx + vtkStaticEnSightGoldBinaryReader.cxx + vtkStaticEnSightGoldReader.cxx + vtkStaticMeshObjectFactory.cxx + vtkStaticPlaneCutter.cxx +) +if (PARAVIEW_USE_MPI) + list(APPEND Module_SRCS vtkStaticPUnstructuredGridGhostCellsGenerator.cxx) +endif() + +add_paraview_plugin( + StaticMesh "0.1" + SERVER_MANAGER_XML + StaticMeshSM.xml + SERVER_MANAGER_SOURCES + vtkTemporalUGWavelet.cxx + SOURCES + ${Module_SRCS} + ) + +if (BUILD_TESTING) + add_subdirectory(Testing) +endif() +INSTALL(TARGETS StaticMesh RUNTIME DESTINATION lib/paraview LIBRARY DESTINATION lib/paraview ARCHIVE DESTINATION lib/paraview) diff --git a/src/Plugins/StaticMesh/StaticMeshSM.xml b/src/Plugins/StaticMesh/StaticMeshSM.xml new file mode 100644 index 00000000..bf089f3b --- /dev/null +++ b/src/Plugins/StaticMesh/StaticMeshSM.xml @@ -0,0 +1,154 @@ + + + + + + See Wavelet documentation. + This is a wavelet converted into an unstructured grid. + It contains timesteps yet the mesh is statis over time. + It also contains tPoint and tCell wich are scalar array varying over time. + + + + + + This property specifies the number of time steps + + + + + The six values in this property indicate the X, Y, and Z + extent of the output data. The first two values represent the minimum + and maximum X indices, the next two are the minimum and maximum Y + indices, and the last two are the minimum and maximum Z + indices. + + + + This property specifies the 3D coordinates of the center of + the periodic function (x, y and z in the equation). + + + + This parameter specifies the maximum value (M) of the + function. + + + + This property specifies the natural frequency in X (XF + in the equation). + + + + This property specifies the natural frequency in Y (YF + in the equation). + + + + This property specifies the natural frequency in Z (ZF + in the equation). + + + + This property specifies the wave amplitude in X (XM in + the equation). + + + + This property specifies the wave amplitude in Y (YM in + the equation). + + + + This property specifies the wave amplitude in Z (ZM in + the equation). + + + + This property specifies the standard deviation of the + Gaussian used in computing this function. + + + + This property specifies the rate at which to subsample + the volume. The extent of the dataset in each dimension will be divided + by this value. (See the Whole Extent property.) + + + + + + + + + + + + + + + + + + diff --git a/src/Plugins/StaticMesh/Testing/CMakeLists.txt b/src/Plugins/StaticMesh/Testing/CMakeLists.txt new file mode 100644 index 00000000..936f7ae1 --- /dev/null +++ b/src/Plugins/StaticMesh/Testing/CMakeLists.txt @@ -0,0 +1,39 @@ +include(ParaViewTestingMacros) + +if (PARAVIEW_BUILD_QT_GUI) + add_client_tests("pv" + LOAD_PLUGIN "StaticMesh" + BASELINE_DIR ${PARAVIEW_TEST_BASELINE_DIR} + TEST_SCRIPTS ${CMAKE_CURRENT_SOURCE_DIR}/StaticDataSetSurfaceFilter.xml) + add_client_tests("pv" + LOAD_PLUGIN "StaticMesh" + BASELINE_DIR ${PARAVIEW_TEST_BASELINE_DIR} + TEST_SCRIPTS ${CMAKE_CURRENT_SOURCE_DIR}/StaticPlaneCutter.xml) + add_client_tests("pvcs" + LOAD_PLUGIN "StaticMesh" + BASELINE_DIR ${PARAVIEW_TEST_BASELINE_DIR} + TEST_SCRIPTS ${CMAKE_CURRENT_SOURCE_DIR}/StaticDataSetSurfaceFilter.xml) + add_client_tests("pvcs" + LOAD_PLUGIN "StaticMesh" + BASELINE_DIR ${PARAVIEW_TEST_BASELINE_DIR} + TEST_SCRIPTS ${CMAKE_CURRENT_SOURCE_DIR}/StaticPlaneCutter.xml) + add_client_tests("pvcrs" + LOAD_PLUGIN "StaticMesh" + BASELINE_DIR ${PARAVIEW_TEST_BASELINE_DIR} + TEST_SCRIPTS ${CMAKE_CURRENT_SOURCE_DIR}/StaticDataSetSurfaceFilter.xml) + add_client_tests("pvcrs" + LOAD_PLUGIN "StaticMesh" + BASELINE_DIR ${PARAVIEW_TEST_BASELINE_DIR} + TEST_SCRIPTS ${CMAKE_CURRENT_SOURCE_DIR}/StaticPlaneCutter.xml) + + if (PARAVIEW_USE_MPI) + add_client_server_tests("pvcs" + LOAD_PLUGIN "StaticMesh" + BASELINE_DIR ${PARAVIEW_TEST_BASELINE_DIR} + TEST_SCRIPTS ${CMAKE_CURRENT_SOURCE_DIR}/StaticGhostCellGenerator.xml) + add_client_render_server_tests("pvcrs" + LOAD_PLUGIN "StaticMesh" + BASELINE_DIR ${PARAVIEW_TEST_BASELINE_DIR} + TEST_SCRIPTS ${CMAKE_CURRENT_SOURCE_DIR}/StaticGhostCellGenerator.xml) + endif() +endif() diff --git a/src/Plugins/StaticMesh/Testing/StaticDataSetSurfaceFilter.xml b/src/Plugins/StaticMesh/Testing/StaticDataSetSurfaceFilter.xml new file mode 100644 index 00000000..16baa350 --- /dev/null +++ b/src/Plugins/StaticMesh/Testing/StaticDataSetSurfaceFilter.xml @@ -0,0 +1,21 @@ + + + + + + + + + + + + + + + + + + + + + diff --git a/src/Plugins/StaticMesh/Testing/StaticGhostCellGenerator.xml b/src/Plugins/StaticMesh/Testing/StaticGhostCellGenerator.xml new file mode 100644 index 00000000..655ea949 --- /dev/null +++ b/src/Plugins/StaticMesh/Testing/StaticGhostCellGenerator.xml @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/Plugins/StaticMesh/Testing/StaticPlaneCutter.xml b/src/Plugins/StaticMesh/Testing/StaticPlaneCutter.xml new file mode 100644 index 00000000..37198c3f --- /dev/null +++ b/src/Plugins/StaticMesh/Testing/StaticPlaneCutter.xml @@ -0,0 +1,24 @@ + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/Plugins/StaticMesh/plugin.cmake b/src/Plugins/StaticMesh/plugin.cmake new file mode 100644 index 00000000..fb8142ac --- /dev/null +++ b/src/Plugins/StaticMesh/plugin.cmake @@ -0,0 +1,3 @@ +pv_plugin(StaticMesh + DESCRIPTION "Extend unstructured dataset and selected filters to take benefit of static meshes in transient data" + DEFAULT_ENABLED) diff --git a/src/Plugins/StaticMesh/vtkStaticDataSetSurfaceFilter.cxx b/src/Plugins/StaticMesh/vtkStaticDataSetSurfaceFilter.cxx new file mode 100644 index 00000000..065447e3 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticDataSetSurfaceFilter.cxx @@ -0,0 +1,172 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticDataSetSurfaceFilter.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +#include "vtkStaticDataSetSurfaceFilter.h" + +#include +#include +#include +#include +#include +#include + +vtkStandardNewMacro(vtkStaticDataSetSurfaceFilter); + +//---------------------------------------------------------------------------- +vtkStaticDataSetSurfaceFilter::vtkStaticDataSetSurfaceFilter() +{ + this->InputMeshTime = 0; + this->FilterMTime = 0; +} + +//---------------------------------------------------------------------------- +vtkStaticDataSetSurfaceFilter::~vtkStaticDataSetSurfaceFilter() +{ +} + +//----------------------------------------------------------------------------- +int vtkStaticDataSetSurfaceFilter::UnstructuredGridExecute(vtkDataSet* input, vtkPolyData* output) +{ + vtkUnstructuredGrid* inputUG = vtkUnstructuredGrid::SafeDownCast(input); + if (!inputUG) + { + // Rely on superclass for any input which is not a vtkUnstructuredGrid + return this->Superclass::UnstructuredGridExecute(input, output); + } + + // Check is cache is still valid + if (this->InputMeshTime == inputUG->GetMeshMTime() && this->FilterMTime == this->GetMTime()) + { + // Use cache as base + output->ShallowCopy(this->Cache.Get()); + + // Recover original ids + vtkPointData* outPD = output->GetPointData(); + vtkCellData* outCD = output->GetCellData(); + vtkIdTypeArray* origPointArray = + vtkIdTypeArray::SafeDownCast(outPD->GetArray(this->GetOriginalPointIdsName())); + vtkIdTypeArray* origCellArray = + vtkIdTypeArray::SafeDownCast(outCD->GetArray(this->GetOriginalCellIdsName())); + if (!origPointArray || !origCellArray) + { + vtkErrorMacro( + "OriginalPointIds or OriginalCellIds are missing, cannot use static mesh cache"); + return this->Superclass::UnstructuredGridExecute(input, output); + } + + // Recover input point and cell data + vtkPointData* inPD = input->GetPointData(); + vtkCellData* inCD = input->GetCellData(); + + // Update output point data + vtkIdType* tmpIds = new vtkIdType[origPointArray->GetNumberOfTuples()]; + memcpy(tmpIds, reinterpret_cast(origPointArray->GetVoidPointer(0)), + sizeof(vtkIdType) * origPointArray->GetNumberOfTuples()); + vtkNew pointIds; + pointIds->SetArray(tmpIds, origPointArray->GetNumberOfTuples()); + + // Remove array that have disappeared from input + for (int iArr = outPD->GetNumberOfArrays() - 1; iArr >= 0; iArr--) + { + vtkAbstractArray* inArr = inPD->GetAbstractArray(outPD->GetArrayName(iArr)); + if (!inArr) + { + outPD->RemoveArray(iArr); + } + } + + // Update or create arrays present in input + for (int iArr = 0; iArr < inPD->GetNumberOfArrays(); iArr++) + { + vtkAbstractArray* outArr = outPD->GetAbstractArray(inPD->GetArrayName(iArr)); + if (outArr) + { + inPD->GetAbstractArray(iArr)->GetTuples(pointIds.Get(), outArr); + } + else + { + // New array in input, create it in output + vtkAbstractArray* inArr = inPD->GetAbstractArray(iArr); + outArr = inArr->NewInstance(); + outArr->SetName(inArr->GetName()); + outArr->SetNumberOfTuples(output->GetNumberOfPoints()); + outArr->SetNumberOfComponents(inArr->GetNumberOfComponents()); + inArr->GetTuples(pointIds.Get(), outArr); + outPD->AddArray(outArr); + } + } + + // Update output cell data + tmpIds = new vtkIdType[origCellArray->GetNumberOfTuples()]; + memcpy(tmpIds, reinterpret_cast(origCellArray->GetVoidPointer(0)), + sizeof(vtkIdType) * origCellArray->GetNumberOfTuples()); + vtkNew cellIds; + cellIds->SetArray(tmpIds, origCellArray->GetNumberOfTuples()); + + // Remove array that have disappeared from input + for (int iArr = outCD->GetNumberOfArrays() - 1; iArr >= 0; iArr--) + { + vtkAbstractArray* inArr = inCD->GetAbstractArray(outCD->GetArrayName(iArr)); + if (!inArr) + { + outCD->RemoveArray(iArr); + } + } + + for (int iArr = 0; iArr < inCD->GetNumberOfArrays(); iArr++) + { + vtkAbstractArray* outArr = outCD->GetAbstractArray(inCD->GetArrayName(iArr)); + if (outArr) + { + inCD->GetAbstractArray(iArr)->GetTuples(cellIds.Get(), outArr); + } + else + { + // New array in input, create it in output + vtkAbstractArray* inArr = inCD->GetAbstractArray(iArr); + outArr = inArr->NewInstance(); + outArr->SetName(inArr->GetName()); + outArr->SetNumberOfTuples(output->GetNumberOfCells()); + outArr->SetNumberOfComponents(inArr->GetNumberOfComponents()); + inArr->GetTuples(cellIds.Get(), outArr); + outCD->AddArray(outArr); + } + + } + + // Update output field data + output->GetFieldData()->ShallowCopy(input->GetFieldData()); + return 1; + } + else + { + // Cache is not valid, Execute supercall algorithm + int ret = this->Superclass::UnstructuredGridExecute(input, output); + + // Update the cache with superclass output + this->Cache->ShallowCopy(output); + this->InputMeshTime = inputUG->GetMeshMTime(); + this->FilterMTime = this->GetMTime(); + return ret; + } +} + +//---------------------------------------------------------------------------- +void vtkStaticDataSetSurfaceFilter::PrintSelf(ostream& os, vtkIndent indent) +{ + this->Superclass::PrintSelf(os, indent); + os << indent << "Cache: " << this->Cache << endl; + os << indent << "Input Mesh Time: " << this->InputMeshTime << endl; + os << indent << "Filter mTime: " << this->FilterMTime << endl; +} diff --git a/src/Plugins/StaticMesh/vtkStaticDataSetSurfaceFilter.h b/src/Plugins/StaticMesh/vtkStaticDataSetSurfaceFilter.h new file mode 100644 index 00000000..4158f246 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticDataSetSurfaceFilter.h @@ -0,0 +1,62 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticDataSetSurfaceFilter.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/** + * @class vtkStaticDataSetSurfaceFilter + * @brief Extract the surface of a dataset, optimized for static unstructured grid + * + * vtkStaticDataSetSurfaceFilter is a specialization of vtkDataSetSurfaceFilter + * that uses a cache to store the surface output and reuses it when associated data + * changes over the time, but the geometry of a unstructured grid is static. + * It is to be noted that, since ParaView use the same surface filter + * for each block of a MultiBlock, this filter is not effective with multiblock + * dataset. + * + * @sa + * vtkStaticMeshObjectFactory +*/ + +#ifndef vtkStaticDataSetSurfaceFilter_h +#define vtkStaticDataSetSurfaceFilter_h + +#include "vtkDataSetSurfaceFilter.h" +#include "vtkNew.h" + +class vtkPolyData; + +class vtkStaticDataSetSurfaceFilter : public vtkDataSetSurfaceFilter +{ +public: + static vtkStaticDataSetSurfaceFilter* New(); + typedef vtkDataSetSurfaceFilter + Superclass; // vtkTypeMacro can't be used with a factory built object + void PrintSelf(ostream& os, vtkIndent indent) VTK_OVERRIDE; + + int UnstructuredGridExecute(vtkDataSet* input, vtkPolyData* output) VTK_OVERRIDE; + +protected: + vtkStaticDataSetSurfaceFilter(); + ~vtkStaticDataSetSurfaceFilter() VTK_OVERRIDE; + + vtkNew Cache; + vtkMTimeType InputMeshTime; + vtkMTimeType FilterMTime; + +private: + // Hide these from the user and the compiler. + vtkStaticDataSetSurfaceFilter(const vtkStaticDataSetSurfaceFilter&) VTK_DELETE_FUNCTION; + void operator=(const vtkStaticDataSetSurfaceFilter&) VTK_DELETE_FUNCTION; +}; + +#endif diff --git a/src/Plugins/StaticMesh/vtkStaticEnSight6BinaryReader.cxx b/src/Plugins/StaticMesh/vtkStaticEnSight6BinaryReader.cxx new file mode 100644 index 00000000..aef35889 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticEnSight6BinaryReader.cxx @@ -0,0 +1,275 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticEnSight6BinaryReader.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +#include "vtkStaticEnSight6BinaryReader.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +vtkStandardNewMacro(vtkStaticEnSight6BinaryReader); + +//---------------------------------------------------------------------------- +int vtkStaticEnSight6BinaryReader::RequestData( + vtkInformation *vtkNotUsed(request), + vtkInformationVector **vtkNotUsed(inputVector), + vtkInformationVector *outputVector) +{ + vtkDebugMacro("In execute "); + + vtkInformation *outInfo = outputVector->GetInformationObject(0); + vtkMultiBlockDataSet *output = vtkMultiBlockDataSet::SafeDownCast( + outInfo->Get(vtkDataObject::DATA_OBJECT())); + + int tsLength = + outInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS()); + double* steps = + outInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS()); + + this->ActualTimeValue = this->TimeValue; + + // Check if a particular time was requested by the pipeline. + // This overrides the ivar. + if(outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()) && tsLength>0) + { + // Get the requested time step. We only support requests of a single time + // step in this reader right now + double requestedTimeStep = + outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()); + + // find the first time value larger than requested time value + // this logic could be improved + int cnt = 0; + while (cnt < tsLength-1 && steps[cnt] < requestedTimeStep) + { + cnt++; + } + this->ActualTimeValue = steps[cnt]; + } + + vtkDebugMacro("Executing with: " << this->ActualTimeValue); + + if (this->CacheMTime < this->GetMTime()) + { + int i, timeSet, fileSet, timeStep, timeStepInFile, fileNum; + vtkDataArray *times; + vtkIdList *numStepsList, *filenameNumbers; + float newTime; + int numSteps; + char* fileName; + int filenameNum; + + if ( ! this->CaseFileRead) + { + vtkErrorMacro("error reading case file"); + return 0; + } + + this->NumberOfNewOutputs = 0; + this->NumberOfGeometryParts = 0; + if (this->GeometryFileName) + { + timeStep = timeStepInFile = 1; + fileNum = 1; + fileName = new char[strlen(this->GeometryFileName) + 10]; + strcpy(fileName, this->GeometryFileName); + + if (this->UseTimeSets) + { + timeSet = this->TimeSetIds->IsId(this->GeometryTimeSet); + if (timeSet >= 0) + { + times = this->TimeSets->GetItem(timeSet); + this->GeometryTimeValue = times->GetComponent(0, 0); + for (i = 1; i < times->GetNumberOfTuples(); i++) + { + newTime = times->GetComponent(i, 0); + if (newTime <= this->ActualTimeValue && + newTime > this->GeometryTimeValue) + { + this->GeometryTimeValue = newTime; + timeStep++; + timeStepInFile++; + } + } + if (this->TimeSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->TimeSetsWithFilenameNumbers-> + IsId(this->GeometryTimeSet); + if (collectionNum > -1) + { + filenameNumbers = + this->TimeSetFileNameNumbers->GetItem(collectionNum); + filenameNum = filenameNumbers->GetId(timeStep-1); + if (! this->UseFileSets) + { + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + + // There can only be file sets if there are also time sets. + if (this->UseFileSets) + { + fileSet = this->FileSets->IsId(this->GeometryFileSet); + numStepsList = static_cast(this->FileSetNumberOfSteps-> + GetItemAsObject(fileSet)); + + if (timeStep > numStepsList->GetId(0)) + { + numSteps = numStepsList->GetId(0); + timeStepInFile -= numSteps; + fileNum = 2; + for (i = 1; i < numStepsList->GetNumberOfIds(); i++) + { + numSteps += numStepsList->GetId(i); + if (timeStep > numSteps) + { + fileNum++; + timeStepInFile -= numStepsList->GetId(i); + } + } + } + if (this->FileSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->FileSetsWithFilenameNumbers-> + IsId(this->GeometryFileSet); + if (collectionNum > -1) + { + filenameNumbers = this->FileSetFileNameNumbers-> + GetItem(collectionNum); + filenameNum = filenameNumbers->GetId(fileNum-1); + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + } + } + + if (!this->ReadGeometryFile(fileName, timeStepInFile, this->Cache)) + { + vtkErrorMacro("error reading geometry file"); + delete [] fileName; + return 0; + } + + delete [] fileName; + } + if (this->MeasuredFileName) + { + timeStep = timeStepInFile = 1; + fileNum = 1; + fileName = new char[strlen(this->MeasuredFileName) + 10]; + strcpy(fileName, this->MeasuredFileName); + + if (this->UseTimeSets) + { + timeSet = this->TimeSetIds->IsId(this->MeasuredTimeSet); + if (timeSet >= 0) + { + times = this->TimeSets->GetItem(timeSet); + this->MeasuredTimeValue = times->GetComponent(0, 0); + for (i = 1; i < times->GetNumberOfTuples(); i++) + { + newTime = times->GetComponent(i, 0); + if (newTime <= this->ActualTimeValue && + newTime > this->MeasuredTimeValue) + { + this->MeasuredTimeValue = newTime; + timeStep++; + timeStepInFile++; + } + } + if (this->TimeSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->TimeSetsWithFilenameNumbers-> + IsId(this->MeasuredTimeSet); + if (collectionNum > -1) + { + filenameNumbers = this->TimeSetFileNameNumbers-> + GetItem(collectionNum); + filenameNum = filenameNumbers->GetId(timeStep-1); + if (! this->UseFileSets) + { + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + + // There can only be file sets if there are also time sets. + if (this->UseFileSets) + { + fileSet = this->FileSets->IsId(this->MeasuredFileSet); + numStepsList = static_cast(this->FileSetNumberOfSteps-> + GetItemAsObject(fileSet)); + + if (timeStep > numStepsList->GetId(0)) + { + numSteps = numStepsList->GetId(0); + timeStepInFile -= numSteps; + fileNum = 2; + for (i = 1; i < numStepsList->GetNumberOfIds(); i++) + { + numSteps += numStepsList->GetId(i); + if (timeStep > numSteps) + { + fileNum++; + timeStepInFile -= numStepsList->GetId(i); + } + } + } + if (this->FileSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->FileSetsWithFilenameNumbers-> + IsId(this->MeasuredFileSet); + if (collectionNum > -1) + { + filenameNumbers = this->FileSetFileNameNumbers-> + GetItem(fileSet); + filenameNum = filenameNumbers->GetId(fileNum-1); + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + } + } + if (!this->ReadMeasuredGeometryFile(fileName, timeStepInFile, this->Cache)) + { + vtkErrorMacro("error reading measured geometry file"); + delete [] fileName; + return 0; + } + delete [] fileName; + } + this->CacheMTime.Modified(); + } + output->ShallowCopy(this->Cache); + + if ((this->NumberOfVariables + this->NumberOfComplexVariables) > 0) + { + if (!this->ReadVariableFiles(output)) + { + vtkErrorMacro("error reading variable files"); + return 0; + } + } + + return 1; +} diff --git a/src/Plugins/StaticMesh/vtkStaticEnSight6BinaryReader.h b/src/Plugins/StaticMesh/vtkStaticEnSight6BinaryReader.h new file mode 100644 index 00000000..58680aa8 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticEnSight6BinaryReader.h @@ -0,0 +1,69 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticEnSight6BinaryReader.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/** + * @class vtkStaticEnSight6BinaryReader + * @brief class to read binary EnSight6 files + * + * vtkStaticEnSight6BinaryReader is a class to read binary EnSight6 files into vtk. + * Because the different parts of the EnSight data can be of various data + * types, this reader produces multiple outputs, one per part in the input + * file. + * All variable information is being stored in field data. The descriptions + * listed in the case file are used as the array names in the field data. + * For complex vector variables, the description is appended with _r (for the + * array of real values) and _i (for the array if imaginary values). Complex + * scalar variables are stored as a single array with 2 components, real and + * imaginary, listed in that order. + * @warning + * You must manually call Update on this reader and then connect the rest + * of the pipeline because (due to the nature of the file format) it is + * not possible to know ahead of time how many outputs you will have or + * what types they will be. + * This reader can only handle static EnSight datasets (both static geometry + * and variables). +*/ + +#ifndef vtkStaticEnSight6BinaryReader_h +#define vtkStaticEnSight6BinaryReader_h + +#include "vtkEnSight6BinaryReader.h" +#include "vtkNew.h" + +class vtkMultiBlockDataSet; + +class vtkStaticEnSight6BinaryReader : public vtkEnSight6BinaryReader +{ +public: + static vtkStaticEnSight6BinaryReader *New(); + vtkTypeMacro(vtkStaticEnSight6BinaryReader, vtkEnSight6BinaryReader); + +protected: + vtkStaticEnSight6BinaryReader() = default; + ~vtkStaticEnSight6BinaryReader() override = default; + + int RequestData(vtkInformation*, + vtkInformationVector**, + vtkInformationVector*) override; + + vtkNew Cache; + vtkTimeStamp CacheMTime; + +private: + vtkStaticEnSight6BinaryReader(const vtkStaticEnSight6BinaryReader&) = delete; + void operator=(const vtkStaticEnSight6BinaryReader&) = delete; +}; + +#endif + diff --git a/src/Plugins/StaticMesh/vtkStaticEnSight6Reader.cxx b/src/Plugins/StaticMesh/vtkStaticEnSight6Reader.cxx new file mode 100644 index 00000000..8ad68cd4 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticEnSight6Reader.cxx @@ -0,0 +1,275 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticEnSight6Reader.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +#include "vtkStaticEnSight6Reader.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +vtkStandardNewMacro(vtkStaticEnSight6Reader); + +//---------------------------------------------------------------------------- +int vtkStaticEnSight6Reader::RequestData( + vtkInformation *vtkNotUsed(request), + vtkInformationVector **vtkNotUsed(inputVector), + vtkInformationVector *outputVector) +{ + vtkDebugMacro("In execute "); + + vtkInformation *outInfo = outputVector->GetInformationObject(0); + vtkMultiBlockDataSet *output = vtkMultiBlockDataSet::SafeDownCast( + outInfo->Get(vtkDataObject::DATA_OBJECT())); + + int tsLength = + outInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS()); + double* steps = + outInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS()); + + this->ActualTimeValue = this->TimeValue; + + // Check if a particular time was requested by the pipeline. + // This overrides the ivar. + if(outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()) && tsLength>0) + { + // Get the requested time step. We only support requests of a single time + // step in this reader right now + double requestedTimeStep = + outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()); + + // find the first time value larger than requested time value + // this logic could be improved + int cnt = 0; + while (cnt < tsLength-1 && steps[cnt] < requestedTimeStep) + { + cnt++; + } + this->ActualTimeValue = steps[cnt]; + } + + vtkDebugMacro("Executing with: " << this->ActualTimeValue); + + if (this->CacheMTime < this->GetMTime()) + { + int i, timeSet, fileSet, timeStep, timeStepInFile, fileNum; + vtkDataArray *times; + vtkIdList *numStepsList, *filenameNumbers; + float newTime; + int numSteps; + char* fileName; + int filenameNum; + + if ( ! this->CaseFileRead) + { + vtkErrorMacro("error reading case file"); + return 0; + } + + this->NumberOfNewOutputs = 0; + this->NumberOfGeometryParts = 0; + if (this->GeometryFileName) + { + timeStep = timeStepInFile = 1; + fileNum = 1; + fileName = new char[strlen(this->GeometryFileName) + 10]; + strcpy(fileName, this->GeometryFileName); + + if (this->UseTimeSets) + { + timeSet = this->TimeSetIds->IsId(this->GeometryTimeSet); + if (timeSet >= 0) + { + times = this->TimeSets->GetItem(timeSet); + this->GeometryTimeValue = times->GetComponent(0, 0); + for (i = 1; i < times->GetNumberOfTuples(); i++) + { + newTime = times->GetComponent(i, 0); + if (newTime <= this->ActualTimeValue && + newTime > this->GeometryTimeValue) + { + this->GeometryTimeValue = newTime; + timeStep++; + timeStepInFile++; + } + } + if (this->TimeSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->TimeSetsWithFilenameNumbers-> + IsId(this->GeometryTimeSet); + if (collectionNum > -1) + { + filenameNumbers = + this->TimeSetFileNameNumbers->GetItem(collectionNum); + filenameNum = filenameNumbers->GetId(timeStep-1); + if (! this->UseFileSets) + { + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + + // There can only be file sets if there are also time sets. + if (this->UseFileSets) + { + fileSet = this->FileSets->IsId(this->GeometryFileSet); + numStepsList = static_cast(this->FileSetNumberOfSteps-> + GetItemAsObject(fileSet)); + + if (timeStep > numStepsList->GetId(0)) + { + numSteps = numStepsList->GetId(0); + timeStepInFile -= numSteps; + fileNum = 2; + for (i = 1; i < numStepsList->GetNumberOfIds(); i++) + { + numSteps += numStepsList->GetId(i); + if (timeStep > numSteps) + { + fileNum++; + timeStepInFile -= numStepsList->GetId(i); + } + } + } + if (this->FileSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->FileSetsWithFilenameNumbers-> + IsId(this->GeometryFileSet); + if (collectionNum > -1) + { + filenameNumbers = this->FileSetFileNameNumbers-> + GetItem(collectionNum); + filenameNum = filenameNumbers->GetId(fileNum-1); + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + } + } + + if (!this->ReadGeometryFile(fileName, timeStepInFile, this->Cache)) + { + vtkErrorMacro("error reading geometry file"); + delete [] fileName; + return 0; + } + + delete [] fileName; + } + if (this->MeasuredFileName) + { + timeStep = timeStepInFile = 1; + fileNum = 1; + fileName = new char[strlen(this->MeasuredFileName) + 10]; + strcpy(fileName, this->MeasuredFileName); + + if (this->UseTimeSets) + { + timeSet = this->TimeSetIds->IsId(this->MeasuredTimeSet); + if (timeSet >= 0) + { + times = this->TimeSets->GetItem(timeSet); + this->MeasuredTimeValue = times->GetComponent(0, 0); + for (i = 1; i < times->GetNumberOfTuples(); i++) + { + newTime = times->GetComponent(i, 0); + if (newTime <= this->ActualTimeValue && + newTime > this->MeasuredTimeValue) + { + this->MeasuredTimeValue = newTime; + timeStep++; + timeStepInFile++; + } + } + if (this->TimeSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->TimeSetsWithFilenameNumbers-> + IsId(this->MeasuredTimeSet); + if (collectionNum > -1) + { + filenameNumbers = this->TimeSetFileNameNumbers-> + GetItem(collectionNum); + filenameNum = filenameNumbers->GetId(timeStep-1); + if (! this->UseFileSets) + { + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + + // There can only be file sets if there are also time sets. + if (this->UseFileSets) + { + fileSet = this->FileSets->IsId(this->MeasuredFileSet); + numStepsList = static_cast(this->FileSetNumberOfSteps-> + GetItemAsObject(fileSet)); + + if (timeStep > numStepsList->GetId(0)) + { + numSteps = numStepsList->GetId(0); + timeStepInFile -= numSteps; + fileNum = 2; + for (i = 1; i < numStepsList->GetNumberOfIds(); i++) + { + numSteps += numStepsList->GetId(i); + if (timeStep > numSteps) + { + fileNum++; + timeStepInFile -= numStepsList->GetId(i); + } + } + } + if (this->FileSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->FileSetsWithFilenameNumbers-> + IsId(this->MeasuredFileSet); + if (collectionNum > -1) + { + filenameNumbers = this->FileSetFileNameNumbers-> + GetItem(fileSet); + filenameNum = filenameNumbers->GetId(fileNum-1); + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + } + } + if (!this->ReadMeasuredGeometryFile(fileName, timeStepInFile, this->Cache)) + { + vtkErrorMacro("error reading measured geometry file"); + delete [] fileName; + return 0; + } + delete [] fileName; + } + this->CacheMTime.Modified(); + } + output->ShallowCopy(this->Cache); + + if ((this->NumberOfVariables + this->NumberOfComplexVariables) > 0) + { + if (!this->ReadVariableFiles(output)) + { + vtkErrorMacro("error reading variable files"); + return 0; + } + } + + return 1; +} diff --git a/src/Plugins/StaticMesh/vtkStaticEnSight6Reader.h b/src/Plugins/StaticMesh/vtkStaticEnSight6Reader.h new file mode 100644 index 00000000..d4a2d6a8 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticEnSight6Reader.h @@ -0,0 +1,69 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticEnSight6Reader.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/** + * @class vtkStaticEnSight6Reader + * @brief class to read binary EnSight6 files + * + * vtkStaticEnSight6Reader is a class to read binary EnSight6 files into vtk. + * Because the different parts of the EnSight data can be of various data + * types, this reader produces multiple outputs, one per part in the input + * file. + * All variable information is being stored in field data. The descriptions + * listed in the case file are used as the array names in the field data. + * For complex vector variables, the description is appended with _r (for the + * array of real values) and _i (for the array if imaginary values). Complex + * scalar variables are stored as a single array with 2 components, real and + * imaginary, listed in that order. + * @warning + * You must manually call Update on this reader and then connect the rest + * of the pipeline because (due to the nature of the file format) it is + * not possible to know ahead of time how many outputs you will have or + * what types they will be. + * This reader can only handle static EnSight datasets (both static geometry + * and variables). +*/ + +#ifndef vtkStaticEnSight6Reader_h +#define vtkStaticEnSight6Reader_h + +#include "vtkEnSight6Reader.h" +#include "vtkNew.h" + +class vtkMultiBlockDataSet; + +class vtkStaticEnSight6Reader : public vtkEnSight6Reader +{ +public: + static vtkStaticEnSight6Reader *New(); + vtkTypeMacro(vtkStaticEnSight6Reader, vtkEnSight6Reader); + +protected: + vtkStaticEnSight6Reader() = default; + ~vtkStaticEnSight6Reader() override = default; + + int RequestData(vtkInformation*, + vtkInformationVector**, + vtkInformationVector*) override; + + vtkNew Cache; + vtkTimeStamp CacheMTime; + +private: + vtkStaticEnSight6Reader(const vtkStaticEnSight6Reader&) = delete; + void operator=(const vtkStaticEnSight6Reader&) = delete; +}; + +#endif + diff --git a/src/Plugins/StaticMesh/vtkStaticEnSightGoldBinaryReader.cxx b/src/Plugins/StaticMesh/vtkStaticEnSightGoldBinaryReader.cxx new file mode 100644 index 00000000..6d336de8 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticEnSightGoldBinaryReader.cxx @@ -0,0 +1,275 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticEnSightGoldBinaryReader.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +#include "vtkStaticEnSightGoldBinaryReader.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +vtkStandardNewMacro(vtkStaticEnSightGoldBinaryReader); + +//---------------------------------------------------------------------------- +int vtkStaticEnSightGoldBinaryReader::RequestData( + vtkInformation *vtkNotUsed(request), + vtkInformationVector **vtkNotUsed(inputVector), + vtkInformationVector *outputVector) +{ + vtkDebugMacro("In execute "); + + vtkInformation *outInfo = outputVector->GetInformationObject(0); + vtkMultiBlockDataSet *output = vtkMultiBlockDataSet::SafeDownCast( + outInfo->Get(vtkDataObject::DATA_OBJECT())); + + int tsLength = + outInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS()); + double* steps = + outInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS()); + + this->ActualTimeValue = this->TimeValue; + + // Check if a particular time was requested by the pipeline. + // This overrides the ivar. + if(outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()) && tsLength>0) + { + // Get the requested time step. We only support requests of a single time + // step in this reader right now + double requestedTimeStep = + outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()); + + // find the first time value larger than requested time value + // this logic could be improved + int cnt = 0; + while (cnt < tsLength-1 && steps[cnt] < requestedTimeStep) + { + cnt++; + } + this->ActualTimeValue = steps[cnt]; + } + + vtkDebugMacro("Executing with: " << this->ActualTimeValue); + + if (this->CacheMTime < this->GetMTime()) + { + int i, timeSet, fileSet, timeStep, timeStepInFile, fileNum; + vtkDataArray *times; + vtkIdList *numStepsList, *filenameNumbers; + float newTime; + int numSteps; + char* fileName; + int filenameNum; + + if ( ! this->CaseFileRead) + { + vtkErrorMacro("error reading case file"); + return 0; + } + + this->NumberOfNewOutputs = 0; + this->NumberOfGeometryParts = 0; + if (this->GeometryFileName) + { + timeStep = timeStepInFile = 1; + fileNum = 1; + fileName = new char[strlen(this->GeometryFileName) + 10]; + strcpy(fileName, this->GeometryFileName); + + if (this->UseTimeSets) + { + timeSet = this->TimeSetIds->IsId(this->GeometryTimeSet); + if (timeSet >= 0) + { + times = this->TimeSets->GetItem(timeSet); + this->GeometryTimeValue = times->GetComponent(0, 0); + for (i = 1; i < times->GetNumberOfTuples(); i++) + { + newTime = times->GetComponent(i, 0); + if (newTime <= this->ActualTimeValue && + newTime > this->GeometryTimeValue) + { + this->GeometryTimeValue = newTime; + timeStep++; + timeStepInFile++; + } + } + if (this->TimeSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->TimeSetsWithFilenameNumbers-> + IsId(this->GeometryTimeSet); + if (collectionNum > -1) + { + filenameNumbers = + this->TimeSetFileNameNumbers->GetItem(collectionNum); + filenameNum = filenameNumbers->GetId(timeStep-1); + if (! this->UseFileSets) + { + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + + // There can only be file sets if there are also time sets. + if (this->UseFileSets) + { + fileSet = this->FileSets->IsId(this->GeometryFileSet); + numStepsList = static_cast(this->FileSetNumberOfSteps-> + GetItemAsObject(fileSet)); + + if (timeStep > numStepsList->GetId(0)) + { + numSteps = numStepsList->GetId(0); + timeStepInFile -= numSteps; + fileNum = 2; + for (i = 1; i < numStepsList->GetNumberOfIds(); i++) + { + numSteps += numStepsList->GetId(i); + if (timeStep > numSteps) + { + fileNum++; + timeStepInFile -= numStepsList->GetId(i); + } + } + } + if (this->FileSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->FileSetsWithFilenameNumbers-> + IsId(this->GeometryFileSet); + if (collectionNum > -1) + { + filenameNumbers = this->FileSetFileNameNumbers-> + GetItem(collectionNum); + filenameNum = filenameNumbers->GetId(fileNum-1); + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + } + } + + if (!this->ReadGeometryFile(fileName, timeStepInFile, this->Cache)) + { + vtkErrorMacro("error reading geometry file"); + delete [] fileName; + return 0; + } + + delete [] fileName; + } + if (this->MeasuredFileName) + { + timeStep = timeStepInFile = 1; + fileNum = 1; + fileName = new char[strlen(this->MeasuredFileName) + 10]; + strcpy(fileName, this->MeasuredFileName); + + if (this->UseTimeSets) + { + timeSet = this->TimeSetIds->IsId(this->MeasuredTimeSet); + if (timeSet >= 0) + { + times = this->TimeSets->GetItem(timeSet); + this->MeasuredTimeValue = times->GetComponent(0, 0); + for (i = 1; i < times->GetNumberOfTuples(); i++) + { + newTime = times->GetComponent(i, 0); + if (newTime <= this->ActualTimeValue && + newTime > this->MeasuredTimeValue) + { + this->MeasuredTimeValue = newTime; + timeStep++; + timeStepInFile++; + } + } + if (this->TimeSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->TimeSetsWithFilenameNumbers-> + IsId(this->MeasuredTimeSet); + if (collectionNum > -1) + { + filenameNumbers = this->TimeSetFileNameNumbers-> + GetItem(collectionNum); + filenameNum = filenameNumbers->GetId(timeStep-1); + if (! this->UseFileSets) + { + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + + // There can only be file sets if there are also time sets. + if (this->UseFileSets) + { + fileSet = this->FileSets->IsId(this->MeasuredFileSet); + numStepsList = static_cast(this->FileSetNumberOfSteps-> + GetItemAsObject(fileSet)); + + if (timeStep > numStepsList->GetId(0)) + { + numSteps = numStepsList->GetId(0); + timeStepInFile -= numSteps; + fileNum = 2; + for (i = 1; i < numStepsList->GetNumberOfIds(); i++) + { + numSteps += numStepsList->GetId(i); + if (timeStep > numSteps) + { + fileNum++; + timeStepInFile -= numStepsList->GetId(i); + } + } + } + if (this->FileSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->FileSetsWithFilenameNumbers-> + IsId(this->MeasuredFileSet); + if (collectionNum > -1) + { + filenameNumbers = this->FileSetFileNameNumbers-> + GetItem(fileSet); + filenameNum = filenameNumbers->GetId(fileNum-1); + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + } + } + if (!this->ReadMeasuredGeometryFile(fileName, timeStepInFile, this->Cache)) + { + vtkErrorMacro("error reading measured geometry file"); + delete [] fileName; + return 0; + } + delete [] fileName; + } + this->CacheMTime.Modified(); + } + output->ShallowCopy(this->Cache); + + if ((this->NumberOfVariables + this->NumberOfComplexVariables) > 0) + { + if (!this->ReadVariableFiles(output)) + { + vtkErrorMacro("error reading variable files"); + return 0; + } + } + + return 1; +} diff --git a/src/Plugins/StaticMesh/vtkStaticEnSightGoldBinaryReader.h b/src/Plugins/StaticMesh/vtkStaticEnSightGoldBinaryReader.h new file mode 100644 index 00000000..7a815743 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticEnSightGoldBinaryReader.h @@ -0,0 +1,69 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticEnSightGoldBinaryReader.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/** + * @class vtkStaticEnSightGoldBinaryReader + * @brief class to read binary EnSight6 files + * + * vtkStaticEnSightGoldBinaryReader is a class to read binary EnSight6 files into vtk. + * Because the different parts of the EnSight data can be of various data + * types, this reader produces multiple outputs, one per part in the input + * file. + * All variable information is being stored in field data. The descriptions + * listed in the case file are used as the array names in the field data. + * For complex vector variables, the description is appended with _r (for the + * array of real values) and _i (for the array if imaginary values). Complex + * scalar variables are stored as a single array with 2 components, real and + * imaginary, listed in that order. + * @warning + * You must manually call Update on this reader and then connect the rest + * of the pipeline because (due to the nature of the file format) it is + * not possible to know ahead of time how many outputs you will have or + * what types they will be. + * This reader can only handle static EnSight datasets (both static geometry + * and variables). +*/ + +#ifndef vtkStaticEnSightGoldBinaryReader_h +#define vtkStaticEnSightGoldBinaryReader_h + +#include "vtkEnSightGoldBinaryReader.h" +#include "vtkNew.h" + +class vtkMultiBlockDataSet; + +class vtkStaticEnSightGoldBinaryReader : public vtkEnSightGoldBinaryReader +{ +public: + static vtkStaticEnSightGoldBinaryReader *New(); + vtkTypeMacro(vtkStaticEnSightGoldBinaryReader, vtkEnSightGoldBinaryReader); + +protected: + vtkStaticEnSightGoldBinaryReader() = default; + ~vtkStaticEnSightGoldBinaryReader() override = default; + + int RequestData(vtkInformation*, + vtkInformationVector**, + vtkInformationVector*) override; + + vtkNew Cache; + vtkTimeStamp CacheMTime; + +private: + vtkStaticEnSightGoldBinaryReader(const vtkStaticEnSightGoldBinaryReader&) = delete; + void operator=(const vtkStaticEnSightGoldBinaryReader&) = delete; +}; + +#endif + diff --git a/src/Plugins/StaticMesh/vtkStaticEnSightGoldReader.cxx b/src/Plugins/StaticMesh/vtkStaticEnSightGoldReader.cxx new file mode 100644 index 00000000..11a62154 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticEnSightGoldReader.cxx @@ -0,0 +1,275 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticEnSightGoldReader.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +#include "vtkStaticEnSightGoldReader.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +vtkStandardNewMacro(vtkStaticEnSightGoldReader); + +//---------------------------------------------------------------------------- +int vtkStaticEnSightGoldReader::RequestData( + vtkInformation *vtkNotUsed(request), + vtkInformationVector **vtkNotUsed(inputVector), + vtkInformationVector *outputVector) +{ + vtkDebugMacro("In execute "); + + vtkInformation *outInfo = outputVector->GetInformationObject(0); + vtkMultiBlockDataSet *output = vtkMultiBlockDataSet::SafeDownCast( + outInfo->Get(vtkDataObject::DATA_OBJECT())); + + int tsLength = + outInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS()); + double* steps = + outInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS()); + + this->ActualTimeValue = this->TimeValue; + + // Check if a particular time was requested by the pipeline. + // This overrides the ivar. + if(outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()) && tsLength>0) + { + // Get the requested time step. We only support requests of a single time + // step in this reader right now + double requestedTimeStep = + outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()); + + // find the first time value larger than requested time value + // this logic could be improved + int cnt = 0; + while (cnt < tsLength-1 && steps[cnt] < requestedTimeStep) + { + cnt++; + } + this->ActualTimeValue = steps[cnt]; + } + + vtkDebugMacro("Executing with: " << this->ActualTimeValue); + + if (this->CacheMTime < this->GetMTime()) + { + int i, timeSet, fileSet, timeStep, timeStepInFile, fileNum; + vtkDataArray *times; + vtkIdList *numStepsList, *filenameNumbers; + float newTime; + int numSteps; + char* fileName; + int filenameNum; + + if ( ! this->CaseFileRead) + { + vtkErrorMacro("error reading case file"); + return 0; + } + + this->NumberOfNewOutputs = 0; + this->NumberOfGeometryParts = 0; + if (this->GeometryFileName) + { + timeStep = timeStepInFile = 1; + fileNum = 1; + fileName = new char[strlen(this->GeometryFileName) + 10]; + strcpy(fileName, this->GeometryFileName); + + if (this->UseTimeSets) + { + timeSet = this->TimeSetIds->IsId(this->GeometryTimeSet); + if (timeSet >= 0) + { + times = this->TimeSets->GetItem(timeSet); + this->GeometryTimeValue = times->GetComponent(0, 0); + for (i = 1; i < times->GetNumberOfTuples(); i++) + { + newTime = times->GetComponent(i, 0); + if (newTime <= this->ActualTimeValue && + newTime > this->GeometryTimeValue) + { + this->GeometryTimeValue = newTime; + timeStep++; + timeStepInFile++; + } + } + if (this->TimeSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->TimeSetsWithFilenameNumbers-> + IsId(this->GeometryTimeSet); + if (collectionNum > -1) + { + filenameNumbers = + this->TimeSetFileNameNumbers->GetItem(collectionNum); + filenameNum = filenameNumbers->GetId(timeStep-1); + if (! this->UseFileSets) + { + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + + // There can only be file sets if there are also time sets. + if (this->UseFileSets) + { + fileSet = this->FileSets->IsId(this->GeometryFileSet); + numStepsList = static_cast(this->FileSetNumberOfSteps-> + GetItemAsObject(fileSet)); + + if (timeStep > numStepsList->GetId(0)) + { + numSteps = numStepsList->GetId(0); + timeStepInFile -= numSteps; + fileNum = 2; + for (i = 1; i < numStepsList->GetNumberOfIds(); i++) + { + numSteps += numStepsList->GetId(i); + if (timeStep > numSteps) + { + fileNum++; + timeStepInFile -= numStepsList->GetId(i); + } + } + } + if (this->FileSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->FileSetsWithFilenameNumbers-> + IsId(this->GeometryFileSet); + if (collectionNum > -1) + { + filenameNumbers = this->FileSetFileNameNumbers-> + GetItem(collectionNum); + filenameNum = filenameNumbers->GetId(fileNum-1); + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + } + } + + if (!this->ReadGeometryFile(fileName, timeStepInFile, this->Cache)) + { + vtkErrorMacro("error reading geometry file"); + delete [] fileName; + return 0; + } + + delete [] fileName; + } + if (this->MeasuredFileName) + { + timeStep = timeStepInFile = 1; + fileNum = 1; + fileName = new char[strlen(this->MeasuredFileName) + 10]; + strcpy(fileName, this->MeasuredFileName); + + if (this->UseTimeSets) + { + timeSet = this->TimeSetIds->IsId(this->MeasuredTimeSet); + if (timeSet >= 0) + { + times = this->TimeSets->GetItem(timeSet); + this->MeasuredTimeValue = times->GetComponent(0, 0); + for (i = 1; i < times->GetNumberOfTuples(); i++) + { + newTime = times->GetComponent(i, 0); + if (newTime <= this->ActualTimeValue && + newTime > this->MeasuredTimeValue) + { + this->MeasuredTimeValue = newTime; + timeStep++; + timeStepInFile++; + } + } + if (this->TimeSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->TimeSetsWithFilenameNumbers-> + IsId(this->MeasuredTimeSet); + if (collectionNum > -1) + { + filenameNumbers = this->TimeSetFileNameNumbers-> + GetItem(collectionNum); + filenameNum = filenameNumbers->GetId(timeStep-1); + if (! this->UseFileSets) + { + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + + // There can only be file sets if there are also time sets. + if (this->UseFileSets) + { + fileSet = this->FileSets->IsId(this->MeasuredFileSet); + numStepsList = static_cast(this->FileSetNumberOfSteps-> + GetItemAsObject(fileSet)); + + if (timeStep > numStepsList->GetId(0)) + { + numSteps = numStepsList->GetId(0); + timeStepInFile -= numSteps; + fileNum = 2; + for (i = 1; i < numStepsList->GetNumberOfIds(); i++) + { + numSteps += numStepsList->GetId(i); + if (timeStep > numSteps) + { + fileNum++; + timeStepInFile -= numStepsList->GetId(i); + } + } + } + if (this->FileSetFileNameNumbers->GetNumberOfItems() > 0) + { + int collectionNum = this->FileSetsWithFilenameNumbers-> + IsId(this->MeasuredFileSet); + if (collectionNum > -1) + { + filenameNumbers = this->FileSetFileNameNumbers-> + GetItem(fileSet); + filenameNum = filenameNumbers->GetId(fileNum-1); + this->ReplaceWildcards(fileName, filenameNum); + } + } + } + } + } + if (!this->ReadMeasuredGeometryFile(fileName, timeStepInFile, this->Cache)) + { + vtkErrorMacro("error reading measured geometry file"); + delete [] fileName; + return 0; + } + delete [] fileName; + } + this->CacheMTime.Modified(); + } + output->ShallowCopy(this->Cache); + + if ((this->NumberOfVariables + this->NumberOfComplexVariables) > 0) + { + if (!this->ReadVariableFiles(output)) + { + vtkErrorMacro("error reading variable files"); + return 0; + } + } + + return 1; +} diff --git a/src/Plugins/StaticMesh/vtkStaticEnSightGoldReader.h b/src/Plugins/StaticMesh/vtkStaticEnSightGoldReader.h new file mode 100644 index 00000000..1760e19f --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticEnSightGoldReader.h @@ -0,0 +1,69 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticEnSightGoldReader.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/** + * @class vtkStaticEnSightGoldReader + * @brief class to read binary EnSight6 files + * + * vtkStaticEnSightGoldReader is a class to read binary EnSight6 files into vtk. + * Because the different parts of the EnSight data can be of various data + * types, this reader produces multiple outputs, one per part in the input + * file. + * All variable information is being stored in field data. The descriptions + * listed in the case file are used as the array names in the field data. + * For complex vector variables, the description is appended with _r (for the + * array of real values) and _i (for the array if imaginary values). Complex + * scalar variables are stored as a single array with 2 components, real and + * imaginary, listed in that order. + * @warning + * You must manually call Update on this reader and then connect the rest + * of the pipeline because (due to the nature of the file format) it is + * not possible to know ahead of time how many outputs you will have or + * what types they will be. + * This reader can only handle static EnSight datasets (both static geometry + * and variables). +*/ + +#ifndef vtkStaticEnSightGoldReader_h +#define vtkStaticEnSightGoldReader_h + +#include "vtkEnSightGoldReader.h" +#include "vtkNew.h" + +class vtkMultiBlockDataSet; + +class vtkStaticEnSightGoldReader : public vtkEnSightGoldReader +{ +public: + static vtkStaticEnSightGoldReader *New(); + vtkTypeMacro(vtkStaticEnSightGoldReader, vtkEnSightGoldReader); + +protected: + vtkStaticEnSightGoldReader() = default; + ~vtkStaticEnSightGoldReader() override = default; + + int RequestData(vtkInformation*, + vtkInformationVector**, + vtkInformationVector*) override; + + vtkNew Cache; + vtkTimeStamp CacheMTime; + +private: + vtkStaticEnSightGoldReader(const vtkStaticEnSightGoldReader&) = delete; + void operator=(const vtkStaticEnSightGoldReader&) = delete; +}; + +#endif + diff --git a/src/Plugins/StaticMesh/vtkStaticMeshObjectFactory.cxx b/src/Plugins/StaticMesh/vtkStaticMeshObjectFactory.cxx new file mode 100644 index 00000000..4c6f1e6e --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticMeshObjectFactory.cxx @@ -0,0 +1,132 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticDataSetSurfaceFilter.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +#include "vtkStaticMeshObjectFactory.h" + +#include +#include +#include +#include + +#include "vtkStaticDataSetSurfaceFilter.h" +#include "vtkStaticEnSight6BinaryReader.h" +#include "vtkStaticEnSight6Reader.h" +#include "vtkStaticEnSightGoldBinaryReader.h" +#include "vtkStaticEnSightGoldReader.h" +#include "vtkStaticPlaneCutter.h" + +#ifdef PARAVIEW_USE_MPI +#include "vtkStaticPUnstructuredGridGhostCellsGenerator.h" +#endif + +vtkStandardNewMacro(vtkStaticMeshObjectFactory); + +VTK_CREATE_CREATE_FUNCTION(vtkStaticDataSetSurfaceFilter); +VTK_CREATE_CREATE_FUNCTION(vtkStaticPlaneCutter); +VTK_CREATE_CREATE_FUNCTION(vtkStaticEnSight6BinaryReader); +VTK_CREATE_CREATE_FUNCTION(vtkStaticEnSight6Reader); +VTK_CREATE_CREATE_FUNCTION(vtkStaticEnSightGoldReader); +VTK_CREATE_CREATE_FUNCTION(vtkStaticEnSightGoldBinaryReader); + +#ifdef PARAVIEW_USE_MPI +VTK_CREATE_CREATE_FUNCTION(vtkStaticPUnstructuredGridGhostCellsGenerator); +#endif + +vtkStaticMeshObjectFactory::vtkStaticMeshObjectFactory() +{ + vtkDebugMacro("Create vtkStaticMeshObjectFactory"); + + this->RegisterOverride("vtkDataSetSurfaceFilter", "vtkStaticDataSetSurfaceFilter", + "StaticDataSetSurfaceFilter", 1, vtkObjectFactoryCreatevtkStaticDataSetSurfaceFilter); + this->RegisterOverride("vtkPlaneCutter", "vtkStaticPlaneCutter", "StaticPlaneCutter", 1, + vtkObjectFactoryCreatevtkStaticPlaneCutter); + this->RegisterOverride("vtkEnSight6BinaryReader", "vtkStaticEnSight6BinaryReader", "StaticEnSight6BinaryReader", 1, + vtkObjectFactoryCreatevtkStaticEnSight6BinaryReader); + this->RegisterOverride("vtkEnSight6Reader", "vtkStaticEnSight6Reader", "StaticEnSight6Reader", 1, + vtkObjectFactoryCreatevtkStaticEnSight6Reader); + this->RegisterOverride("vtkEnSightGoldReader", "vtkStaticEnSight6BinaryReader", "StaticEnSight6BinaryReader", 1, + vtkObjectFactoryCreatevtkStaticEnSightGoldReader); + this->RegisterOverride("vtkEnSightGoldBinaryReader", "vtkStaticEnSightGoldBinaryReader", "StaticEnSightGoldBinaryReader", 1, + vtkObjectFactoryCreatevtkStaticEnSightGoldBinaryReader); + +#ifdef PARAVIEW_USE_MPI + this->RegisterOverride("vtkPUnstructuredGridGhostCellsGenerator", + "vtkStaticPUnstructuredGridGhostCellsGenerator", "StaticPUnstructuredGridGhostCellsGenerator", + 1, vtkObjectFactoryCreatevtkStaticPUnstructuredGridGhostCellsGenerator); +#endif +} + +vtkStaticMeshObjectFactory::~vtkStaticMeshObjectFactory() +{ + vtkDebugMacro("Delete vtkStaticMeshObjectFactory"); +} + +void vtkStaticMeshObjectFactory::PrintSelf(ostream& os, vtkIndent indent) +{ + os << indent << "VTK Static Mesh Extension Factory" << endl; +} + +const char* vtkStaticMeshObjectFactory::GetVTKSourceVersion() +{ + return VTK_SOURCE_VERSION; +} + +const char* vtkStaticMeshObjectFactory::GetDescription() +{ + return "VTK Static Mesh Extension Factory"; +} + +class StaticFactoryInitialize +{ +public: + StaticFactoryInitialize() + { + bool hasStaticPluginFactory = false; + vtkObjectFactoryCollection* collection = vtkObjectFactory::GetRegisteredFactories(); + collection->InitTraversal(); + vtkObjectFactory* f = collection->GetNextItem(); + while (f) + { + if (f->IsA("vtkStaticMeshObjectFactory")) + { + hasStaticPluginFactory = true; + break; + } + f = collection->GetNextItem(); + } + if (!hasStaticPluginFactory) + { + vtkStaticMeshObjectFactory* instance = vtkStaticMeshObjectFactory::New(); + vtkObjectFactory::RegisterFactory(instance); + instance->Delete(); + } + } + + virtual ~StaticFactoryInitialize() + { + vtkObjectFactoryCollection* collection = vtkObjectFactory::GetRegisteredFactories(); + collection->InitTraversal(); + vtkObjectFactory* f; + while ((f = collection->GetNextItem())) + { + if (f->IsA("vtkStaticMeshObjectFactory")) + { + vtkObjectFactory::UnRegisterFactory(f); + break; + } + } + } +}; + +static StaticFactoryInitialize StaticFactory; diff --git a/src/Plugins/StaticMesh/vtkStaticMeshObjectFactory.h b/src/Plugins/StaticMesh/vtkStaticMeshObjectFactory.h new file mode 100644 index 00000000..8af46d5f --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticMeshObjectFactory.h @@ -0,0 +1,58 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticDataSetSurfaceFilter.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/** + * @class vtkStaticMeshObjectFactory + * @brief Generate static version of dataset and filter for statix mesh plugin + * + * vtkStaticMeshObjectFactory is a vtk object factory, instantiating static version + * of some dataset and filters. +*/ + +#ifndef vtkStaticMeshObjectFactory_h +#define vtkStaticMeshObjectFactory_h + +#include "vtkObjectFactory.h" // Must be included before singletons + +class vtkStaticMeshObjectFactory : public vtkObjectFactory +{ +public: + vtkTypeMacro(vtkStaticMeshObjectFactory, vtkObjectFactory); + static vtkStaticMeshObjectFactory* New(); + void PrintSelf(ostream& os, vtkIndent indent) VTK_OVERRIDE; + + /** + * All sub-classes of vtkObjectFactory must return the version of + * VTK they were built with. This should be implemented with the macro + * VTK_SOURCE_VERSION and NOT a call to vtkVersion::GetVTKSourceVersion. + * As the version needs to be compiled into the file as a string constant. + * This is critical to determine possible incompatible dynamic factory loads. + */ + const char* GetVTKSourceVersion() VTK_OVERRIDE; + + /** + * Return a descriptive string describing the factory. + */ + const char* GetDescription() VTK_OVERRIDE; + +protected: + vtkStaticMeshObjectFactory(); + ~vtkStaticMeshObjectFactory() VTK_OVERRIDE; + +private: + vtkStaticMeshObjectFactory(const vtkStaticMeshObjectFactory&) VTK_DELETE_FUNCTION; + void operator=(const vtkStaticMeshObjectFactory&) VTK_DELETE_FUNCTION; +}; + +#endif diff --git a/src/Plugins/StaticMesh/vtkStaticPUnstructuredGridGhostCellsGenerator.cxx b/src/Plugins/StaticMesh/vtkStaticPUnstructuredGridGhostCellsGenerator.cxx new file mode 100644 index 00000000..b4ab07c9 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticPUnstructuredGridGhostCellsGenerator.cxx @@ -0,0 +1,517 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticPUnstructuredGridGhostCellsGenerator.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +#include "vtkStaticPUnstructuredGridGhostCellsGenerator.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +static const int SUGGCG_SIZE_EXCHANGE_TAG = 9002; +static const int SUGGCG_DATA_EXCHANGE_TAG = 9003; + +vtkStandardNewMacro(vtkStaticPUnstructuredGridGhostCellsGenerator); + +//---------------------------------------------------------------------------- +vtkStaticPUnstructuredGridGhostCellsGenerator::vtkStaticPUnstructuredGridGhostCellsGenerator() +{ + this->InputMeshTime = 0; + this->FilterMTime = 0; + + vtkMPIController* controller = + vtkMPIController::SafeDownCast(vtkMultiProcessController::GetGlobalController()); + if (controller) + { + // Initialise vtkIdList vectors + this->GhostPointsToReceive.resize(controller->GetNumberOfProcesses()); + this->GhostPointsToSend.resize(controller->GetNumberOfProcesses()); + this->GhostCellsToReceive.resize(controller->GetNumberOfProcesses()); + this->GhostCellsToSend.resize(controller->GetNumberOfProcesses()); + + int nProc = controller->GetNumberOfProcesses(); + + for (int i = 0; i < nProc; i++) + { + this->GhostCellsToReceive[i] = vtkSmartPointer::New(); + this->GhostCellsToSend[i] = vtkSmartPointer::New(); + this->GhostPointsToReceive[i] = vtkSmartPointer::New(); + this->GhostPointsToSend[i] = vtkSmartPointer::New(); + } + } +} + +//---------------------------------------------------------------------------- +vtkStaticPUnstructuredGridGhostCellsGenerator::~vtkStaticPUnstructuredGridGhostCellsGenerator() +{ +} + +//----------------------------------------------------------------------------- +int vtkStaticPUnstructuredGridGhostCellsGenerator::RequestData( + vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector) +{ + // get the inputs and outputs + vtkInformation* inInfo = inputVector[0]->GetInformationObject(0); + vtkInformation* outInfo = outputVector->GetInformationObject(0); + vtkUnstructuredGridBase* input = + vtkUnstructuredGridBase::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT())); + vtkUnstructuredGrid* output = + vtkUnstructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())); + + // Recover the static unstructured grid + vtkUnstructuredGrid* inputUG = vtkUnstructuredGrid::SafeDownCast(input); + if (!inputUG) + { + // For any other types of input, fall back to superclass implementation + return this->Superclass::RequestData(request, inputVector, outputVector); + } + + // Check cache validity + if (this->InputMeshTime == inputUG->GetMeshMTime() && this->FilterMTime == this->GetMTime()) + { + // Cache mesh is up to date, use it to generate data + // Update the cache data + this->UpdateCacheData(input); + + // Copy the updated cache into the output + output->ShallowCopy(this->Cache.Get()); + return 1; + } + else + { + // Add Arrays Ids needed + vtkNew tmpInput; + this->AddIdsArrays(input, tmpInput.Get()); + + // Create an input vector to pass the completed input to the superclass + // RequestData method + vtkNew tmpInputVec; + tmpInputVec->Copy(inputVector[0], 1); + vtkInformation* tmpInInfo = tmpInputVec->GetInformationObject(0); + tmpInInfo->Set(vtkDataObject::DATA_OBJECT(), tmpInput.Get()); + vtkInformationVector* tmpInputVecPt = tmpInputVec.Get(); + int ret = this->Superclass::RequestData(request, &tmpInputVecPt, outputVector); + + // Update the cache with superclass output + this->Cache->ShallowCopy(output); + this->InputMeshTime = inputUG->GetMeshMTime(); + this->FilterMTime = this->GetMTime(); + + this->ProcessGhostIds(); + + return ret; + } +} + +//----------------------------------------------------------------------------- +void vtkStaticPUnstructuredGridGhostCellsGenerator::ProcessGhostIds() +{ + vtkMPIController* controller = + vtkMPIController::SafeDownCast(vtkMultiProcessController::GetGlobalController()); + if (controller) + { + int nProc = controller->GetNumberOfProcesses(); + int rank = controller->GetLocalProcessId(); + + // Clear ghost vectors + for (int i = 0; i < nProc; i++) + { + this->GhostCellsToReceive[i]->SetNumberOfIds(0); + this->GhostCellsToSend[i]->SetNumberOfIds(0); + this->GhostPointsToReceive[i]->SetNumberOfIds(0); + this->GhostPointsToSend[i]->SetNumberOfIds(0); + } + + // Recover all needed data arrays + vtkPointData* cachePD = this->Cache->GetPointData(); + vtkCellData* cacheCD = this->Cache->GetCellData(); + vtkUnsignedCharArray* pointGhostArray = this->Cache->GetPointGhostArray(); + vtkUnsignedCharArray* cellGhostArray = this->Cache->GetCellGhostArray(); + vtkIdTypeArray* pointIds = vtkIdTypeArray::SafeDownCast(cachePD->GetAbstractArray("Ids")); + vtkIdTypeArray* cellIds = vtkIdTypeArray::SafeDownCast(cacheCD->GetAbstractArray("Ids")); + vtkIntArray* pointProcIds = vtkIntArray::SafeDownCast(cachePD->GetAbstractArray("ProcessId")); + vtkIntArray* cellProcIds = vtkIntArray::SafeDownCast(cacheCD->GetAbstractArray("ProcessId")); + if (!pointGhostArray || !pointIds || !pointProcIds || !cellGhostArray || !cellIds || + !cellProcIds) + { + // Sanity check + vtkWarningMacro("Arrays are missing from cache, cache is discarded"); + this->InputMeshTime = 0; + this->FilterMTime = 0; + } + else + { + // Compute list of remote ghost point ids + // and corresponding local point ids. + std::vector > remoteGhostPoints; + remoteGhostPoints.resize(nProc); + for (vtkIdType i = 0; i < pointGhostArray->GetNumberOfTuples(); i++) + { + if (pointGhostArray->GetValue(i) != 0) + { + this->GhostPointsToReceive[pointProcIds->GetValue(i)]->InsertNextId(i); + remoteGhostPoints[pointProcIds->GetValue(i)].push_back(pointIds->GetValue(i)); + } + } + + // Compute list of remote ghost cell ids + // and corresponding local cell ids. + std::vector > remoteGhostCells; + remoteGhostCells.resize(nProc); + for (vtkIdType i = 0; i < cellGhostArray->GetNumberOfTuples(); i++) + { + if (cellGhostArray->GetValue(i) != 0) + { + this->GhostCellsToReceive[cellProcIds->GetValue(i)]->InsertNextId(i); + remoteGhostCells[cellProcIds->GetValue(i)].push_back(cellIds->GetValue(i)); + } + } + + // Send requested ghost point ids to their own rank + vtkMPICommunicator::Request pointSizeReqs[nProc]; + vtkMPICommunicator::Request pointIdsReqs[nProc]; + vtkIdType lengths[nProc]; + for (int i = 0; i < nProc; i++) + { + if (i != rank) + { + lengths[i] = remoteGhostPoints[i].size(); + controller->NoBlockSend(&lengths[i], 1, i, SUGGCG_SIZE_EXCHANGE_TAG, pointSizeReqs[i]); + controller->NoBlockSend(&remoteGhostPoints[i][0], remoteGhostPoints[i].size(), i, + SUGGCG_DATA_EXCHANGE_TAG, pointIdsReqs[i]); + } + } + + // Receive and store requested ghost point ids. + for (int i = 0; i < nProc; i++) + { + if (i != rank) + { + vtkIdType length; + controller->Receive(&length, 1, i, SUGGCG_SIZE_EXCHANGE_TAG); + this->GhostPointsToSend[i]->SetNumberOfIds(length); + controller->Receive( + this->GhostPointsToSend[i]->GetPointer(0), length, i, SUGGCG_DATA_EXCHANGE_TAG); + } + } + controller->Barrier(); + + // Send requested ghost cell ids to their own rank + vtkMPICommunicator::Request cellSizeReqs[nProc]; + vtkMPICommunicator::Request cellIdsReqs[nProc]; + for (int i = 0; i < nProc; i++) + { + if (i != rank) + { + lengths[i] = remoteGhostCells[i].size(); + controller->NoBlockSend(&lengths[i], 1, i, SUGGCG_SIZE_EXCHANGE_TAG, cellSizeReqs[i]); + controller->NoBlockSend( + &remoteGhostCells[i][0], lengths[i], i, SUGGCG_DATA_EXCHANGE_TAG, cellIdsReqs[i]); + } + } + // Receive and store requested ghost cell ids. + for (int i = 0; i < nProc; i++) + { + if (i != rank) + { + vtkIdType length; + controller->Receive(&length, 1, i, SUGGCG_SIZE_EXCHANGE_TAG); + this->GhostCellsToSend[i]->SetNumberOfIds(length); + controller->Receive( + this->GhostCellsToSend[i]->GetPointer(0), length, i, SUGGCG_DATA_EXCHANGE_TAG); + } + } + controller->Barrier(); + } + } +} + +//----------------------------------------------------------------------------- +void vtkStaticPUnstructuredGridGhostCellsGenerator::AddIdsArrays( + vtkDataSet* input, vtkDataSet* output) +{ + vtkDataSet* tmpInput; + tmpInput = input; + vtkNew generateIdScalars; + vtkNew processPointIdScalars; + vtkNew processCellIdScalars; + + // Check for Ids array + vtkAbstractArray* pointIdsTmp = input->GetPointData()->GetAbstractArray("Ids"); + vtkAbstractArray* cellIdsTmp = input->GetCellData()->GetAbstractArray("Ids"); + if (!pointIdsTmp || !cellIdsTmp) + { + // Create Ids array + generateIdScalars->SetInputData(tmpInput); + generateIdScalars->SetIdsArrayName("Ids"); + generateIdScalars->Update(); + tmpInput = generateIdScalars->GetOutput(); + } + + // Check for ProcessId point array + vtkAbstractArray* procIdsTmp = input->GetPointData()->GetAbstractArray("ProcessId"); + if (!procIdsTmp) + { + // Create ProcessId point Array + processPointIdScalars->SetInputData(tmpInput); + processPointIdScalars->SetScalarModeToPointData(); + processPointIdScalars->Update(); + tmpInput = processPointIdScalars->GetOutput(); + } + + // Check for ProcessId Cell Array + procIdsTmp = input->GetCellData()->GetAbstractArray("ProcessId"); + if (!procIdsTmp) + { + // Create ProcessId Cell array + vtkNew processIdScalars; + processCellIdScalars->SetInputData(tmpInput); + processCellIdScalars->SetScalarModeToCellData(); + processCellIdScalars->Update(); + tmpInput = processCellIdScalars->GetOutput(); + } + output->ShallowCopy(tmpInput); +} + +//----------------------------------------------------------------------------- +void vtkStaticPUnstructuredGridGhostCellsGenerator::UpdateCacheData(vtkDataSet* input) +{ + this->UpdateCacheDataWithInput(input); + this->UpdateCacheGhostCellAndPointData(input); +} + +//----------------------------------------------------------------------------- +void vtkStaticPUnstructuredGridGhostCellsGenerator::UpdateCacheDataWithInput(vtkDataSet* input) +{ + // Recover point and cell data + vtkPointData* cachePD = this->Cache->GetPointData(); + vtkCellData* cacheCD = this->Cache->GetCellData(); + vtkPointData* inPD = input->GetPointData(); + vtkCellData* inCD = input->GetCellData(); + + // Update cache point data using input point data + // Of course this concerns only non-ghost points + for (int i = 0; i < inPD->GetNumberOfArrays(); i++) + { + vtkAbstractArray* cacheArray = cachePD->GetAbstractArray(inPD->GetArrayName(i)); + if (cacheArray) + { + cacheArray->InsertTuples(0, input->GetNumberOfPoints(), 0, inPD->GetAbstractArray(i)); + } + } + + // Update cache cell data using input cell data + // Of course this concerns only non-ghost cells + for (int i = 0; i < inCD->GetNumberOfArrays(); i++) + { + vtkAbstractArray* cacheArray = cacheCD->GetAbstractArray(inCD->GetArrayName(i)); + if (cacheArray) + { + cacheArray->InsertTuples(0, input->GetNumberOfCells(), 0, inCD->GetAbstractArray(i)); + } + } + + // Update field data + this->Cache->GetFieldData()->ShallowCopy(input->GetFieldData()); +} + +//----------------------------------------------------------------------------- +void vtkStaticPUnstructuredGridGhostCellsGenerator::UpdateCacheGhostCellAndPointData( + vtkDataSet* input) +{ + vtkMPIController* controller = + vtkMPIController::SafeDownCast(vtkMultiProcessController::GetGlobalController()); + if (controller) + { + vtkPointData* cachePD = this->Cache->GetPointData(); + vtkCellData* cacheCD = this->Cache->GetCellData(); + vtkPointData* inPD = input->GetPointData(); + vtkCellData* inCD = input->GetCellData(); + + int nProc = controller->GetNumberOfProcesses(); + int rank = controller->GetLocalProcessId(); + + vtkNew buffers[nProc]; + vtkIdType lengths[nProc]; + vtkMPICommunicator::Request sizeReqs[nProc]; + vtkMPICommunicator::Request dataReqs[nProc]; + + // For each rank + for (int i = 0; i < nProc; i++) + { + if (i != rank && this->GhostPointsToSend[i]->GetNumberOfIds() > 0) + { + // Prepare ghost points point data to send it as a table + vtkNew ghostPointData; + ghostPointData->CopyAllocate(inPD); + + // Prepare a list of iterating ids and copy all ghost point ids + // for this rank into the ghostPointData + vtkNew dumStaticPointIds; + dumStaticPointIds->SetNumberOfIds(this->GhostPointsToSend[i]->GetNumberOfIds()); + for (vtkIdType id = 0; id < this->GhostPointsToSend[i]->GetNumberOfIds(); id++) + { + dumStaticPointIds->SetId(id, id); + } + ghostPointData->CopyData(inPD, this->GhostPointsToSend[i].Get(), dumStaticPointIds.Get()); + + // Add each point data array to a dumStatic table + vtkNew pointDataTable; + for (int iArr = 0; iArr < ghostPointData->GetNumberOfArrays(); iArr++) + { + pointDataTable->AddColumn(ghostPointData->GetArray(iArr)); + } + + // Marshall the table and transfer it to rank + vtkCommunicator::MarshalDataObject(pointDataTable.Get(), buffers[i].Get()); + lengths[i] = buffers[i]->GetNumberOfTuples(); + controller->NoBlockSend(&lengths[i], 1, i, SUGGCG_SIZE_EXCHANGE_TAG, sizeReqs[i]); + controller->NoBlockSend((char*)(buffers[i]->GetVoidPointer(0)), lengths[i], i, + SUGGCG_DATA_EXCHANGE_TAG, dataReqs[i]); + } + } + // Foe each rank + for (int i = 0; i < nProc; i++) + { + if (i != rank && this->GhostPointsToReceive[i]->GetNumberOfIds() > 0) + { + // Receive dumStatic table to unmarshall + vtkIdType length; + controller->Receive(&length, 1, i, SUGGCG_SIZE_EXCHANGE_TAG); + + vtkNew recvBuffer; + recvBuffer->SetNumberOfValues(length); + controller->Receive( + (char*)(recvBuffer->GetVoidPointer(0)), length, i, SUGGCG_DATA_EXCHANGE_TAG); + vtkNew pointDataTable; + vtkCommunicator::UnMarshalDataObject(recvBuffer.Get(), pointDataTable.Get()); + + // Create a dumStatic iterating point ids + vtkNew dumStaticPointIds; + dumStaticPointIds->SetNumberOfIds(this->GhostPointsToReceive[i]->GetNumberOfIds()); + for (vtkIdType id = 0; id < this->GhostPointsToReceive[i]->GetNumberOfIds(); id++) + { + dumStaticPointIds->SetId(id, id); + } + + // Copy the tuples of each array from the dumStatic table + // into the ghost point data + for (int iArr = 0; iArr < pointDataTable->GetNumberOfColumns(); iArr++) + { + vtkAbstractArray* arrayToCopyIn = + cachePD->GetAbstractArray(pointDataTable->GetColumnName(iArr)); + if (arrayToCopyIn) + { + arrayToCopyIn->InsertTuples(this->GhostPointsToReceive[i].Get(), + dumStaticPointIds.Get(), pointDataTable->GetColumn(iArr)); + } + } + } + } + // Make sure all rank finished + controller->Barrier(); + + for (int i = 0; i < nProc; i++) + { + if (i != rank && this->GhostCellsToSend[i]->GetNumberOfIds() > 0) + { + // Prepare ghost cells data to send it as a table + vtkNew ghostCellData; + ghostCellData->CopyAllocate(inCD); + + // Prepare a list of iterating ids and copy all ghost point ids + // for this rank into the ghostPointData + vtkNew dumStaticCellIds; + dumStaticCellIds->SetNumberOfIds(this->GhostCellsToSend[i]->GetNumberOfIds()); + for (vtkIdType id = 0; id < this->GhostCellsToSend[i]->GetNumberOfIds(); id++) + { + dumStaticCellIds->SetId(id, id); + } + ghostCellData->CopyData(inCD, this->GhostCellsToSend[i].Get(), dumStaticCellIds.Get()); + + // Add each point data array to a dumStatic table + vtkNew cellDataTable; + for (int iArr = 0; iArr < ghostCellData->GetNumberOfArrays(); iArr++) + { + cellDataTable->AddColumn(ghostCellData->GetArray(iArr)); + } + + // Marshall the table and transfer it to rank + vtkCommunicator::MarshalDataObject(cellDataTable.Get(), buffers[i].Get()); + lengths[i] = buffers[i]->GetNumberOfTuples(); + controller->NoBlockSend(&lengths[i], 1, i, SUGGCG_SIZE_EXCHANGE_TAG, sizeReqs[i]); + controller->NoBlockSend((char*)(buffers[i]->GetVoidPointer(0)), lengths[i], i, + SUGGCG_DATA_EXCHANGE_TAG, dataReqs[i]); + } + } + for (int i = 0; i < nProc; i++) + { + if (i != rank && this->GhostCellsToReceive[i]->GetNumberOfIds() > 0) + { + // Receive dumStatic table to unmarshall + vtkIdType length; + controller->Receive(&length, 1, i, SUGGCG_SIZE_EXCHANGE_TAG); + + vtkNew recvBuffer; + recvBuffer->SetNumberOfValues(length); + controller->Receive( + (char*)(recvBuffer->GetVoidPointer(0)), length, i, SUGGCG_DATA_EXCHANGE_TAG); + vtkNew cellDataTable; + vtkCommunicator::UnMarshalDataObject(recvBuffer.Get(), cellDataTable.Get()); + + // Create a dumStatic iterating point ids + vtkNew dumStaticCellIds; + dumStaticCellIds->SetNumberOfIds(this->GhostCellsToReceive[i]->GetNumberOfIds()); + for (vtkIdType id = 0; id < this->GhostCellsToReceive[i]->GetNumberOfIds(); id++) + { + dumStaticCellIds->SetId(id, id); + } + + // Copy the tuples of each array from the dumStatic table + // into the ghost point data + for (int iArr = 0; iArr < cellDataTable->GetNumberOfColumns(); iArr++) + { + vtkAbstractArray* arrayToCopyIn = + cacheCD->GetAbstractArray(cellDataTable->GetColumnName(iArr)); + if (arrayToCopyIn) + { + arrayToCopyIn->InsertTuples(this->GhostCellsToReceive[i].Get(), dumStaticCellIds.Get(), + cellDataTable->GetColumn(iArr)); + } + } + } + } + // Make sure all rank finished + controller->Barrier(); + } +} + +//---------------------------------------------------------------------------- +void vtkStaticPUnstructuredGridGhostCellsGenerator::PrintSelf(ostream& os, vtkIndent indent) +{ + this->Superclass::PrintSelf(os, indent); + os << indent << "Cache: " << this->Cache << endl; + os << indent << "Input Mesh Time: " << this->InputMeshTime << endl; + os << indent << "Filter mTime: " << this->FilterMTime << endl; +} diff --git a/src/Plugins/StaticMesh/vtkStaticPUnstructuredGridGhostCellsGenerator.h b/src/Plugins/StaticMesh/vtkStaticPUnstructuredGridGhostCellsGenerator.h new file mode 100644 index 00000000..f40653aa --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticPUnstructuredGridGhostCellsGenerator.h @@ -0,0 +1,102 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticPUnstructuredGridGhostCellsGenerator.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/** + * @class vtkStaticPUnstructuredGridGhostCellsGenerator + * @brief StaticMesh aware GhostCellGenerator implementation + * + * This class specialize vtkPUnstructuredGridGhostCellGenerator + * This class improves it by making it static mesh aware. + * The first time this filter is executed it will store its output + * in a cache as well as a list of ghost and point ids to request from other rank + * On next execution, if the mesh is static, it will uses the list of ids to request + * only point and cell data for the ghost point and cell from other + * allowing to update the output without needing to recompute everything + * + * @sa + * vtkPUnstructuredGridGhostCellsGenerator +*/ + +#ifndef vtkStaticPUnstructuredGridGhostCellsGenerator_h +#define vtkStaticPUnstructuredGridGhostCellsGenerator_h + +#include +#include +#include +#include + +#include + +class vtkUnstructuredGrid; + +class vtkStaticPUnstructuredGridGhostCellsGenerator : public vtkPUnstructuredGridGhostCellsGenerator +{ +public: + static vtkStaticPUnstructuredGridGhostCellsGenerator* New(); + typedef vtkPUnstructuredGridGhostCellsGenerator + Superclass; // vtkTypeMacro can't be used with a factory built object + void PrintSelf(ostream& os, vtkIndent indent) VTK_OVERRIDE; + + int RequestData(vtkInformation*, vtkInformationVector**, vtkInformationVector*) VTK_OVERRIDE; + +protected: + vtkStaticPUnstructuredGridGhostCellsGenerator(); + ~vtkStaticPUnstructuredGridGhostCellsGenerator() VTK_OVERRIDE; + + /** + * Check input for "ProcessId" and "Ids" point and cell array, + * If absent, compute and add them. + */ + static void AddIdsArrays(vtkDataSet* input, vtkDataSet* output); + + /** + * Using the Cache, exchange and update ghost point and cell + * ids between ranks + */ + virtual void ProcessGhostIds(); + + /** + * Update cache with input and with ghost cells info + */ + virtual void UpdateCacheData(vtkDataSet* input); + + /** + * Copy input point, cell and field data into cache + */ + virtual void UpdateCacheDataWithInput(vtkDataSet* input); + + /** + * Using Cached ghost cell and points info + * Update ghost cell and point data in cache + * by sending input point and cell data to other ranks + */ + virtual void UpdateCacheGhostCellAndPointData(vtkDataSet* input); + + vtkNew Cache; + vtkMTimeType InputMeshTime; + vtkMTimeType FilterMTime; + + std::vector > GhostCellsToReceive; + std::vector > GhostCellsToSend; + std::vector > GhostPointsToReceive; + std::vector > GhostPointsToSend; + +private: + // Hide these from the user and the compiler. + vtkStaticPUnstructuredGridGhostCellsGenerator( + const vtkStaticPUnstructuredGridGhostCellsGenerator&) VTK_DELETE_FUNCTION; + void operator=(const vtkStaticPUnstructuredGridGhostCellsGenerator&) VTK_DELETE_FUNCTION; +}; + +#endif diff --git a/src/Plugins/StaticMesh/vtkStaticPlaneCutter.cxx b/src/Plugins/StaticMesh/vtkStaticPlaneCutter.cxx new file mode 100644 index 00000000..68fe3e3d --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticPlaneCutter.cxx @@ -0,0 +1,250 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticPlaneCutter.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +#include "vtkStaticPlaneCutter.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +vtkStandardNewMacro(vtkStaticPlaneCutter); + +//---------------------------------------------------------------------------- +vtkStaticPlaneCutter::vtkStaticPlaneCutter() +{ + this->InputMeshTime = 0; + this->FilterMTime = 0; +} + +//---------------------------------------------------------------------------- +vtkStaticPlaneCutter::~vtkStaticPlaneCutter() +{ +} + +//----------------------------------------------------------------------------- +int vtkStaticPlaneCutter::RequestData( + vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector) +{ + // get the inputs and outputs + vtkInformation* inInfo = inputVector[0]->GetInformationObject(0); + vtkInformation* outInfo = outputVector->GetInformationObject(0); + vtkUnstructuredGrid* input = vtkUnstructuredGrid::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT())); + vtkMultiBlockDataSet* inputMB = vtkMultiBlockDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT())); + vtkMultiBlockDataSet* mb = + vtkMultiBlockDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())); + if (!mb) + { + vtkErrorMacro("Ouput information does not contain expected type of data object"); + return 0; + } + + // Recover the fisrt and only block so this works with single block mb + if (inputMB && inputMB->GetNumberOfBlocks() == 1) + { + input = vtkUnstructuredGrid::SafeDownCast(inputMB->GetBlock(0)); + } + + // Recover the static unstructured grid + if (!input) + { + // For any other type of input, fall back to superclass implementation + return this->Superclass::RequestData(request, inputVector, outputVector); + } + + // Check cache validity + if (this->InputMeshTime == input->GetMeshMTime() && this->FilterMTime == this->GetMTime()) + { + // Cache mesh is up to date, use it to generate data + if (this->InterpolateAttributes) + { + // Update the cache data + this->UpdateCacheData(input); + } + + // Copy the updated cache into the output + mb->SetBlock(0, this->Cache.Get()); + return 1; + } + else + { + // Cache is invalid + // Add needed Arrays + vtkNew tmpInput; + this->AddIdsArray(input, tmpInput.Get()); + + // Create an input vector to pass the completed input to the superclass + // RequestData method + vtkNew tmpInputVec; + tmpInputVec->Copy(inputVector[0], 1); + vtkInformation* tmpInInfo = tmpInputVec->GetInformationObject(0); + tmpInInfo->Set(vtkDataObject::DATA_OBJECT(), tmpInput.Get()); + vtkInformationVector* tmpInputVecPt = tmpInputVec.Get(); + int ret = this->Superclass::RequestData(request, &tmpInputVecPt, outputVector); + + // Update the cache with superclass output + vtkMultiPieceDataSet* output = vtkMultiPieceDataSet::SafeDownCast(mb->GetBlock(0)); + if (!output) + { + vtkErrorMacro("Output is not of expected type"); + return 0; + } + this->RemovePointData(output); + + this->Cache->ShallowCopy(output); + this->InputMeshTime = input->GetMeshMTime(); + this->FilterMTime = this->GetMTime(); + + // Compute the cell ids to be passed from the input to the cache + this->ComputeCellIds(); + + return ret; + } +} + +//----------------------------------------------------------------------------- +void vtkStaticPlaneCutter::AddIdsArray(vtkDataSet* input, vtkDataSet* output) +{ + vtkDataSet* tmpInput; + tmpInput = input; + vtkNew generateIdScalars; + + // Check for Ids array + vtkAbstractArray* cellIdsTmp = input->GetCellData()->GetAbstractArray("Ids"); + if (!cellIdsTmp) + { + // Create Ids array + generateIdScalars->SetInputData(tmpInput); + generateIdScalars->SetIdsArrayName("Ids"); + generateIdScalars->Update(); + tmpInput = generateIdScalars->GetOutput(); + } + output->ShallowCopy(tmpInput); +} + +//----------------------------------------------------------------------------- +void vtkStaticPlaneCutter::RemovePointData(vtkMultiPieceDataSet* dataset) +{ + // Iterate over each piece of the multipiece output + vtkCompositeDataIterator* iter = dataset->NewIterator(); + iter->SkipEmptyNodesOn(); + for (iter->GoToFirstItem(); !iter->IsDoneWithTraversal(); iter->GoToNextItem()) + { + vtkPolyData* slice = vtkPolyData::SafeDownCast(iter->GetCurrentDataObject()); + if (slice) + { + // Remove all point data array on each slice + vtkPointData* pd = slice->GetPointData(); + for (int i = pd->GetNumberOfArrays() - 1; i >= 0; i--) + { + pd->RemoveArray(i); + } + } + } + iter->Delete(); +} + +//----------------------------------------------------------------------------- +void vtkStaticPlaneCutter::ComputeCellIds() +{ + this->CellToCopyFrom.clear(); + this->CellToCopyTo.clear(); + + // Iterate over each piece of the multipiece output + vtkCompositeDataIterator* iter = this->Cache->NewIterator(); + iter->SkipEmptyNodesOn(); + for (iter->GoToFirstItem(); !iter->IsDoneWithTraversal(); iter->GoToNextItem()) + { + vtkPolyData* slice = vtkPolyData::SafeDownCast(iter->GetCurrentDataObject()); + if (slice) + { + // For each piece, recover the Ids of the cells sliced and the corresponding + // cellId in the slice + vtkSmartPointer cellIdsFrom = vtkSmartPointer::New(); + vtkSmartPointer cellIdsTo = vtkSmartPointer::New(); + this->CellToCopyFrom.push_back(cellIdsFrom); + this->CellToCopyTo.push_back(cellIdsTo); + + vtkIdTypeArray* ids = vtkIdTypeArray::SafeDownCast(slice->GetCellData()->GetArray("Ids")); + cellIdsFrom->SetNumberOfIds(ids->GetNumberOfValues()); + cellIdsTo->SetNumberOfIds(ids->GetNumberOfValues()); + for (vtkIdType i = 0; i < ids->GetNumberOfValues(); i++) + { + cellIdsFrom->SetId(i, ids->GetValue(i)); + cellIdsTo->SetId(i, i); + } + } + } + iter->Delete(); +} + +//----------------------------------------------------------------------------- +void vtkStaticPlaneCutter::UpdateCacheData(vtkDataSet* input) +{ + // Remove useless FieldData Array from multipiece + // Created by automatic pass data in pipeline + vtkFieldData* mpFieldData = this->Cache->GetFieldData(); + for (int i = mpFieldData->GetNumberOfArrays() - 1; i >= 0; i--) + { + mpFieldData->RemoveArray(i); + } + + // Recover cell data + // We may want to improve this and support PointData at some point + vtkCellData* inCD = input->GetCellData(); + + vtkCompositeDataIterator* iter = this->Cache->NewIterator(); + iter->SkipEmptyNodesOn(); + int i = 0; + for (iter->GoToFirstItem(); !iter->IsDoneWithTraversal(); iter->GoToNextItem()) + { + vtkPolyData* slice = vtkPolyData::SafeDownCast(iter->GetCurrentDataObject()); + if (slice) + { + vtkCellData* sliceCD = slice->GetCellData(); + for (int iArr = 0; iArr < inCD->GetNumberOfArrays(); iArr++) + { + // For each array of the cell data of each slice + vtkAbstractArray* arrayToCopyIn = sliceCD->GetAbstractArray(inCD->GetArrayName(iArr)); + if (arrayToCopyIn) + { + // Copy the tuples from the input cell ids to the slice cell ids + arrayToCopyIn->InsertTuples(this->CellToCopyTo[i].Get(), this->CellToCopyFrom[i].Get(), + inCD->GetAbstractArray(iArr)); + } + } + // Update field data + slice->GetFieldData()->PassData(input->GetFieldData()); + } + i++; + } + iter->Delete(); +} + +//---------------------------------------------------------------------------- +void vtkStaticPlaneCutter::PrintSelf(ostream& os, vtkIndent indent) +{ + this->Superclass::PrintSelf(os, indent); + os << indent << "Cache: " << this->Cache << endl; + os << indent << "Input Mesh Time: " << this->InputMeshTime << endl; + os << indent << "Filter mTime: " << this->FilterMTime << endl; +} diff --git a/src/Plugins/StaticMesh/vtkStaticPlaneCutter.h b/src/Plugins/StaticMesh/vtkStaticPlaneCutter.h new file mode 100644 index 00000000..badea399 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkStaticPlaneCutter.h @@ -0,0 +1,88 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkStaticPlaneCutter.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/** + * @class vtkStaticPlaneCutter + * @brief StaticMesh aware implementation of vtkPlaneCutter vtk vtkUnstructuredGrid + * + * This class specialize vtkPlaneCutter for vtkUnstructuredGrid input. + * It uses a cache when the associated data chage over time but not the geometry. + * In order to be able to update the cache, we keep track of cells ids + * when the cache is computed. + * Contrary to its parent, this class does not interpolate point data, + * only transmit cell data. + * + * + * @sa + * vtkPlaneCutter +*/ + +#ifndef vtkStaticPlaneCutter_h +#define vtkStaticPlaneCutter_h + +#include +#include +#include +#include + +#include + +class vtkMultiPieceDataSet; + +class vtkStaticPlaneCutter : public vtkPlaneCutter +{ +public: + static vtkStaticPlaneCutter* New(); + typedef vtkPlaneCutter Superclass; // vtkTypeMacro can't be used with a factory built object + void PrintSelf(ostream& os, vtkIndent indent) VTK_OVERRIDE; + + int RequestData(vtkInformation*, vtkInformationVector**, vtkInformationVector*) VTK_OVERRIDE; + +protected: + vtkStaticPlaneCutter(); + ~vtkStaticPlaneCutter() VTK_OVERRIDE; + + /** + * Check input for "Ids" cell array, if absent, compute and add them. + */ + static void AddIdsArray(vtkDataSet* input, vtkDataSet* output); + + /** + * Update cache point, cell and field data using input + */ + void UpdateCacheData(vtkDataSet* input); + + /** + * Compute the cells ids to be used when updating the cache later + */ + void ComputeCellIds(); + + /** + * Remove all point data array of a multipiece input with polydata leafs + */ + static void RemovePointData(vtkMultiPieceDataSet* dataset); + + vtkNew Cache; + std::vector > CellToCopyFrom; + std::vector > CellToCopyTo; + vtkMTimeType InputMeshTime; + vtkMTimeType FilterMTime; + +private: + // Hide these from the user and the compiler. + vtkStaticPlaneCutter(const vtkStaticPlaneCutter&) VTK_DELETE_FUNCTION; + void operator=(const vtkStaticPlaneCutter&) VTK_DELETE_FUNCTION; +}; + +#endif diff --git a/src/Plugins/StaticMesh/vtkTemporalUGWavelet.cxx b/src/Plugins/StaticMesh/vtkTemporalUGWavelet.cxx new file mode 100644 index 00000000..404ba4d0 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkTemporalUGWavelet.cxx @@ -0,0 +1,140 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkRTAnalyticSource.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +#include "vtkTemporalUGWavelet.h" + +#include "vtkCellData.h" +#include "vtkDataSetTriangleFilter.h" +#include "vtkFloatArray.h" +#include "vtkImageData.h" +#include "vtkInformation.h" +#include "vtkInformationVector.h" +#include "vtkMultiProcessController.h" +#include "vtkNew.h" +#include "vtkObjectFactory.h" +#include "vtkPointData.h" +#include "vtkRTAnalyticSource.h" +#include "vtkStreamingDemandDrivenPipeline.h" +#include "vtkUnstructuredGrid.h" + +vtkStandardNewMacro(vtkTemporalUGWavelet); + +// ---------------------------------------------------------------------------- +vtkTemporalUGWavelet::vtkTemporalUGWavelet() +{ + this->Cache = vtkUnstructuredGrid::New(); + this->NumberOfTimeSteps = 10; +} + +// ---------------------------------------------------------------------------- +vtkTemporalUGWavelet::~vtkTemporalUGWavelet() +{ + this->Cache->Delete(); +} + +//---------------------------------------------------------------------------- +int vtkTemporalUGWavelet::FillOutputPortInformation(int vtkNotUsed(port), vtkInformation* info) +{ + info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkUnstructuredGrid"); + return 1; +} + +// ---------------------------------------------------------------------------- +int vtkTemporalUGWavelet::RequestInformation( + vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector) +{ + vtkInformation* outInfo = outputVector->GetInformationObject(0); + double range[2] = { 0, static_cast(this->NumberOfTimeSteps - 1) }; + outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), range, 2); + double* outTimes = new double[this->NumberOfTimeSteps]; + for (int i = 0; i < this->NumberOfTimeSteps; i++) + { + outTimes[i] = i; + } + outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), outTimes, this->NumberOfTimeSteps); + outInfo->Set(CAN_HANDLE_PIECE_REQUEST(), 1); + return Superclass::RequestInformation(request, inputVector, outputVector); +} + +// ---------------------------------------------------------------------------- +int vtkTemporalUGWavelet::RequestData( + vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector) +{ + vtkInformation* outInfo = outputVector->GetInformationObject(0); + vtkUnstructuredGrid* data = vtkUnstructuredGrid::GetData(outInfo); + if (this->CacheMTime < this->GetMTime()) + { + // Create an output vector to recover the output image data + // RequestData method + vtkNew tmpOutputVec; + vtkNew image; + tmpOutputVec->Copy(outputVector, 1); + vtkInformation* tmpOutInfo = tmpOutputVec->GetInformationObject(0); + tmpOutInfo->Set(vtkDataObject::DATA_OBJECT(), image.Get()); + + // Generate wavelet + int ret = this->Superclass::RequestData(request, inputVector, tmpOutputVec.Get()); + if (ret == 0) + { + return ret; + } + + // Transform it to unstructured grid + vtkNew tetra; + tetra->SetInputData(image.Get()); + tetra->Update(); + + // Create the cache + this->Cache->ShallowCopy(tetra->GetOutput()); + this->CacheMTime.Modified(); + } + + // Use the cache + data->ShallowCopy(this->Cache); + + double t = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()); + + // Generate tPoint array + vtkIdType nbPoints = data->GetNumberOfPoints(); + vtkNew pointArray; + pointArray->SetName("tPoint"); + pointArray->SetNumberOfValues(nbPoints); + data->GetPointData()->AddArray(pointArray.Get()); + for (vtkIdType i = 0; i < nbPoints; i++) + { + pointArray->SetValue( + i, static_cast(i + t * (nbPoints / this->NumberOfTimeSteps)) % nbPoints); + } + + // Generate tCell array + vtkIdType nbCells = data->GetNumberOfCells(); + vtkNew cellArray; + cellArray->SetName("tCell"); + cellArray->SetNumberOfValues(nbCells); + data->GetCellData()->AddArray(cellArray.Get()); + + for (vtkIdType i = 0; i < nbCells; i++) + { + cellArray->SetValue(i, static_cast(i + t * (nbCells / this->NumberOfTimeSteps)) % nbCells); + } + + return 1; +} + +// ---------------------------------------------------------------------------- +void vtkTemporalUGWavelet::PrintSelf(ostream& os, vtkIndent indent) +{ + this->Superclass::PrintSelf(os, indent); + os << indent << "NumberOfTimeSteps: " << this->NumberOfTimeSteps << endl; +} diff --git a/src/Plugins/StaticMesh/vtkTemporalUGWavelet.h b/src/Plugins/StaticMesh/vtkTemporalUGWavelet.h new file mode 100644 index 00000000..e1a11965 --- /dev/null +++ b/src/Plugins/StaticMesh/vtkTemporalUGWavelet.h @@ -0,0 +1,67 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkTemporalUGWavelet.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/** + * @class vtkTemporalUGWavelet + * @brief Create a Temporal Unstructured Grid Wavelet with a static mesh + * + * vtkTemporalUGWavelet specialize vtkRTAnalyticSource to create + * a wavelet converted to vtkUnstructuredGrid, with timesteps. + * The "tPoint" and "tCell" arrays are the only data actually changing over time + * make the output a static mesh with data evolving over time. +*/ + +#ifndef vtkTemporalUGWavelet_h +#define vtkTemporalUGWavelet_h + +#include "vtkRTAnalyticSource.h" + +class vtkUnstructuredGrid; + +class vtkTemporalUGWavelet : public vtkRTAnalyticSource +{ +public: + static vtkTemporalUGWavelet* New(); + vtkTypeMacro(vtkTemporalUGWavelet, vtkRTAnalyticSource); + void PrintSelf(ostream& os, vtkIndent indent) VTK_OVERRIDE; + + //@{ + /** + * Set/Get the number of time steps. Initial value is 10. + */ + vtkSetMacro(NumberOfTimeSteps, int); + vtkGetMacro(NumberOfTimeSteps, int); + //@} + +protected: + vtkTemporalUGWavelet(); + ~vtkTemporalUGWavelet(); + + int FillOutputPortInformation(int vtkNotUsed(port), vtkInformation* info) VTK_OVERRIDE; + + int RequestInformation(vtkInformation* request, vtkInformationVector** inputVector, + vtkInformationVector* outputVector) VTK_OVERRIDE; + + int RequestData(vtkInformation*, vtkInformationVector**, vtkInformationVector*) VTK_OVERRIDE; + + int NumberOfTimeSteps; + vtkUnstructuredGrid* Cache; + vtkTimeStamp CacheMTime; + +private: + vtkTemporalUGWavelet(const vtkTemporalUGWavelet&) VTK_DELETE_FUNCTION; + void operator=(const vtkTemporalUGWavelet&) VTK_DELETE_FUNCTION; +}; + +#endif diff --git a/src/Plugins/VoroGauss/CMakeLists.txt b/src/Plugins/VoroGauss/CMakeLists.txt new file mode 100644 index 00000000..ac198030 --- /dev/null +++ b/src/Plugins/VoroGauss/CMakeLists.txt @@ -0,0 +1,66 @@ +# Copyright (C) 2017 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +PROJECT(VoroGauss) +CMAKE_MINIMUM_REQUIRED(VERSION 2.8) + +CMAKE_POLICY(SET CMP0003 NEW) +IF(${CMAKE_VERSION} VERSION_GREATER "3.0.0") + CMAKE_POLICY(SET CMP0022 OLD) + CMAKE_POLICY(SET CMP0023 OLD) +ENDIF() + +SET(MED_READER_VERSION "0.0.0") + +# Common CMake macros +# =================== +SET(CONFIGURATION_ROOT_DIR $ENV{CONFIGURATION_ROOT_DIR} CACHE PATH "Path to the Salome CMake configuration files") +IF(EXISTS ${CONFIGURATION_ROOT_DIR}) + LIST(APPEND CMAKE_MODULE_PATH "${CONFIGURATION_ROOT_DIR}/cmake") + INCLUDE(SalomeMacros) +ELSE() + MESSAGE(FATAL_ERROR "We absolutely need the Salome CMake configuration files, please define CONFIGURATION_ROOT_DIR !") +ENDIF() +FIND_PACKAGE(SalomePythonInterp REQUIRED) +FIND_PACKAGE(SalomePythonLibs REQUIRED) + +FIND_PACKAGE(ParaView REQUIRED) +IF(NOT ParaView_FOUND) + MESSAGE(FATAL_ERROR "Please locate ParaView." ) +ENDIF(NOT ParaView_FOUND) +INCLUDE(${PARAVIEW_USE_FILE}) + +SET(MEDCOUPLING_ROOT_DIR $ENV{MEDCOUPLING_ROOT_DIR} CACHE PATH "MEDCOUPLING_ROOT_DIR") +LIST(APPEND CMAKE_MODULE_PATH "${MEDCOUPLING_ROOT_DIR}/cmake_files") +FIND_PACKAGE(SalomeMEDCoupling REQUIRED) + +OPTION(BUILD_SHARED_LIBS "Build with shared libraries." ${VTK_BUILD_SHARED_LIBS}) + +SET(VTK_INSTALL_RUNTIME_DIR lib) +SET(VTK_INSTALL_LIBRARY_DIR lib) +SET(VTK_INSTALL_ARCHIVE_DIR lib) + +PV_PROCESS_MODULES() + +INCLUDE_DIRECTORIES( + ${MEDCOUPLING_INCLUDE_DIRS} +) + +ADD_SUBDIRECTORY(ParaViewPlugin) diff --git a/src/Plugins/VoroGauss/IO/vtkVoroGauss.cxx b/src/Plugins/VoroGauss/IO/vtkVoroGauss.cxx new file mode 100644 index 00000000..0cfc9e06 --- /dev/null +++ b/src/Plugins/VoroGauss/IO/vtkVoroGauss.cxx @@ -0,0 +1,898 @@ +// Copyright (C) 2017 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 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 +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay (EDF R&D) + +#include "vtkVoroGauss.h" + +#include "vtkAdjacentVertexIterator.h" +#include "vtkIntArray.h" +#include "vtkCellData.h" +#include "vtkPointData.h" +#include "vtkCellType.h" +#include "vtkCell.h" +#include "vtkCellArray.h" +#include "vtkIdTypeArray.h" + +#include "vtkStreamingDemandDrivenPipeline.h" +#include "vtkUnstructuredGrid.h" +#include "vtkMultiBlockDataSet.h" + +#include "vtkInformationStringKey.h" +#include "vtkAlgorithmOutput.h" +#include "vtkObjectFactory.h" +#include "vtkMutableDirectedGraph.h" +#include "vtkMultiBlockDataSet.h" +#include "vtkDataSet.h" +#include "vtkInformationVector.h" +#include "vtkInformation.h" +#include "vtkDataArraySelection.h" +#include "vtkTimeStamp.h" +#include "vtkInEdgeIterator.h" +#include "vtkInformationDataObjectKey.h" +#include "vtkInformationDataObjectMetaDataKey.h" +#include "vtkInformationDoubleVectorKey.h" +#include "vtkExecutive.h" +#include "vtkVariantArray.h" +#include "vtkStringArray.h" +#include "vtkDoubleArray.h" +#include "vtkFloatArray.h" +#include "vtkCharArray.h" +#include "vtkLongArray.h" +#include "vtkUnsignedCharArray.h" +#include "vtkDataSetAttributes.h" +#include "vtkDemandDrivenPipeline.h" +#include "vtkDataObjectTreeIterator.h" +#include "vtkWarpScalar.h" +#include "vtkQuadratureSchemeDefinition.h" +#include "vtkInformationQuadratureSchemeDefinitionVectorKey.h" +#include "vtkCompositeDataToUnstructuredGridFilter.h" +#include "vtkMultiBlockDataGroupFilter.h" + +#include "MEDCouplingMemArray.hxx" +#include "MEDCouplingUMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" +#include "InterpKernelAutoPtr.hxx" +#include "InterpKernelGaussCoords.hxx" + +#include +#include +#include +#include + +using MEDCoupling::DataArray; +using MEDCoupling::DataArrayInt; +using MEDCoupling::DataArrayDouble; +using MEDCoupling::MEDCouplingMesh; +using MEDCoupling::MEDCouplingUMesh; +using MEDCoupling::DynamicCastSafe; +using MEDCoupling::MEDCouplingFieldDouble; +using MEDCoupling::ON_GAUSS_PT; +using MEDCoupling::MCAuto; + +vtkStandardNewMacro(vtkVoroGauss); +/////////////////// + +std::map ComputeMapOfType() +{ + std::map ret; + int nbOfTypesInMC(sizeof(MEDCouplingUMesh::MEDCOUPLING2VTKTYPETRADUCER)/sizeof(int)); + for(int i=0;i ComputeRevMapOfType() +{ + std::map ret; + int nbOfTypesInMC(sizeof(MEDCouplingUMesh::MEDCOUPLING2VTKTYPETRADUCER)/sizeof(int)); + for(int i=0;ihasKey(ZE_KEY)) + return 0; + std::string ptSt(gd->value(ZE_KEY)); + void *pt(0); + std::istringstream iss(ptSt); iss >> pt; + return reinterpret_cast(pt); +} + +bool IsInformationOK(vtkInformation *info, std::vector& data) +{ + vtkInformationDoubleVectorKey *key(GetMEDReaderMetaDataIfAny()); + if(!key) + return false; + // Check the information contain meta data key + if(!info->Has(key)) + return false; + int lgth(key->Length(info)); + const double *data2(info->Get(key)); + data.insert(data.end(),data2,data2+lgth); + return true; +} + +void ExtractInfo(vtkInformationVector *inputVector, vtkUnstructuredGrid *& usgIn) +{ + vtkInformation *inputInfo(inputVector->GetInformationObject(0)); + vtkDataSet *input(0); + vtkDataSet *input0(vtkDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT()))); + vtkMultiBlockDataSet *input1(vtkMultiBlockDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT()))); + if(input0) + input=input0; + else + { + if(!input1) + throw INTERP_KERNEL::Exception("Input dataSet must be a DataSet or single elt multi block dataset expected !"); + if(input1->GetNumberOfBlocks()!=1) + throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with not exactly one block ! Use MergeBlocks or ExtractBlocks filter before calling this filter !"); + vtkDataObject *input2(input1->GetBlock(0)); + if(!input2) + throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with exactly one block but this single element is NULL !"); + vtkDataSet *input2c(vtkDataSet::SafeDownCast(input2)); + if(!input2c) + throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with exactly one block but this single element is not a dataset ! Use MergeBlocks or ExtractBlocks filter before calling this filter !"); + input=input2c; + } + if(!input) + throw INTERP_KERNEL::Exception("Input data set is NULL !"); + usgIn=vtkUnstructuredGrid::SafeDownCast(input); + if(!usgIn) + throw INTERP_KERNEL::Exception("Input data set is not an unstructured mesh ! This filter works only on unstructured meshes !"); +} + +DataArrayInt *ConvertVTKArrayToMCArrayInt(vtkDataArray *data) +{ + if(!data) + throw INTERP_KERNEL::Exception("ConvertVTKArrayToMCArrayInt : internal error !"); + int nbTuples(data->GetNumberOfTuples()),nbComp(data->GetNumberOfComponents()); + std::size_t nbElts(nbTuples*nbComp); + MCAuto ret(DataArrayInt::New()); + ret->alloc(nbTuples,nbComp); + for(int i=0;iGetComponentName(i)); + if(comp) + ret->setInfoOnComponent(i,comp); + } + int *ptOut(ret->getPointer()); + vtkIntArray *d0(vtkIntArray::SafeDownCast(data)); + if(d0) + { + const int *pt(d0->GetPointer(0)); + std::copy(pt,pt+nbElts,ptOut); + return ret.retn(); + } + vtkLongArray *d1(vtkLongArray::SafeDownCast(data)); + if(d1) + { + const long *pt(d1->GetPointer(0)); + std::copy(pt,pt+nbElts,ptOut); + return ret.retn(); + } + vtkIdTypeArray *d2(vtkIdTypeArray::SafeDownCast(data)); + if(d2) + { + const int *pt(d2->GetPointer(0)); + std::copy(pt,pt+nbElts,ptOut); + return ret.retn(); + } + std::ostringstream oss; + oss << "ConvertVTKArrayToMCArrayInt : unrecognized array \"" << typeid(*data).name() << "\" type !"; + throw INTERP_KERNEL::Exception(oss.str()); +} + +DataArrayDouble *ConvertVTKArrayToMCArrayDouble(vtkDataArray *data) +{ + if(!data) + throw INTERP_KERNEL::Exception("ConvertVTKArrayToMCArrayDouble : internal error !"); + int nbTuples(data->GetNumberOfTuples()),nbComp(data->GetNumberOfComponents()); + std::size_t nbElts(nbTuples*nbComp); + MCAuto ret(DataArrayDouble::New()); + ret->alloc(nbTuples,nbComp); + for(int i=0;iGetComponentName(i)); + if(comp) + ret->setInfoOnComponent(i,comp); + } + double *ptOut(ret->getPointer()); + vtkFloatArray *d0(vtkFloatArray::SafeDownCast(data)); + if(d0) + { + const float *pt(d0->GetPointer(0)); + for(std::size_t i=0;iGetPointer(0)); + std::copy(pt,pt+nbElts,ptOut); + return ret.retn(); + } + std::ostringstream oss; + oss << "ConvertVTKArrayToMCArrayDouble : unrecognized array \"" << typeid(*data).name() << "\" type !"; + throw INTERP_KERNEL::Exception(oss.str()); +} + +DataArray *ConvertVTKArrayToMCArray(vtkDataArray *data) +{ + if(!data) + throw INTERP_KERNEL::Exception("ConvertVTKArrayToMCArray : internal error !"); + vtkFloatArray *d0(vtkFloatArray::SafeDownCast(data)); + vtkDoubleArray *d1(vtkDoubleArray::SafeDownCast(data)); + if(d0 || d1) + return ConvertVTKArrayToMCArrayDouble(data); + vtkIntArray *d2(vtkIntArray::SafeDownCast(data)); + vtkLongArray *d3(vtkLongArray::SafeDownCast(data)); + if(d2 || d3) + return ConvertVTKArrayToMCArrayInt(data); + std::ostringstream oss; + oss << "ConvertVTKArrayToMCArray : unrecognized array \"" << typeid(*data).name() << "\" type !"; + throw INTERP_KERNEL::Exception(oss.str()); +} + +DataArrayDouble *BuildCoordsFrom(vtkPointSet *ds) +{ + if(!ds) + throw INTERP_KERNEL::Exception("BuildCoordsFrom : internal error !"); + vtkPoints *pts(ds->GetPoints()); + if(!pts) + throw INTERP_KERNEL::Exception("BuildCoordsFrom : internal error 2 !"); + vtkDataArray *data(pts->GetData()); + if(!data) + throw INTERP_KERNEL::Exception("BuildCoordsFrom : internal error 3 !"); + MCAuto coords(ConvertVTKArrayToMCArrayDouble(data)); + return coords.retn(); +} + +void ConvertFromUnstructuredGrid(vtkUnstructuredGrid *ds, std::vector< MCAuto >& ms, std::vector< MCAuto >& ids) +{ + MCAuto coords(BuildCoordsFrom(ds)); + vtkIdType nbCells(ds->GetNumberOfCells()); + vtkCellArray *ca(ds->GetCells()); + if(!ca) + return ; + vtkIdType nbEnt(ca->GetNumberOfConnectivityEntries()); + vtkIdType *caPtr(ca->GetPointer()); + vtkUnsignedCharArray *ct(ds->GetCellTypesArray()); + if(!ct) + throw INTERP_KERNEL::Exception("ConvertFromUnstructuredGrid : internal error"); + vtkIdTypeArray *cla(ds->GetCellLocationsArray()); + const vtkIdType *claPtr(cla->GetPointer(0)); + if(!cla) + throw INTERP_KERNEL::Exception("ConvertFromUnstructuredGrid : internal error 2"); + const unsigned char *ctPtr(ct->GetPointer(0)); + std::map m(ComputeMapOfType()); + MCAuto lev(DataArrayInt::New()) ; lev->alloc(nbCells,1); + int *levPtr(lev->getPointer()); + for(vtkIdType i=0;i::iterator it(m.find(ctPtr[i])); + if(it!=m.end()) + { + const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)(*it).second)); + levPtr[i]=cm.getDimension(); + } + else + { + std::ostringstream oss; oss << "ConvertFromUnstructuredGrid : at pos #" << i << " unrecognized VTK cell with type =" << ctPtr[i]; + throw INTERP_KERNEL::Exception(oss.str()); + } + } + MCAuto levs(lev->getDifferentValues()); + vtkIdTypeArray *faces(ds->GetFaces()),*faceLoc(ds->GetFaceLocations()); + for(const int *curLev=levs->begin();curLev!=levs->end();curLev++) + { + MCAuto m0(MEDCouplingUMesh::New("",*curLev)); + m0->setCoords(coords); m0->allocateCells(); + MCAuto cellIdsCurLev(lev->findIdsEqual(*curLev)); + for(const int *cellId=cellIdsCurLev->begin();cellId!=cellIdsCurLev->end();cellId++) + { + std::map::iterator it(m.find(ctPtr[*cellId])); + vtkIdType offset(claPtr[*cellId]); + vtkIdType sz(caPtr[offset]); + INTERP_KERNEL::NormalizedCellType ct((INTERP_KERNEL::NormalizedCellType)(*it).second); + if(ct!=INTERP_KERNEL::NORM_POLYHED) + { + std::vector conn2(sz); + for(int kk=0;kkinsertNextCell(ct,sz,&conn2[0]); + } + else + { + if(!faces || !faceLoc) + throw INTERP_KERNEL::Exception("ConvertFromUnstructuredGrid : faces are expected when there are polyhedra !"); + const vtkIdType *facPtr(faces->GetPointer(0)),*facLocPtr(faceLoc->GetPointer(0)); + std::vector conn; + int off0(facLocPtr[*cellId]); + int nbOfFaces(facPtr[off0++]); + for(int k=0;kinsertNextCell(ct,conn.size(),&conn[0]); + } + } + ms.push_back(m0); ids.push_back(cellIdsCurLev); + } +} + +vtkSmartPointer ConvertUMeshFromMCToVTK(const MEDCouplingUMesh *mVor) +{ + std::map zeMapRev(ComputeRevMapOfType()); + int nbCells(mVor->getNumberOfCells()); + vtkSmartPointer ret(vtkSmartPointer::New()); + ret->Initialize(); + ret->Allocate(); + vtkSmartPointer points(vtkSmartPointer::New()); + { + const DataArrayDouble *vorCoords(mVor->getCoords()); + vtkSmartPointer da(vtkSmartPointer::New()); + da->SetNumberOfComponents(vorCoords->getNumberOfComponents()); + da->SetNumberOfTuples(vorCoords->getNumberOfTuples()); + std::copy(vorCoords->begin(),vorCoords->end(),da->GetPointer(0)); + points->SetData(da); + } + mVor->checkConsistencyLight(); + switch(mVor->getMeshDimension()) + { + case 3: + { + int *cPtr(nullptr),*dPtr(nullptr); + unsigned char *aPtr(nullptr); + vtkSmartPointer cellTypes(vtkSmartPointer::New()); + { + cellTypes->SetNumberOfComponents(1); + cellTypes->SetNumberOfTuples(nbCells); + aPtr=cellTypes->GetPointer(0); + } + vtkSmartPointer cellLocations(vtkSmartPointer::New()); + { + cellLocations->SetNumberOfComponents(1); + cellLocations->SetNumberOfTuples(nbCells); + cPtr=cellLocations->GetPointer(0); + } + vtkSmartPointer cells(vtkSmartPointer::New()); + { + MCAuto tmp2(mVor->computeEffectiveNbOfNodesPerCell()); + cells->SetNumberOfComponents(1); + cells->SetNumberOfTuples(tmp2->accumulate(0)+nbCells); + dPtr=cells->GetPointer(0); + } + const int *connPtr(mVor->getNodalConnectivity()->begin()),*connIPtr(mVor->getNodalConnectivityIndex()->begin()); + int k(0),kk(0); + std::vector ee,ff; + for(int i=0;i(connPtr[connIPtr[0]])); + *aPtr++=zeMapRev[connPtr[connIPtr[0]]]; + if(ct!=INTERP_KERNEL::NORM_POLYHED) + { + int sz(connIPtr[1]-connIPtr[0]-1); + *dPtr++=sz; + dPtr=std::copy(connPtr+connIPtr[0]+1,connPtr+connIPtr[1],dPtr); + *cPtr++=k; k+=sz+1; + ee.push_back(kk); + } + else + { + std::set s(connPtr+connIPtr[0]+1,connPtr+connIPtr[1]); s.erase(-1); + int nbFace(std::count(connPtr+connIPtr[0]+1,connPtr+connIPtr[1],-1)+1); + ff.push_back(nbFace); + const int *work(connPtr+connIPtr[0]+1); + for(int j=0;j faceLocations(vtkSmartPointer::New()); + { + faceLocations->SetNumberOfComponents(1); + faceLocations->SetNumberOfTuples(ee.size()); + std::copy(ee.begin(),ee.end(),faceLocations->GetPointer(0)); + } + vtkSmartPointer faces(vtkSmartPointer::New()); + { + faces->SetNumberOfComponents(1); + faces->SetNumberOfTuples(ff.size()); + std::copy(ff.begin(),ff.end(),faces->GetPointer(0)); + } + vtkSmartPointer cells2(vtkSmartPointer::New()); + cells2->SetCells(nbCells,cells); + ret->SetCells(cellTypes,cellLocations,cells2,faceLocations,faces); + break; + } + case 2: + { + vtkSmartPointer cellTypes(vtkSmartPointer::New()); + { + cellTypes->SetNumberOfComponents(1); + cellTypes->SetNumberOfTuples(nbCells); + unsigned char *ptr(cellTypes->GetPointer(0)); + std::fill(ptr,ptr+nbCells,zeMapRev[(int)INTERP_KERNEL::NORM_POLYGON]); + } + int *cPtr(0),*dPtr(0); + vtkSmartPointer cellLocations(vtkSmartPointer::New()); + { + cellLocations->SetNumberOfComponents(1); + cellLocations->SetNumberOfTuples(nbCells); + cPtr=cellLocations->GetPointer(0); + } + vtkSmartPointer cells(vtkSmartPointer::New()); + { + cells->SetNumberOfComponents(1); + cells->SetNumberOfTuples(mVor->getNodalConnectivity()->getNumberOfTuples()); + dPtr=cells->GetPointer(0); + } + const int *connPtr(mVor->getNodalConnectivity()->begin()),*connIPtr(mVor->getNodalConnectivityIndex()->begin()); + int k(0); + for(int i=0;i cells2(vtkSmartPointer::New()); + cells2->SetCells(nbCells,cells); + ret->SetCells(cellTypes,cellLocations,cells2); + break; + } + case 1: + { + vtkSmartPointer cellTypes(vtkSmartPointer::New()); + { + cellTypes->SetNumberOfComponents(1); + cellTypes->SetNumberOfTuples(nbCells); + unsigned char *ptr(cellTypes->GetPointer(0)); + std::fill(ptr,ptr+nbCells,zeMapRev[(int)INTERP_KERNEL::NORM_SEG2]); + } + int *cPtr(0),*dPtr(0); + vtkSmartPointer cellLocations(vtkSmartPointer::New()); + { + cellLocations->SetNumberOfComponents(1); + cellLocations->SetNumberOfTuples(nbCells); + cPtr=cellLocations->GetPointer(0); + } + vtkSmartPointer cells(vtkSmartPointer::New()); + { + cells->SetNumberOfComponents(1); + cells->SetNumberOfTuples(mVor->getNodalConnectivity()->getNumberOfTuples()); + dPtr=cells->GetPointer(0); + } + const int *connPtr(mVor->getNodalConnectivity()->begin()),*connIPtr(mVor->getNodalConnectivityIndex()->begin()); + for(int i=0;i cells2(vtkSmartPointer::New()); + cells2->SetCells(nbCells,cells); + ret->SetCells(cellTypes,cellLocations,cells2); + break; + } + default: + throw INTERP_KERNEL::Exception("Not implemented yet !"); + } + ret->SetPoints(points); + return ret; +} + +class OffsetKeeper +{ +public: + OffsetKeeper():_vtk_arr(0) { } + void pushBack(vtkDataArray *da) { _da_on.push_back(da); } + void setVTKArray(vtkIdTypeArray *arr) { + MCAuto offmc(ConvertVTKArrayToMCArrayInt(arr)); + _off_arr=offmc; _vtk_arr=arr; } + const std::vector& getArrayGauss() const { return _da_on; } + const DataArrayInt *getOffsets() const { return _off_arr; } + vtkIdTypeArray *getVTKOffsets() const { return _vtk_arr; } +private: + std::vector _da_on; + MCAuto _off_arr; + vtkIdTypeArray *_vtk_arr; +}; + +void FillAdvInfoFrom(int vtkCT, const std::vector& GaussAdvData, int nbGaussPt, int nbNodesPerCell, std::vector& refCoo,std::vector& posInRefCoo) +{ + int nbOfCTS((int)GaussAdvData[0]),pos(1); + for(int i=0;i(std::cerr," ")); std::cerr << std::endl; + //std::copy(posInRefCoo.begin(),posInRefCoo.end(),std::ostream_iterator(std::cerr," ")); std::cerr << std::endl; + return ; + } + std::ostringstream oss; oss << "FillAdvInfoFrom : Internal error ! Not found cell type " << vtkCT << " in advanced Gauss info !"; + throw INTERP_KERNEL::Exception(oss.str()); +} + +template +vtkSmartPointer ExtractFieldFieldArr(T *elt2, int sizeOfOutArr, int nbOfCellsOfInput, const int *offsetsPtr, const int *nbPtsPerCellPtr) +{ + vtkSmartPointer elt3(vtkSmartPointer::New()); + int nbc(elt2->GetNumberOfComponents()); + elt3->SetNumberOfComponents(nbc); + elt3->SetNumberOfTuples(sizeOfOutArr); + for(int i=0;iGetComponentName(i)); + if(name) + elt3->SetComponentName(i,name); + } + elt3->SetName(elt2->GetName()); + // + U *ptr(elt3->GetPointer(0)); + const U *srcPtr(elt2->GetPointer(0)); + for(int i=0;i +vtkSmartPointer ExtractCellFieldArr(T *elt2, int sizeOfOutArr, int nbOfCellsOfInput, const int *idsPtr, const int *nbPtsPerCellPtr) +{ + vtkSmartPointer elt3(vtkSmartPointer::New()); + int nbc(elt2->GetNumberOfComponents()); + elt3->SetNumberOfComponents(nbc); + elt3->SetNumberOfTuples(sizeOfOutArr); + for(int i=0;iGetComponentName(i)); + if(name) + elt3->SetComponentName(i,name); + } + elt3->SetName(elt2->GetName()); + // + U *ptr(elt3->GetPointer(0)); + const U *srcPtr(elt2->GetPointer(0)); + for(int i=0;i Voronize(const MEDCouplingUMesh *m, const DataArrayInt *ids, vtkIdTypeArray *vtkOff, const DataArrayInt *offsetsBase, const std::vector& arrGauss, const std::vector& GaussAdvData, const std::vector& arrsOnCells) +{ + if(arrGauss.empty()) + throw INTERP_KERNEL::Exception("Voronize : no Gauss array !"); + int nbTuples(arrGauss[0]->GetNumberOfTuples()); + for(std::vector::const_iterator it=arrGauss.begin();it!=arrGauss.end();it++) + { + if((*it)->GetNumberOfTuples()!=nbTuples) + { + std::ostringstream oss; oss << "Mismatch of number of tuples in Gauss arrays for array \"" << (*it)->GetName() << "\""; + throw INTERP_KERNEL::Exception(oss.str()); + } + } + // Look at vtkOff has in the stomac + vtkInformation *info(vtkOff->GetInformation()); + if(!info) + throw INTERP_KERNEL::Exception("info is null ! Internal error ! Looks bad !"); + vtkInformationQuadratureSchemeDefinitionVectorKey *key(vtkQuadratureSchemeDefinition::DICTIONARY()); + if(!key->Has(info)) + throw INTERP_KERNEL::Exception("No quadrature key in info included in offets array ! Internal error ! Looks bad !"); + int dictSize(key->Size(info)); + INTERP_KERNEL::AutoPtr dict(new vtkQuadratureSchemeDefinition *[dictSize]); + key->GetRange(info,dict,0,0,dictSize); + // Voronoize + MCAuto field(MEDCouplingFieldDouble::New(ON_GAUSS_PT)); + field->setMesh(m); + // Gauss Part + int nbOfCellsOfInput(m->getNumberOfCells()); + MCAuto nbPtsPerCellArr(DataArrayInt::New()); nbPtsPerCellArr->alloc(nbOfCellsOfInput,1); + std::map zeMapRev(ComputeRevMapOfType()),zeMap(ComputeMapOfType()); + std::set agt(m->getAllGeoTypes()); + for(std::set::const_iterator it=agt.begin();it!=agt.end();it++) + { + const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(*it)); + std::map::const_iterator it2(zeMapRev.find((int)*it)); + if(it2==zeMapRev.end()) + throw INTERP_KERNEL::Exception("Internal error ! no type conversion available !"); + vtkQuadratureSchemeDefinition *gaussLoc(dict[(*it2).second]); + if(!gaussLoc) + { + std::ostringstream oss; oss << "For cell type " << cm.getRepr() << " no Gauss info !"; + throw INTERP_KERNEL::Exception(oss.str()); + } + int np(gaussLoc->GetNumberOfQuadraturePoints()),nbPtsPerCell((int)cm.getNumberOfNodes()); + const double *sfw(gaussLoc->GetShapeFunctionWeights()),*w(gaussLoc->GetQuadratureWeights());; + std::vector refCoo,posInRefCoo,wCpp(w,w+np); + FillAdvInfoFrom((*it2).second,GaussAdvData,np,nbPtsPerCell,refCoo,posInRefCoo); + field->setGaussLocalizationOnType(*it,refCoo,posInRefCoo,wCpp); + MCAuto ids2(m->giveCellsWithType(*it)); + nbPtsPerCellArr->setPartOfValuesSimple3(np,ids2->begin(),ids2->end(),0,1,1); + } + int zeSizeOfOutCellArr(nbPtsPerCellArr->accumulate(0)); + { MCAuto fakeArray(DataArrayDouble::New()); fakeArray->alloc(zeSizeOfOutCellArr,1); field->setArray(fakeArray); } + field->checkConsistencyLight(); + MCAuto vor(field->voronoize(1e-12));// The key is here ! + MEDCouplingUMesh *mVor(dynamic_cast(vor->getMesh())); + // + vtkSmartPointer ret(ConvertUMeshFromMCToVTK(mVor)); + // now fields... + MCAuto myOffsets(offsetsBase->selectByTupleIdSafe(ids->begin(),ids->end())); + const int *myOffsetsPtr(myOffsets->begin()),*nbPtsPerCellArrPtr(nbPtsPerCellArr->begin()); + for(std::vector::const_iterator it=arrGauss.begin();it!=arrGauss.end();it++) + { + vtkDataArray *elt(*it); + vtkDoubleArray *elt2(vtkDoubleArray::SafeDownCast(elt)); + vtkIntArray *elt3(vtkIntArray::SafeDownCast(elt)); + vtkIdTypeArray *elt4(vtkIdTypeArray::SafeDownCast(elt)); + if(elt2) + { + vtkSmartPointer arr(ExtractFieldFieldArr(elt2,zeSizeOfOutCellArr,nbOfCellsOfInput,myOffsetsPtr,nbPtsPerCellArrPtr)); + ret->GetCellData()->AddArray(arr); + continue; + } + if(elt3) + { + vtkSmartPointer arr(ExtractFieldFieldArr(elt3,zeSizeOfOutCellArr,nbOfCellsOfInput,myOffsetsPtr,nbPtsPerCellArrPtr)); + ret->GetCellData()->AddArray(arr); + continue; + } + if(elt4) + { + vtkSmartPointer arr(ExtractFieldFieldArr(elt4,zeSizeOfOutCellArr,nbOfCellsOfInput,myOffsetsPtr,nbPtsPerCellArrPtr)); + ret->GetCellData()->AddArray(arr); + continue; + } + } + for(std::vector::const_iterator it=arrsOnCells.begin();it!=arrsOnCells.end();it++) + { + vtkDataArray *elt(*it); + vtkDoubleArray *elt2(vtkDoubleArray::SafeDownCast(elt)); + vtkIntArray *elt3(vtkIntArray::SafeDownCast(elt)); + vtkIdTypeArray *elt4(vtkIdTypeArray::SafeDownCast(elt)); + if(elt2) + { + vtkSmartPointer arr(ExtractCellFieldArr(elt2,zeSizeOfOutCellArr,nbOfCellsOfInput,ids->begin(),nbPtsPerCellArrPtr)); + ret->GetCellData()->AddArray(arr); + continue; + } + if(elt3) + { + vtkSmartPointer arr(ExtractCellFieldArr(elt3,zeSizeOfOutCellArr,nbOfCellsOfInput,ids->begin(),nbPtsPerCellArrPtr)); + ret->GetCellData()->AddArray(arr); + continue; + } + if(elt4) + { + vtkSmartPointer arr(ExtractCellFieldArr(elt4,zeSizeOfOutCellArr,nbOfCellsOfInput,ids->begin(),nbPtsPerCellArrPtr)); + ret->GetCellData()->AddArray(arr); + continue; + } + } + return ret; +} + +vtkSmartPointer ComputeVoroGauss(vtkUnstructuredGrid *usgIn, const std::vector& GaussAdvData) +{ + OffsetKeeper zeOffsets; + std::string zeArrOffset; + std::vector cellFieldNamesToDestroy; + { + int nArrays(usgIn->GetFieldData()->GetNumberOfArrays()); + for(int i=0;iGetFieldData()->GetArray(i)); + if(!array) + continue; + const char* arrayOffsetName(array->GetInformation()->Get(vtkQuadratureSchemeDefinition::QUADRATURE_OFFSET_ARRAY_NAME())); + if(!arrayOffsetName) + continue; + std::string arrOffsetNameCpp(arrayOffsetName); + cellFieldNamesToDestroy.push_back(arrOffsetNameCpp); + if(arrOffsetNameCpp.find("ELGA@")==std::string::npos) + continue; + if(zeArrOffset.empty()) + zeArrOffset=arrOffsetNameCpp; + else + if(zeArrOffset!=arrOffsetNameCpp) + { + throw INTERP_KERNEL::Exception("ComputeVoroGauss : error in QUADRATURE_OFFSET_ARRAY_NAME for Gauss fields array !"); + } + zeOffsets.pushBack(array); + } + if(zeArrOffset.empty()) + throw INTERP_KERNEL::Exception("ComputeVoroGauss : no Gauss points fields in DataSet !"); + } + std::vector< MCAuto > ms; + std::vector< MCAuto > ids; + ConvertFromUnstructuredGrid(usgIn,ms,ids); + { + vtkDataArray *offTmp(usgIn->GetCellData()->GetArray(zeArrOffset.c_str())); + if(!offTmp) + { + std::ostringstream oss; oss << "ComputeVoroGauss : cell field " << zeArrOffset << " not found !"; + throw INTERP_KERNEL::Exception(oss.str()); + } + vtkIdTypeArray *offsets(vtkIdTypeArray::SafeDownCast(offTmp)); + if(!offsets) + { + std::ostringstream oss; oss << "ComputeVoroGauss : cell field " << zeArrOffset << " exists but not with the right type of data !"; + throw INTERP_KERNEL::Exception(oss.str()); + } + /// + zeOffsets.setVTKArray(offsets); + } + // + std::vector arrsOnCells; + { + int nArrays(usgIn->GetCellData()->GetNumberOfArrays()); + for(int i=0;iGetCellData()->GetArray(i)); + if(!array) + continue; + std::string name(array->GetName()); + if(std::find(cellFieldNamesToDestroy.begin(),cellFieldNamesToDestroy.end(),name)==cellFieldNamesToDestroy.end()) + { + arrsOnCells.push_back(array); + } + } + { + vtkDataArray *doNotKeepThis(zeOffsets.getVTKOffsets()); + std::vector::iterator it2(std::find(arrsOnCells.begin(),arrsOnCells.end(),doNotKeepThis)); + if(it2!=arrsOnCells.end()) + arrsOnCells.erase(it2); + } + } + // + std::size_t sz(ms.size()); + std::vector< vtkSmartPointer > res; + for(std::size_t i=0;i mmc(ms[i]); + MCAuto myIds(ids[i]); + vtkSmartPointer vor(Voronize(mmc,myIds,zeOffsets.getVTKOffsets(),zeOffsets.getOffsets(),zeOffsets.getArrayGauss(),GaussAdvData,arrsOnCells)); + res.push_back(vor); + } + if(res.empty()) + throw INTERP_KERNEL::Exception("Dataset is empty !"); + vtkSmartPointer mb(vtkSmartPointer::New()); + vtkSmartPointer cd(vtkSmartPointer::New()); + for(std::vector< vtkSmartPointer >::const_iterator it=res.begin();it!=res.end();it++) + mb->AddInputData(*it); + cd->SetInputConnection(mb->GetOutputPort()); + cd->SetMergePoints(0); + cd->Update(); + vtkSmartPointer ret; + ret=cd->GetOutput(); + return ret; +} + +//////////////////// + +vtkVoroGauss::vtkVoroGauss() +{ + this->SetNumberOfInputPorts(1); + this->SetNumberOfOutputPorts(1); +} + +vtkVoroGauss::~vtkVoroGauss() +{ +} + +int vtkVoroGauss::RequestInformation(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) +{ + //std::cerr << "########################################## vtkVoroGauss::RequestInformation ##########################################" << std::endl; + try + { + vtkUnstructuredGrid *usgIn(0); + ExtractInfo(inputVector[0],usgIn); + } + catch(INTERP_KERNEL::Exception& e) + { + std::ostringstream oss; + oss << "Exception has been thrown in vtkVoroGauss::RequestInformation : " << e.what() << std::endl; + if(this->HasObserver("ErrorEvent") ) + this->InvokeEvent("ErrorEvent",const_cast(oss.str().c_str())); + else + vtkOutputWindowDisplayErrorText(const_cast(oss.str().c_str())); + vtkObject::BreakOnError(); + return 0; + } + return 1; +} + +int vtkVoroGauss::RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) +{ + //std::cerr << "########################################## vtkVoroGauss::RequestData ##########################################" << std::endl; + try + { + + std::vector GaussAdvData; + bool isOK(IsInformationOK(inputVector[0]->GetInformationObject(0),GaussAdvData)); + if(!isOK) + throw INTERP_KERNEL::Exception("Sorry but no advanced gauss info found ! Expect to be called right after a MEDReader containing Gauss Points !"); + vtkUnstructuredGrid *usgIn(0); + ExtractInfo(inputVector[0],usgIn); + // + vtkSmartPointer ret(ComputeVoroGauss(usgIn,GaussAdvData)); + vtkInformation *inInfo(inputVector[0]->GetInformationObject(0)); + vtkInformation *outInfo(outputVector->GetInformationObject(0)); + vtkUnstructuredGrid *output(vtkUnstructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()))); + output->ShallowCopy(ret); + } + catch(INTERP_KERNEL::Exception& e) + { + std::ostringstream oss; + oss << "Exception has been thrown in vtkVoroGauss::RequestData : " << e.what() << std::endl; + if(this->HasObserver("ErrorEvent") ) + this->InvokeEvent("ErrorEvent",const_cast(oss.str().c_str())); + else + vtkOutputWindowDisplayErrorText(const_cast(oss.str().c_str())); + vtkObject::BreakOnError(); + return 0; + } + return 1; +} + +void vtkVoroGauss::PrintSelf(ostream& os, vtkIndent indent) +{ + this->Superclass::PrintSelf(os, indent); +} diff --git a/src/Plugins/VoroGauss/IO/vtkVoroGauss.h b/src/Plugins/VoroGauss/IO/vtkVoroGauss.h new file mode 100644 index 00000000..d344cb1d --- /dev/null +++ b/src/Plugins/VoroGauss/IO/vtkVoroGauss.h @@ -0,0 +1,52 @@ +// Copyright (C) 2017 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 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 +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay (EDF R&D) + +#ifndef vtkVoroGauss_h__ +#define vtkVoroGauss_h__ + +#include "vtkUnstructuredGridAlgorithm.h" + +class vtkMutableDirectedGraph; + +class VTK_EXPORT vtkVoroGauss : public vtkUnstructuredGridAlgorithm +{ +public: + static vtkVoroGauss* New(); + vtkTypeMacro(vtkVoroGauss, vtkUnstructuredGridAlgorithm) + void PrintSelf(ostream& os, vtkIndent indent); + +protected: + vtkVoroGauss(); + ~vtkVoroGauss(); + + int RequestInformation(vtkInformation *request, + vtkInformationVector **inputVector, vtkInformationVector *outputVector); + + int RequestData(vtkInformation *request, vtkInformationVector **inputVector, + vtkInformationVector *outputVector); +private: + vtkVoroGauss(const vtkVoroGauss&); + void operator=(const vtkVoroGauss&); // Not implemented. + private: + //BTX + //ETX +}; + +#endif diff --git a/src/Plugins/VoroGauss/PG_3D.med b/src/Plugins/VoroGauss/PG_3D.med new file mode 100644 index 00000000..9c7e51e7 Binary files /dev/null and b/src/Plugins/VoroGauss/PG_3D.med differ diff --git a/src/Plugins/VoroGauss/ParaViewPlugin/CMakeLists.txt b/src/Plugins/VoroGauss/ParaViewPlugin/CMakeLists.txt new file mode 100644 index 00000000..1e5cfe5d --- /dev/null +++ b/src/Plugins/VoroGauss/ParaViewPlugin/CMakeLists.txt @@ -0,0 +1,26 @@ +# Copyright (C) 2016 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +INCLUDE_DIRECTORIES( ${CMAKE_CURRENT_SOURCE_DIR}/../IO ) +ADD_PARAVIEW_PLUGIN(VoroGaussPlugin "1.0" + SERVER_MANAGER_SOURCES ${SM_SRCS} ${PROJECT_SOURCE_DIR}/IO/vtkVoroGauss.cxx + SERVER_MANAGER_XML Resources/VoroGaussServer.xml) +TARGET_LINK_LIBRARIES(VoroGaussPlugin ${MEDCoupling_medcoupling}) +INSTALL(TARGETS VoroGaussPlugin RUNTIME DESTINATION lib/paraview LIBRARY DESTINATION lib/paraview ARCHIVE DESTINATION lib/paraview) diff --git a/src/Plugins/VoroGauss/ParaViewPlugin/Resources/VoroGaussServer.xml b/src/Plugins/VoroGauss/ParaViewPlugin/Resources/VoroGaussServer.xml new file mode 100644 index 00000000..0a7b10b7 --- /dev/null +++ b/src/Plugins/VoroGauss/ParaViewPlugin/Resources/VoroGaussServer.xml @@ -0,0 +1,21 @@ + + + + + + + + + + + + + This property specifies the input to the Level Scalars filter. + + + + + + + + diff --git a/src/Plugins/VoroGauss/TestCase.py b/src/Plugins/VoroGauss/TestCase.py new file mode 100644 index 00000000..de5c5049 --- /dev/null +++ b/src/Plugins/VoroGauss/TestCase.py @@ -0,0 +1,74 @@ +# Copyright (C) 2017 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 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 +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +from MEDLoader import * + +fname="VoroGauss1.med" +meshName="mesh" +mm=MEDFileUMesh() +coords=DataArrayDouble([0,0, 1,0, 2,0, 3,0, 4,0, 5,0, 0,1, 1,1, 2,1, 0,2, 1,2, 3,1, 4,1],13,2) +m0=MEDCouplingUMesh(meshName,2) +m0.setCoords(coords) +m0.allocateCells() +m0.insertNextCell(NORM_TRI3,[2,3,8]) +m0.insertNextCell(NORM_TRI3,[3,4,11]) +m0.insertNextCell(NORM_TRI3,[4,5,12]) +m0.insertNextCell(NORM_TRI3,[6,7,9]) +m0.insertNextCell(NORM_TRI3,[7,8,10]) +m0.insertNextCell(NORM_QUAD4,[0,1,7,6]) +m0.insertNextCell(NORM_QUAD4,[1,2,8,7]) +mm[0]=m0 +m1=MEDCouplingUMesh(meshName,1) +m1.setCoords(coords) +m1.allocateCells() +m1.insertNextCell(NORM_SEG2,[0,1]) +m1.insertNextCell(NORM_SEG2,[1,2]) +m1.insertNextCell(NORM_SEG2,[2,3]) +m1.insertNextCell(NORM_SEG2,[3,4]) +m1.insertNextCell(NORM_SEG2,[4,5]) +mm[-1]=m1 +mm.setFamilyFieldArr(0,DataArrayInt([-1,-1,-2,-3,-3,-1,-3])) +mm.setFamilyFieldArr(-1,DataArrayInt([-1,-4,-4,-4,-1])) +for i in [-1,-2,-3,-4]: + mm.setFamilyId("Fam_%d"%i,i) + mm.setFamiliesOnGroup("G%d"%(abs(i)),["Fam_%d"%i]) + pass +mm.write(fname,2) +# +f0=MEDCouplingFieldDouble(ON_GAUSS_PT) +f0.setMesh(m0) +f0.setName("MyFieldPG") ; f0.setMesh(m0) +f0.setGaussLocalizationOnType(NORM_TRI3,[0,0, 1,0, 0,1],[0.1,0.1, 0.8,0.1, 0.1,0.8],[0.3,0.3,0.4]) +f0.setGaussLocalizationOnType(NORM_QUAD4,[-1,-1, 1,-1, 1,1, -1,1],[-0.57735,-0.57735,0.57735,-0.57735,0.57735,0.57735,-0.57735,0.57735],[0.25,0.25,0.25,0.25]) +arr=DataArrayDouble(f0.getNumberOfTuplesExpected()) ; arr.iota() +arr=DataArrayDouble.Meld(arr,arr) +arr.setInfoOnComponents(["comp0","comp1"]) +f0.setArray(arr) +WriteFieldUsingAlreadyWrittenMesh(fname,f0) +# +f1=MEDCouplingFieldDouble(ON_CELLS) +f1.setMesh(m0) +f1.setName("MyFieldCell") ; f1.setMesh(m0) +arr=DataArrayDouble(f1.getNumberOfTuplesExpected()) ; arr.iota() +arr=DataArrayDouble.Meld(arr,arr) +arr.setInfoOnComponents(["comp2","comp3"]) +f1.setArray(arr) +WriteFieldUsingAlreadyWrittenMesh(fname,f1) + diff --git a/src/Plugins/VoroGauss/testMEDReader14.med b/src/Plugins/VoroGauss/testMEDReader14.med new file mode 100644 index 00000000..c7ff867f Binary files /dev/null and b/src/Plugins/VoroGauss/testMEDReader14.med differ