1 // Copyright (C) 2017-2020 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (EDF R&D)
21 #include "vtkDevelopedSurface.h"
22 #include "VTKToMEDMem.h"
24 #include "vtkAdjacentVertexIterator.h"
25 #include "vtkIntArray.h"
26 #include "vtkLongArray.h"
27 #include "vtkCellData.h"
28 #include "vtkPointData.h"
29 #include "vtkCylinder.h"
31 #include "vtkCutter.h"
32 #include "vtkTransform.h"
34 #include "vtkStreamingDemandDrivenPipeline.h"
35 #include "vtkUnstructuredGrid.h"
36 #include "vtkMultiBlockDataSet.h"
38 #include "vtkInformationStringKey.h"
39 #include "vtkAlgorithmOutput.h"
40 #include "vtkObjectFactory.h"
41 #include "vtkMutableDirectedGraph.h"
42 #include "vtkMultiBlockDataSet.h"
43 #include "vtkDataSet.h"
44 #include "vtkInformationVector.h"
45 #include "vtkInformation.h"
46 #include "vtkDataArraySelection.h"
47 #include "vtkTimeStamp.h"
48 #include "vtkInEdgeIterator.h"
49 #include "vtkInformationDataObjectKey.h"
50 #include "vtkExecutive.h"
51 #include "vtkVariantArray.h"
52 #include "vtkStringArray.h"
53 #include "vtkDoubleArray.h"
54 #include "vtkFloatArray.h"
55 #include "vtkCharArray.h"
56 #include "vtkUnsignedCharArray.h"
57 #include "vtkDataSetAttributes.h"
58 #include "vtkDemandDrivenPipeline.h"
59 #include "vtkDataObjectTreeIterator.h"
60 #include "vtkWarpScalar.h"
62 #include "MEDCouplingMemArray.hxx"
64 #include "VTKMEDTraits.hxx"
67 #define _USE_MATH_DEFINES
76 vtkStandardNewMacro(vtkDevelopedSurface);
86 struct VTKTraits<double>
88 typedef vtkDoubleArray ArrayType;
92 struct VTKTraits<float>
94 typedef vtkFloatArray ArrayType;
97 void ExtractInfo(vtkInformationVector *inputVector, vtkDataSet *& usgIn)
99 vtkInformation *inputInfo(inputVector->GetInformationObject(0));
100 vtkDataSet *input(0);
101 vtkDataSet *input0(vtkDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
102 vtkMultiBlockDataSet *input1(vtkMultiBlockDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
108 throw MZCException("Input dataSet must be a DataSet or single elt multi block dataset expected !");
109 if(input1->GetNumberOfBlocks()!=1)
110 throw MZCException("Input dataSet is a multiblock dataset with not exactly one block ! Use MergeBlocks or ExtractBlocks filter before calling this filter !");
111 vtkDataObject *input2(input1->GetBlock(0));
113 throw MZCException("Input dataSet is a multiblock dataset with exactly one block but this single element is NULL !");
114 vtkDataSet *input2c(vtkDataSet::SafeDownCast(input2));
116 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 !");
120 throw MZCException("Input data set is NULL !");
121 vtkPointData *att(input->GetPointData());
123 throw MZCException("Input dataset has no point data attribute ! Impossible to deduce a developed surface on it !");
127 class vtkDevelopedSurface::vtkInternals
130 vtkNew<vtkCutter> Cutter;
135 vtkDevelopedSurface::vtkDevelopedSurface():_cyl(nullptr),Internal(new vtkInternals),InvertStatus(false),OffsetInRad(0.)
137 //this->RegisterFilter(this->Internal->Cutter.GetPointer());
140 vtkDevelopedSurface::~vtkDevelopedSurface()
142 delete this->Internal;
145 void vtkDevelopedSurface::SetInvertWay(bool invertStatus)
147 this->InvertStatus=invertStatus;
151 void vtkDevelopedSurface::SetThetaOffset(double offsetInDegrees)
153 double tmp(std::min(offsetInDegrees,180.));
154 tmp=std::max(tmp,-180.);
155 this->OffsetInRad=tmp/180.*M_PI;
159 int vtkDevelopedSurface::RequestInformation(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
161 //std::cerr << "########################################## vtkDevelopedSurface::RequestInformation ##########################################" << std::endl;
164 vtkDataSet *usgIn(0);
165 ExtractInfo(inputVector[0],usgIn);
167 catch(MZCException& e)
169 std::ostringstream oss;
170 oss << "Exception has been thrown in vtkDevelopedSurface::RequestInformation : " << e.what() << std::endl;
171 if(this->HasObserver("ErrorEvent") )
172 this->InvokeEvent("ErrorEvent",const_cast<char *>(oss.str().c_str()));
174 vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
175 vtkObject::BreakOnError();
181 std::vector<mcIdType> UnWrapByDuplicatingNodes(vtkCellArray *ca, vtkIdType& offset, const MEDCoupling::DataArrayDouble *thetas)
183 std::vector<mcIdType> ret;
184 vtkIdType nbCells(ca->GetNumberOfCells());
185 vtkIdType *conn(ca->GetData()->GetPointer(0));
186 const double *tptr(thetas->begin());
187 for(vtkIdType i=0;i<nbCells;i++)
189 vtkIdType sz(*conn++);
190 double ma(-std::numeric_limits<double>::max()),mi(std::numeric_limits<double>::max());
191 for(vtkIdType j=0;j<sz;j++)
193 double tmp_theta(tptr[conn[j]]);
194 ma=std::max(ma,tmp_theta);
195 mi=std::min(mi,tmp_theta);
199 for(vtkIdType j=0;j<sz;j++)
201 double tmp_theta(tptr[conn[j]]);
204 ret.push_back(conn[j]);
215 void DealArray(vtkDataSetAttributes *pd, int pos, typename MEDFileVTKTraits<T>::VtkType *arr, std::vector<int>& nodeSel)
217 int nbc(arr->GetNumberOfComponents());
218 std::size_t nbt(nodeSel.size());
219 vtkSmartPointer< typename MEDFileVTKTraits<T>::VtkType > newArr;
220 newArr.TakeReference(MEDFileVTKTraits<T>::VtkType::New());
221 newArr->SetNumberOfComponents(nbc);
222 newArr->SetNumberOfTuples(nbt);
223 T *ptr(newArr->GetPointer(0));
224 const T *inPtr(arr->GetPointer(0));
225 for(std::size_t i=0;i<nbt;i++)
227 ptr=std::copy(inPtr+nbc*nodeSel[i],inPtr+nbc*(nodeSel[i]+1),ptr);
229 newArr->SetName(arr->GetName());
230 arr->DeepCopy(newArr);
233 void ToDouble(vtkDataArray *coords, vtkSmartPointer<vtkDoubleArray>& coordsOut)
235 vtkDoubleArray *coords2(vtkDoubleArray::SafeDownCast(coords));
236 vtkFloatArray *coords3(vtkFloatArray::SafeDownCast(coords));
237 if(!coords2 && !coords3)
238 throw MZCException("Input coordinates are neither float64 or float32 !");
242 coordsOut.TakeReference(coords2);
243 coords2->Register(0);
247 coordsOut.TakeReference(vtkDoubleArray::New());
248 coordsOut->SetNumberOfComponents(3);
249 vtkIdType nbTuples(coords3->GetNumberOfTuples());
250 coordsOut->SetNumberOfTuples(nbTuples);
251 std::copy(coords3->GetPointer(0),coords3->GetPointer(0)+3*nbTuples,coordsOut->GetPointer(0));
255 void dealWith(vtkPolyData *outdata, const double center[3], const double axis[3], double radius, double eps, bool invertThetaInc, double offsetInRad)
257 vtkDataArray *coords(outdata->GetPoints()->GetData());
258 if(coords->GetNumberOfComponents()!=3)
259 throw MZCException("Input coordinates are expected to have 3 components !");
261 vtkIdType nbNodes(coords->GetNumberOfTuples());
263 throw MZCException("No points -> impossible to develop anything !");
265 vtkSmartPointer<vtkDoubleArray> zeCoords;
266 ToDouble(coords,zeCoords);
268 double axis_cross_Z[3]={axis[1],-axis[0],0.};
269 double n_axis(sqrt(axis_cross_Z[0]*axis_cross_Z[0]+axis_cross_Z[1]*axis_cross_Z[1]));
272 double ang(asin(n_axis));
275 MEDCoupling::DataArrayDouble::Rotate3DAlg(center,axis_cross_Z,ang,nbNodes,zeCoords->GetPointer(0),zeCoords->GetPointer(0));
278 MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> c_cyl;
280 MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> cc(MEDCoupling::DataArrayDouble::New()); cc->alloc(nbNodes,3);
281 double *ccPtr(cc->getPointer());
282 const double *zeCoordsPtr(zeCoords->GetPointer(0));
283 for(vtkIdType i=0;i<nbNodes;i++,ccPtr+=3,zeCoordsPtr+=3)
285 std::transform(zeCoordsPtr,zeCoordsPtr+3,center,ccPtr,std::minus<double>());
287 c_cyl=cc->fromCartToCyl();
289 MEDCoupling::MCAuto<MEDCoupling::MEDFileData> mfd(MEDCoupling::MEDFileData::New());
290 WriteMEDFileFromVTKDataSet(mfd,outdata,{},0.,0);
293 MEDCoupling::MEDFileMeshes *ms(mfd->getMeshes());
294 if(ms->getNumberOfMeshes()!=1)
295 throw MZCException("Unexpected number of meshes !");
296 MEDCoupling::MEDFileMesh *mm(ms->getMeshAtPos(0));
297 MEDCoupling::MEDFileUMesh *mmu(dynamic_cast<MEDCoupling::MEDFileUMesh *>(mm));
299 throw MZCException("Expecting unstructured one !");
300 MEDCoupling::MCAuto<MEDCoupling::MEDCouplingUMesh> m0(mmu->getMeshAtLevel(0));
303 MEDCoupling::MCAuto<MEDCoupling::DataArrayIdType> c0s(m0->getCellIdsLyingOnNodes(&v,&v+1,false));
305 throw MZCException("Orphan node 0 !");
306 std::vector<mcIdType> nodes0;
307 m0->getNodeIdsOfCell(c0s->getIJ(0,0),nodes0);
308 MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> tmp0(c_cyl->selectByTupleIdSafe(nodes0.data(),nodes0.data()+nodes0.size()));
309 tmp0=tmp0->keepSelectedComponents({1});
310 double tmp(tmp0->getMaxAbsValueInArray());
315 constexpr double EPS_FOR_RADIUS=1e-2;
316 MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> rs(c_cyl->keepSelectedComponents({0}));
317 if(!rs->isUniform(radius,radius*EPS_FOR_RADIUS))
319 double mi(rs->getMinValueInArray()),ma(rs->getMaxValueInArray());
320 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 << " !";
321 throw MZCException(oss.str());
323 double tetha0(c_cyl->getIJ(0,1));
325 double *ccylptr(c_cyl->getPointer()+1);
326 double mi02(std::numeric_limits<double>::max());
327 for(vtkIdType i=0;i<nbNodes;i++,ccylptr+=3)
329 *ccylptr-=tetha0+offsetInRad;
330 mi02=std::min(mi02,*ccylptr);
334 ccylptr=c_cyl->getPointer()+1;
335 for(vtkIdType i=0;i<nbNodes;i++,ccylptr+=3)
345 MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> c_cyl_2(c_cyl->keepSelectedComponents({1}));
347 MEDCoupling::MCAuto<MEDCoupling::DataArrayIdType> poses(c_cyl_2->findIdsInRange(0.,eps));
348 c_cyl->setPartOfValuesSimple3(0.,poses->begin(),poses->end(),1,2,1);
351 if(a ^ (!invertThetaInc))
353 MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> tmp(c_cyl->keepSelectedComponents({1}));
355 std::for_each(tmp->getPointer(),tmp->getPointer()+tmp->getNumberOfTuples(),[](double& v) { if(v==-0.) v=0.; });
356 c_cyl->setPartOfValues1(tmp,0,nbNodes,1,1,2,1);
358 MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> c_cyl_post(c_cyl->keepSelectedComponents({1}));
360 double *c_cyl_post_ptr(c_cyl_post->getPointer());
361 for(vtkIdType i=0;i<nbNodes;i++,c_cyl_post_ptr++)
363 if(*c_cyl_post_ptr<0.)
364 *c_cyl_post_ptr+=2*M_PI;
367 vtkCellArray *cb(outdata->GetPolys());
368 vtkIdType offset(nbNodes);
369 std::vector<mcIdType> dupNodes(UnWrapByDuplicatingNodes(cb,offset,c_cyl_post));
371 MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> c_cyl_post2(c_cyl_post->selectByTupleId(dupNodes.data(),dupNodes.data()+dupNodes.size()));
372 c_cyl_post2->applyLin(1.,2*M_PI);
373 c_cyl_post=MEDCoupling::DataArrayDouble::Aggregate(c_cyl_post,c_cyl_post2);
374 MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> z0(c_cyl->keepSelectedComponents({2}));
375 MEDCoupling::MCAuto<MEDCoupling::DataArrayDouble> z1(z0->selectByTupleId(dupNodes.data(),dupNodes.data()+dupNodes.size()));
376 z0=MEDCoupling::DataArrayDouble::Aggregate(z0,z1);
378 std::size_t outNbNodes(z0->getNumberOfTuples());
379 vtkSmartPointer<vtkDoubleArray> zeCoords2;
380 zeCoords2.TakeReference(vtkDoubleArray::New());
381 zeCoords2->SetNumberOfComponents(3);
382 zeCoords2->SetNumberOfTuples(outNbNodes);
384 const double *tptr(c_cyl_post->begin()),*zptr(z0->begin());
385 double *outPtr(zeCoords2->GetPointer(0));
386 for(std::size_t i=0;i<outNbNodes;i++,tptr++,zptr++,outPtr+=3)
388 outPtr[0]=radius*(*tptr);
394 outdata->GetPoints()->SetData(zeCoords2);
395 // now post process nodes
396 std::vector<int> nodeSel(nbNodes+dupNodes.size());
399 std::for_each(nodeSel.begin(),nodeSel.begin()+nbNodes,[&cnt](int& v){ v=cnt++; });
400 std::copy(dupNodes.begin(),dupNodes.end(),nodeSel.begin()+nbNodes);
402 vtkDataSetAttributes *pd(outdata->GetPointData());
403 int nba(pd->GetNumberOfArrays());
404 for(int i=0;i<nba;i++)
406 vtkDataArray *arr(pd->GetArray(i));
408 vtkIntArray *arr0(vtkIntArray::SafeDownCast(arr));
411 DealArray<int>(pd,i,arr0,nodeSel);
416 vtkFloatArray *arr0(vtkFloatArray::SafeDownCast(arr));
419 DealArray<float>(pd,i,arr0,nodeSel);
424 vtkDoubleArray *arr0(vtkDoubleArray::SafeDownCast(arr));
427 DealArray<double>(pd,i,arr0,nodeSel);
434 int vtkDevelopedSurface::RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
436 //std::cerr << "########################################## vtkDevelopedSurface::RequestData ##########################################" << std::endl;
440 throw MZCException("No cylinder object as cut function !");
441 double center[3],axis[3],radius;
442 vtkAbstractTransform* trf(_cyl->GetTransform());
444 _cyl->GetCenter(center);
445 _cyl->GetAxis(axis[0],axis[1],axis[2]);
446 radius=_cyl->GetRadius();
450 double axis3[3]={center[0]+0.,center[1]+1.,center[2]+0.},axis4[3];
451 trf->TransformPoint(axis3,axis4);
452 std::transform(axis4,axis4+3,center,axis,[](double a, double b) { return b-a; });
454 if(std::isnan(axis[0]) && std::isnan(axis[1]) && std::isnan(axis[2]))
455 { axis[0]=0.; axis[1]=-1.; axis[2]=0.; }
457 //std::cerr << trf << " jjj " << axis[0] << " " << axis[1] << " " << axis[2] << " : " << center[0] << " " << center[1] << " " << center[2] << " " " " << " -> " << radius << std::endl;
458 vtkDataSet *usgIn(0);
459 ExtractInfo(inputVector[0],usgIn);
460 vtkSmartPointer<vtkPolyData> outData;
462 vtkNew<vtkCutter> Cutter;
463 Cutter->SetInputData(usgIn);
464 Cutter->SetCutFunction(_cyl);
466 vtkDataSet *zeComputedOutput(Cutter->GetOutput());
467 vtkPolyData *zeComputedOutput2(vtkPolyData::SafeDownCast(zeComputedOutput));
468 if(!zeComputedOutput2)
469 throw MZCException("Unexpected output of cutter !");
470 outData.TakeReference(zeComputedOutput2);
471 zeComputedOutput2->Register(0);
473 if(outData->GetNumberOfCells()==0)
474 return 1;// no cells -> nothing to do
476 dealWith(outData,center,axis,radius,1e-7,this->InvertStatus,this->OffsetInRad);
478 vtkInformation *outInfo(outputVector->GetInformationObject(0));
479 vtkPolyData *output(vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())));
480 output->ShallowCopy(outData);
482 catch(MZCException& e)
484 std::ostringstream oss;
485 oss << "Exception has been thrown in vtkDevelopedSurface::RequestInformation : " << e.what() << std::endl;
486 if(this->HasObserver("ErrorEvent") )
487 this->InvokeEvent("ErrorEvent",const_cast<char *>(oss.str().c_str()));
489 vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
490 vtkObject::BreakOnError();
496 int vtkDevelopedSurface::FillOutputPortInformation( int vtkNotUsed(port), vtkInformation* info)
498 info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData");
503 void vtkDevelopedSurface::PrintSelf(ostream& os, vtkIndent indent)
505 this->Superclass::PrintSelf(os, indent);
508 void vtkDevelopedSurface::SetCutFunction(vtkImplicitFunction* func)
510 vtkCylinder *cyl(vtkCylinder::SafeDownCast(func));
518 vtkMTimeType vtkDevelopedSurface::GetMTime()
520 vtkMTimeType maxMTime = this->Superclass::GetMTime(); // My MTime
523 maxMTime=std::max(maxMTime,_cyl->GetMTime());