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 "vtkVoroGauss.h"
23 #include "vtkAdjacentVertexIterator.h"
24 #include "vtkIntArray.h"
25 #include "vtkCellData.h"
26 #include "vtkPointData.h"
27 #include "vtkCellType.h"
29 #include "vtkCellArray.h"
30 #include "vtkIdTypeArray.h"
32 #include "vtkStreamingDemandDrivenPipeline.h"
33 #include "vtkUnstructuredGrid.h"
34 #include "vtkMultiBlockDataSet.h"
36 #include "vtkInformationStringKey.h"
37 #include "vtkAlgorithmOutput.h"
38 #include "vtkObjectFactory.h"
39 #include "vtkMutableDirectedGraph.h"
40 #include "vtkMultiBlockDataSet.h"
41 #include "vtkDataSet.h"
42 #include "vtkInformationVector.h"
43 #include "vtkInformation.h"
44 #include "vtkDataArraySelection.h"
45 #include "vtkTimeStamp.h"
46 #include "vtkInEdgeIterator.h"
47 #include "vtkInformationDataObjectKey.h"
48 #include "vtkInformationDataObjectMetaDataKey.h"
49 #include "vtkInformationDoubleVectorKey.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 "vtkLongArray.h"
57 #include "vtkUnsignedCharArray.h"
58 #include "vtkDataSetAttributes.h"
59 #include "vtkDemandDrivenPipeline.h"
60 #include "vtkDataObjectTreeIterator.h"
61 #include "vtkWarpScalar.h"
62 #include "vtkQuadratureSchemeDefinition.h"
63 #include "vtkInformationQuadratureSchemeDefinitionVectorKey.h"
64 #include "vtkCompositeDataToUnstructuredGridFilter.h"
65 #include "vtkMultiBlockDataGroupFilter.h"
67 #include "MEDCouplingMemArray.hxx"
68 #include "MEDCouplingMemArray.txx"
69 #include "MEDCouplingUMesh.hxx"
70 #include "MEDCouplingFieldDouble.hxx"
71 #include "InterpKernelAutoPtr.hxx"
72 #include "InterpKernelGaussCoords.hxx"
79 using MEDCoupling::DataArray;
80 using MEDCoupling::DataArrayInt32;
81 using MEDCoupling::DataArrayInt64;
82 using MEDCoupling::DataArrayDouble;
83 using MEDCoupling::MEDCouplingMesh;
84 using MEDCoupling::MEDCouplingUMesh;
85 using MEDCoupling::DynamicCastSafe;
86 using MEDCoupling::MEDCouplingFieldDouble;
87 using MEDCoupling::ON_GAUSS_PT;
88 using MEDCoupling::MCAuto;
90 vtkStandardNewMacro(vtkVoroGauss)
93 std::map<int,int> ComputeMapOfType()
95 std::map<int,int> ret;
96 int nbOfTypesInMC(sizeof(MEDCOUPLING2VTKTYPETRADUCER)/sizeof( decltype(MEDCOUPLING2VTKTYPETRADUCER[0]) ));
97 for(int i=0;i<nbOfTypesInMC;i++)
99 auto vtkId(MEDCOUPLING2VTKTYPETRADUCER[i]);
100 if(vtkId!=MEDCOUPLING2VTKTYPETRADUCER_NONE)
106 std::map<int,int> ComputeRevMapOfType()
108 std::map<int,int> ret;
109 int nbOfTypesInMC(sizeof(MEDCOUPLING2VTKTYPETRADUCER)/sizeof( decltype(MEDCOUPLING2VTKTYPETRADUCER[0]) ));
110 for(int i=0;i<nbOfTypesInMC;i++)
112 auto vtkId(MEDCOUPLING2VTKTYPETRADUCER[i]);
113 if(vtkId!=MEDCOUPLING2VTKTYPETRADUCER_NONE)
121 vtkInformationDoubleVectorKey *GetMEDReaderMetaDataIfAny()
123 static const char ZE_KEY[]="vtkMEDReader::GAUSS_DATA";
124 MEDCoupling::GlobalDict *gd(MEDCoupling::GlobalDict::GetInstance());
125 if(!gd->hasKey(ZE_KEY))
127 std::string ptSt(gd->value(ZE_KEY));
129 std::istringstream iss(ptSt); iss >> pt;
130 return reinterpret_cast<vtkInformationDoubleVectorKey *>(pt);
133 bool IsInformationOK(vtkInformation *info, std::vector<double>& data)
135 vtkInformationDoubleVectorKey *key(GetMEDReaderMetaDataIfAny());
138 // Check the information contain meta data key
141 int lgth(key->Length(info));
142 const double *data2(info->Get(key));
143 data.insert(data.end(),data2,data2+lgth);
147 void ExtractInfo(vtkInformationVector *inputVector, vtkUnstructuredGrid *& usgIn)
149 vtkInformation *inputInfo(inputVector->GetInformationObject(0));
150 vtkDataSet *input(0);
151 vtkDataSet *input0(vtkDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
152 vtkMultiBlockDataSet *input1(vtkMultiBlockDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
158 throw INTERP_KERNEL::Exception("Input dataSet must be a DataSet or single elt multi block dataset expected !");
159 if(input1->GetNumberOfBlocks()!=1)
160 throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with not exactly one block ! Use MergeBlocks or ExtractBlocks filter before calling this filter !");
161 vtkDataObject *input2(input1->GetBlock(0));
163 throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with exactly one block but this single element is NULL !");
164 vtkDataSet *input2c(vtkDataSet::SafeDownCast(input2));
166 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 !");
170 throw INTERP_KERNEL::Exception("Input data set is NULL !");
171 usgIn=vtkUnstructuredGrid::SafeDownCast(input);
173 throw INTERP_KERNEL::Exception("Input data set is not an unstructured mesh ! This filter works only on unstructured meshes !");
176 DataArrayIdType *ConvertVTKArrayToMCArrayInt(vtkDataArray *data)
179 throw INTERP_KERNEL::Exception("ConvertVTKArrayToMCArrayInt : internal error !");
180 int nbTuples(data->GetNumberOfTuples()),nbComp(data->GetNumberOfComponents());
181 std::size_t nbElts(nbTuples*nbComp);
182 MCAuto<DataArrayIdType> ret(DataArrayIdType::New());
183 ret->alloc(nbTuples,nbComp);
184 for(int i=0;i<nbComp;i++)
186 const char *comp(data->GetComponentName(i));
188 ret->setInfoOnComponent(i,comp);
190 mcIdType *ptOut(ret->getPointer());
191 vtkIntArray *d0(vtkIntArray::SafeDownCast(data));
194 const int *pt(d0->GetPointer(0));
195 std::copy(pt,pt+nbElts,ptOut);
198 vtkLongArray *d1(vtkLongArray::SafeDownCast(data));
201 const long *pt(d1->GetPointer(0));
202 std::copy(pt,pt+nbElts,ptOut);
205 vtkIdTypeArray *d2(vtkIdTypeArray::SafeDownCast(data));
208 const vtkIdType *pt(d2->GetPointer(0));
209 std::copy(pt,pt+nbElts,ptOut);
212 std::ostringstream oss;
213 oss << "ConvertVTKArrayToMCArrayInt : unrecognized array \"" << typeid(*data).name() << "\" type !";
214 throw INTERP_KERNEL::Exception(oss.str());
217 DataArrayDouble *ConvertVTKArrayToMCArrayDouble(vtkDataArray *data)
220 throw INTERP_KERNEL::Exception("ConvertVTKArrayToMCArrayDouble : internal error !");
221 int nbTuples(data->GetNumberOfTuples()),nbComp(data->GetNumberOfComponents());
222 std::size_t nbElts(nbTuples*nbComp);
223 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
224 ret->alloc(nbTuples,nbComp);
225 for(int i=0;i<nbComp;i++)
227 const char *comp(data->GetComponentName(i));
229 ret->setInfoOnComponent(i,comp);
231 double *ptOut(ret->getPointer());
232 vtkFloatArray *d0(vtkFloatArray::SafeDownCast(data));
235 const float *pt(d0->GetPointer(0));
236 for(std::size_t i=0;i<nbElts;i++)
240 vtkDoubleArray *d1(vtkDoubleArray::SafeDownCast(data));
243 const double *pt(d1->GetPointer(0));
244 std::copy(pt,pt+nbElts,ptOut);
247 std::ostringstream oss;
248 oss << "ConvertVTKArrayToMCArrayDouble : unrecognized array \"" << typeid(*data).name() << "\" type !";
249 throw INTERP_KERNEL::Exception(oss.str());
252 DataArray *ConvertVTKArrayToMCArray(vtkDataArray *data)
255 throw INTERP_KERNEL::Exception("ConvertVTKArrayToMCArray : internal error !");
256 vtkFloatArray *d0(vtkFloatArray::SafeDownCast(data));
257 vtkDoubleArray *d1(vtkDoubleArray::SafeDownCast(data));
259 return ConvertVTKArrayToMCArrayDouble(data);
260 vtkIntArray *d2(vtkIntArray::SafeDownCast(data));
261 vtkLongArray *d3(vtkLongArray::SafeDownCast(data));
263 return ConvertVTKArrayToMCArrayInt(data);
264 std::ostringstream oss;
265 oss << "ConvertVTKArrayToMCArray : unrecognized array \"" << typeid(*data).name() << "\" type !";
266 throw INTERP_KERNEL::Exception(oss.str());
269 DataArrayDouble *BuildCoordsFrom(vtkPointSet *ds)
272 throw INTERP_KERNEL::Exception("BuildCoordsFrom : internal error !");
273 vtkPoints *pts(ds->GetPoints());
275 throw INTERP_KERNEL::Exception("BuildCoordsFrom : internal error 2 !");
276 vtkDataArray *data(pts->GetData());
278 throw INTERP_KERNEL::Exception("BuildCoordsFrom : internal error 3 !");
279 MCAuto<DataArrayDouble> coords(ConvertVTKArrayToMCArrayDouble(data));
280 return coords.retn();
283 void ConvertFromUnstructuredGrid(vtkUnstructuredGrid *ds, std::vector< MCAuto<MEDCouplingUMesh> >& ms, std::vector< MCAuto<DataArrayIdType> >& ids)
285 MCAuto<DataArrayDouble> coords(BuildCoordsFrom(ds));
286 vtkIdType nbCells(ds->GetNumberOfCells());
287 vtkCellArray *ca(ds->GetCells());
290 //vtkIdType nbEnt(ca->GetNumberOfConnectivityEntries()); // todo: unused
291 //vtkIdType *caPtr(ca->GetData()->GetPointer(0)); // todo: unused
292 vtkUnsignedCharArray *ct(ds->GetCellTypesArray());
294 throw INTERP_KERNEL::Exception("ConvertFromUnstructuredGrid : internal error");
295 vtkIdTypeArray *cla(ds->GetCellLocationsArray());
296 //const vtkIdType *claPtr(cla->GetPointer(0)); // todo: unused
298 throw INTERP_KERNEL::Exception("ConvertFromUnstructuredGrid : internal error 2");
299 const unsigned char *ctPtr(ct->GetPointer(0));
300 std::map<int,int> m(ComputeMapOfType());
301 MCAuto<DataArrayInt> lev(DataArrayInt::New()) ; lev->alloc(nbCells,1);
302 int *levPtr(lev->getPointer());
303 for(vtkIdType i=0;i<nbCells;i++)
305 std::map<int,int>::iterator it(m.find(ctPtr[i]));
308 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)(*it).second));
309 levPtr[i]=cm.getDimension();
313 std::ostringstream oss; oss << "ConvertFromUnstructuredGrid : at pos #" << i << " unrecognized VTK cell with type =" << ctPtr[i];
314 throw INTERP_KERNEL::Exception(oss.str());
317 MCAuto<DataArrayInt> levs(lev->getDifferentValues());
318 //vtkIdTypeArray *faces(ds->GetFaces()),*faceLoc(ds->GetFaceLocations()); // todo: unused
319 for(const int *curLev=levs->begin();curLev!=levs->end();curLev++)
321 MCAuto<MEDCouplingUMesh> m0(MEDCouplingUMesh::New("",*curLev));
322 m0->setCoords(coords); m0->allocateCells();
323 MCAuto<DataArrayIdType> cellIdsCurLev(lev->findIdsEqual(*curLev));
324 for(const mcIdType *cellId=cellIdsCurLev->begin();cellId!=cellIdsCurLev->end();cellId++)
326 std::map<int,int>::iterator it(m.find(ds->GetCellType(*cellId)));
327 INTERP_KERNEL::NormalizedCellType ct((INTERP_KERNEL::NormalizedCellType)(*it).second);
328 if(ct!=INTERP_KERNEL::NORM_POLYHED)
331 const vtkIdType *ptszz(nullptr);
332 ds->GetCellPoints(*cellId, szzz, ptszz);
333 std::vector<mcIdType> conn(ptszz,ptszz+szzz);
334 m0->insertNextCell(ct,szzz,conn.data());
338 // # de faces du polyèdre
339 vtkIdType nbOfFaces(0);
340 // connectivé des faces (numFace0Pts, id1, id2, id3, numFace1Pts,id1, id2, id3, ...)
341 const vtkIdType *facPtr(nullptr);
342 ds->GetFaceStream(*cellId, nbOfFaces, facPtr);
343 std::vector<mcIdType> conn;
344 for(vtkIdType k=0;k<nbOfFaces;k++)
346 vtkIdType nbOfNodesInFace(*facPtr++);
347 std::copy(facPtr,facPtr+nbOfNodesInFace,std::back_inserter(conn));
350 facPtr+=nbOfNodesInFace;
352 m0->insertNextCell(ct,ToIdType(conn.size()),&conn[0]);
355 ms.push_back(m0); ids.push_back(cellIdsCurLev);
359 vtkSmartPointer<vtkUnstructuredGrid> ConvertUMeshFromMCToVTK(const MEDCouplingUMesh *mVor)
361 std::map<int,int> zeMapRev(ComputeRevMapOfType());
362 int nbCells(mVor->getNumberOfCells());
363 vtkSmartPointer<vtkUnstructuredGrid> ret(vtkSmartPointer<vtkUnstructuredGrid>::New());
366 vtkSmartPointer<vtkPoints> points(vtkSmartPointer<vtkPoints>::New());
368 const DataArrayDouble *vorCoords(mVor->getCoords());
369 vtkSmartPointer<vtkDoubleArray> da(vtkSmartPointer<vtkDoubleArray>::New());
370 da->SetNumberOfComponents((vtkIdType)vorCoords->getNumberOfComponents());
371 da->SetNumberOfTuples((vtkIdType)vorCoords->getNumberOfTuples());
372 std::copy(vorCoords->begin(),vorCoords->end(),da->GetPointer(0));
375 mVor->checkConsistencyLight();
376 switch(mVor->getMeshDimension())
380 vtkIdType *cPtr(nullptr),*dPtr(nullptr);
381 unsigned char *aPtr(nullptr);
382 vtkSmartPointer<vtkUnsignedCharArray> cellTypes(vtkSmartPointer<vtkUnsignedCharArray>::New());
384 cellTypes->SetNumberOfComponents(1);
385 cellTypes->SetNumberOfTuples(nbCells);
386 aPtr=cellTypes->GetPointer(0);
388 vtkSmartPointer<vtkIdTypeArray> cellLocations(vtkSmartPointer<vtkIdTypeArray>::New());
390 cellLocations->SetNumberOfComponents(1);
391 cellLocations->SetNumberOfTuples(nbCells);
392 cPtr=cellLocations->GetPointer(0);
394 vtkSmartPointer<vtkIdTypeArray> cells(vtkSmartPointer<vtkIdTypeArray>::New());
396 MCAuto<DataArrayIdType> tmp2(mVor->computeEffectiveNbOfNodesPerCell());
397 cells->SetNumberOfComponents(1);
398 cells->SetNumberOfTuples(tmp2->accumulate((std::size_t)0)+nbCells);
399 dPtr=cells->GetPointer(0);
401 const mcIdType *connPtr(mVor->getNodalConnectivity()->begin()),*connIPtr(mVor->getNodalConnectivityIndex()->begin());
403 std::vector<vtkIdType> ee,ff;
404 for(int i=0;i<nbCells;i++,connIPtr++)
406 INTERP_KERNEL::NormalizedCellType ct(static_cast<INTERP_KERNEL::NormalizedCellType>(connPtr[connIPtr[0]]));
407 *aPtr++=(unsigned char)zeMapRev[connPtr[connIPtr[0]]];
408 if(ct!=INTERP_KERNEL::NORM_POLYHED)
410 int sz(connIPtr[1]-connIPtr[0]-1);
412 dPtr=std::copy(connPtr+connIPtr[0]+1,connPtr+connIPtr[1],dPtr);
418 std::set<int> s(connPtr+connIPtr[0]+1,connPtr+connIPtr[1]); s.erase(-1);
419 vtkIdType nbFace((vtkIdType)(std::count(connPtr+connIPtr[0]+1,connPtr+connIPtr[1],-1)+1));
420 ff.push_back(nbFace);
421 const mcIdType *work(connPtr+connIPtr[0]+1);
422 for(int j=0;j<nbFace;j++)
424 const mcIdType *work2=std::find(work,connPtr+connIPtr[1],-1);
425 ff.push_back((vtkIdType)std::distance(work,work2));
426 ff.insert(ff.end(),work,work2);
429 *dPtr++=(int)s.size();
430 dPtr=std::copy(s.begin(),s.end(),dPtr);
431 *cPtr++=k; k+=(int)s.size()+1;
432 ee.push_back(kk); kk+=connIPtr[1]-connIPtr[0]+1;
436 vtkSmartPointer<vtkIdTypeArray> faceLocations(vtkSmartPointer<vtkIdTypeArray>::New());
438 faceLocations->SetNumberOfComponents(1);
439 faceLocations->SetNumberOfTuples((vtkIdType)ee.size());
440 std::copy(ee.begin(),ee.end(),faceLocations->GetPointer(0));
442 vtkSmartPointer<vtkIdTypeArray> faces(vtkSmartPointer<vtkIdTypeArray>::New());
444 faces->SetNumberOfComponents(1);
445 faces->SetNumberOfTuples((vtkIdType)ff.size());
446 std::copy(ff.begin(),ff.end(),faces->GetPointer(0));
448 vtkSmartPointer<vtkCellArray> cells2(vtkSmartPointer<vtkCellArray>::New());
449 cells2->SetCells(nbCells,cells);
450 ret->SetCells(cellTypes,cellLocations,cells2,faceLocations,faces);
455 vtkSmartPointer<vtkUnsignedCharArray> cellTypes(vtkSmartPointer<vtkUnsignedCharArray>::New());
457 cellTypes->SetNumberOfComponents(1);
458 cellTypes->SetNumberOfTuples(nbCells);
459 unsigned char *ptr(cellTypes->GetPointer(0));
460 std::fill(ptr,ptr+nbCells,zeMapRev[(int)INTERP_KERNEL::NORM_POLYGON]);
462 vtkIdType *cPtr(nullptr),*dPtr(nullptr);
463 vtkSmartPointer<vtkIdTypeArray> cellLocations(vtkSmartPointer<vtkIdTypeArray>::New());
465 cellLocations->SetNumberOfComponents(1);
466 cellLocations->SetNumberOfTuples(nbCells);
467 cPtr=cellLocations->GetPointer(0);
469 vtkSmartPointer<vtkIdTypeArray> cells(vtkSmartPointer<vtkIdTypeArray>::New());
471 cells->SetNumberOfComponents(1);
472 cells->SetNumberOfTuples(mVor->getNodalConnectivity()->getNumberOfTuples());
473 dPtr=cells->GetPointer(0);
475 const mcIdType *connPtr(mVor->getNodalConnectivity()->begin()),*connIPtr(mVor->getNodalConnectivityIndex()->begin());
477 for(int i=0;i<nbCells;i++,connIPtr++)
479 *dPtr++=connIPtr[1]-connIPtr[0]-1;
480 dPtr=std::copy(connPtr+connIPtr[0]+1,connPtr+connIPtr[1],dPtr);
481 *cPtr++=k; k+=connIPtr[1]-connIPtr[0];
483 vtkSmartPointer<vtkCellArray> cells2(vtkSmartPointer<vtkCellArray>::New());
484 cells2->SetCells(nbCells,cells);
485 ret->SetCells(cellTypes,cellLocations,cells2);
490 vtkSmartPointer<vtkUnsignedCharArray> cellTypes(vtkSmartPointer<vtkUnsignedCharArray>::New());
492 cellTypes->SetNumberOfComponents(1);
493 cellTypes->SetNumberOfTuples(nbCells);
494 unsigned char *ptr(cellTypes->GetPointer(0));
495 std::fill(ptr,ptr+nbCells,zeMapRev[(int)INTERP_KERNEL::NORM_SEG2]);
497 vtkIdType *cPtr(nullptr),*dPtr(nullptr);
498 vtkSmartPointer<vtkIdTypeArray> cellLocations(vtkSmartPointer<vtkIdTypeArray>::New());
500 cellLocations->SetNumberOfComponents(1);
501 cellLocations->SetNumberOfTuples(nbCells);
502 cPtr=cellLocations->GetPointer(0);
504 vtkSmartPointer<vtkIdTypeArray> cells(vtkSmartPointer<vtkIdTypeArray>::New());
506 cells->SetNumberOfComponents(1);
507 cells->SetNumberOfTuples(mVor->getNodalConnectivity()->getNumberOfTuples());
508 dPtr=cells->GetPointer(0);
510 const mcIdType *connPtr(mVor->getNodalConnectivity()->begin()),*connIPtr(mVor->getNodalConnectivityIndex()->begin());
511 for(int i=0;i<nbCells;i++,connIPtr++)
514 dPtr=std::copy(connPtr+connIPtr[0]+1,connPtr+connIPtr[1],dPtr);
517 vtkSmartPointer<vtkCellArray> cells2(vtkSmartPointer<vtkCellArray>::New());
518 cells2->SetCells(nbCells,cells);
519 ret->SetCells(cellTypes,cellLocations,cells2);
523 throw INTERP_KERNEL::Exception("Not implemented yet !");
525 ret->SetPoints(points);
532 OffsetKeeper():_vtk_arr(0) { }
533 void pushBack(vtkDataArray *da) { _da_on.push_back(da); }
534 void setVTKArray(vtkIdTypeArray *arr) {
535 MCAuto<DataArrayIdType> offmc(ConvertVTKArrayToMCArrayInt(arr));
536 _off_arr=offmc; _vtk_arr=arr; }
537 const std::vector<vtkDataArray *>& getArrayGauss() const { return _da_on; }
538 const DataArrayIdType *getOffsets() const { return _off_arr; }
539 vtkIdTypeArray *getVTKOffsets() const { return _vtk_arr; }
541 std::vector<vtkDataArray *> _da_on;
542 MCAuto<DataArrayIdType> _off_arr;
543 vtkIdTypeArray *_vtk_arr;
546 void FillAdvInfoFrom(int vtkCT, const std::vector<double>& GaussAdvData, int nbGaussPt, int nbNodesPerCell, std::vector<double>& refCoo,std::vector<double>& posInRefCoo)
548 int nbOfCTS((int)GaussAdvData[0]),pos(1);
549 for(int i=0;i<nbOfCTS;i++)
551 int lgth((int)GaussAdvData[pos]);
552 int curCT((int)GaussAdvData[pos+1]),dim((int)GaussAdvData[pos+2]);
558 int lgthExp(nbNodesPerCell*dim+nbGaussPt*dim);
559 if(lgth!=lgthExp+2)//+2 for cell type and dimension !
561 std::ostringstream oss; oss << "FillAdvInfoFrom : Internal error. Unmatch with MEDReader version ? Expect size " << lgthExp << " and have " << lgth << " !";
562 throw INTERP_KERNEL::Exception(oss.str());
564 refCoo.insert(refCoo.end(),GaussAdvData.begin()+pos+3,GaussAdvData.begin()+pos+3+nbNodesPerCell*dim);
565 posInRefCoo.insert(posInRefCoo.end(),GaussAdvData.begin()+pos+3+nbNodesPerCell*dim,GaussAdvData.begin()+pos+3+nbNodesPerCell*dim+nbGaussPt*dim);
566 //std::copy(refCoo.begin(),refCoo.end(),std::ostream_iterator<double>(std::cerr," ")); std::cerr << std::endl;
567 //std::copy(posInRefCoo.begin(),posInRefCoo.end(),std::ostream_iterator<double>(std::cerr," ")); std::cerr << std::endl;
570 std::ostringstream oss; oss << "FillAdvInfoFrom : Internal error ! Not found cell type " << vtkCT << " in advanced Gauss info !";
571 throw INTERP_KERNEL::Exception(oss.str());
574 template<class T, class U>
575 vtkSmartPointer<T> ExtractFieldFieldArr(T *elt2, int sizeOfOutArr, int nbOfCellsOfInput, const mcIdType *offsetsPtr, const mcIdType *nbPtsPerCellPtr)
577 vtkSmartPointer<T> elt3(vtkSmartPointer<T>::New());
578 int nbc(elt2->GetNumberOfComponents());
579 elt3->SetNumberOfComponents(nbc);
580 elt3->SetNumberOfTuples(sizeOfOutArr);
581 for(int i=0;i<nbc;i++)
583 const char *name(elt2->GetComponentName(i));
585 elt3->SetComponentName(i,name);
587 elt3->SetName(elt2->GetName());
589 U *ptr(elt3->GetPointer(0));
590 const U *srcPtr(elt2->GetPointer(0));
591 for(int i=0;i<nbOfCellsOfInput;i++)
592 ptr=std::copy(srcPtr+nbc*offsetsPtr[i],srcPtr+nbc*(offsetsPtr[i]+nbPtsPerCellPtr[i]),ptr);
596 template<class T, class U>
597 vtkSmartPointer<T> ExtractCellFieldArr(T *elt2, int sizeOfOutArr, int nbOfCellsOfInput, const mcIdType *idsPtr, const mcIdType *nbPtsPerCellPtr)
599 vtkSmartPointer<T> elt3(vtkSmartPointer<T>::New());
600 int nbc(elt2->GetNumberOfComponents());
601 elt3->SetNumberOfComponents(nbc);
602 elt3->SetNumberOfTuples(sizeOfOutArr);
603 for(int i=0;i<nbc;i++)
605 const char *name(elt2->GetComponentName(i));
607 elt3->SetComponentName(i,name);
609 elt3->SetName(elt2->GetName());
611 U *ptr(elt3->GetPointer(0));
612 const U *srcPtr(elt2->GetPointer(0));
613 for(int i=0;i<nbOfCellsOfInput;i++)
614 for(int j=0;j<nbPtsPerCellPtr[i];j++)
615 ptr=std::copy(srcPtr+nbc*idsPtr[i],srcPtr+nbc*(idsPtr[i]+1),ptr);
619 vtkSmartPointer<vtkUnstructuredGrid> Voronize(const MEDCouplingUMesh *m, const DataArrayIdType *ids, vtkIdTypeArray *vtkOff, const DataArrayIdType *offsetsBase, const std::vector<vtkDataArray *>& arrGauss, const std::vector<double>& GaussAdvData, const std::vector<vtkDataArray *>& arrsOnCells)
622 throw INTERP_KERNEL::Exception("Voronize : no Gauss array !");
623 int nbTuples(arrGauss[0]->GetNumberOfTuples());
624 for(std::vector<vtkDataArray *>::const_iterator it=arrGauss.begin();it!=arrGauss.end();it++)
626 if((*it)->GetNumberOfTuples()!=nbTuples)
628 std::ostringstream oss; oss << "Mismatch of number of tuples in Gauss arrays for array \"" << (*it)->GetName() << "\"";
629 throw INTERP_KERNEL::Exception(oss.str());
632 // Look at vtkOff has in the stomac
633 vtkInformation *info(vtkOff->GetInformation());
635 throw INTERP_KERNEL::Exception("info is null ! Internal error ! Looks bad !");
636 vtkInformationQuadratureSchemeDefinitionVectorKey *key(vtkQuadratureSchemeDefinition::DICTIONARY());
638 throw INTERP_KERNEL::Exception("No quadrature key in info included in offets array ! Internal error ! Looks bad !");
639 int dictSize(key->Size(info));
640 INTERP_KERNEL::AutoPtr<vtkQuadratureSchemeDefinition *> dict(new vtkQuadratureSchemeDefinition *[dictSize]);
641 key->GetRange(info,dict,0,0,dictSize);
643 MCAuto<MEDCouplingFieldDouble> field(MEDCouplingFieldDouble::New(ON_GAUSS_PT));
646 int nbOfCellsOfInput(m->getNumberOfCells());
647 MCAuto<DataArrayIdType> nbPtsPerCellArr(DataArrayIdType::New()); nbPtsPerCellArr->alloc(nbOfCellsOfInput,1);
648 std::map<int,int> zeMapRev(ComputeRevMapOfType()),zeMap(ComputeMapOfType());
649 std::set<INTERP_KERNEL::NormalizedCellType> agt(m->getAllGeoTypes());
650 for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=agt.begin();it!=agt.end();it++)
652 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(*it));
653 std::map<int,int>::const_iterator it2(zeMapRev.find((int)*it));
654 if(it2==zeMapRev.end())
655 throw INTERP_KERNEL::Exception("Internal error ! no type conversion available !");
656 vtkQuadratureSchemeDefinition *gaussLoc(dict[(*it2).second]);
659 std::ostringstream oss; oss << "For cell type " << cm.getRepr() << " no Gauss info !";
660 throw INTERP_KERNEL::Exception(oss.str());
662 int np(gaussLoc->GetNumberOfQuadraturePoints()),nbPtsPerCell((int)cm.getNumberOfNodes());
663 const double /**sfw(gaussLoc->GetShapeFunctionWeights()),*/*w(gaussLoc->GetQuadratureWeights());; // todo: sfw is unused
664 std::vector<double> refCoo,posInRefCoo,wCpp(w,w+np);
665 FillAdvInfoFrom((*it2).second,GaussAdvData,np,nbPtsPerCell,refCoo,posInRefCoo);
666 field->setGaussLocalizationOnType(*it,refCoo,posInRefCoo,wCpp);
667 MCAuto<DataArrayIdType> ids2(m->giveCellsWithType(*it));
668 nbPtsPerCellArr->setPartOfValuesSimple3(np,ids2->begin(),ids2->end(),0,1,1);
670 int zeSizeOfOutCellArr(nbPtsPerCellArr->accumulate((std::size_t)0));
671 { MCAuto<DataArrayDouble> fakeArray(DataArrayDouble::New()); fakeArray->alloc(zeSizeOfOutCellArr,1); field->setArray(fakeArray); }
672 field->checkConsistencyLight();
673 MCAuto<MEDCouplingFieldDouble> vor(field->voronoize(1e-12));// The key is here !
674 MEDCouplingUMesh *mVor(dynamic_cast<MEDCouplingUMesh *>(vor->getMesh()));
676 vtkSmartPointer<vtkUnstructuredGrid> ret(ConvertUMeshFromMCToVTK(mVor));
678 MCAuto<DataArrayIdType> myOffsets(offsetsBase->selectByTupleIdSafe(ids->begin(),ids->end()));
679 const mcIdType *myOffsetsPtr(myOffsets->begin()),*nbPtsPerCellArrPtr(nbPtsPerCellArr->begin());
680 for(std::vector<vtkDataArray *>::const_iterator it=arrGauss.begin();it!=arrGauss.end();it++)
682 vtkDataArray *elt(*it);
683 vtkDoubleArray *elt2(vtkDoubleArray::SafeDownCast(elt));
684 vtkIntArray *elt3(vtkIntArray::SafeDownCast(elt));
685 vtkIdTypeArray *elt4(vtkIdTypeArray::SafeDownCast(elt));
688 vtkSmartPointer<vtkDoubleArray> arr(ExtractFieldFieldArr<vtkDoubleArray,double>(elt2,zeSizeOfOutCellArr,nbOfCellsOfInput,myOffsetsPtr,nbPtsPerCellArrPtr));
689 ret->GetCellData()->AddArray(arr);
694 vtkSmartPointer<vtkIntArray> arr(ExtractFieldFieldArr<vtkIntArray,int>(elt3,zeSizeOfOutCellArr,nbOfCellsOfInput,myOffsetsPtr,nbPtsPerCellArrPtr));
695 ret->GetCellData()->AddArray(arr);
700 vtkSmartPointer<vtkIdTypeArray> arr(ExtractFieldFieldArr<vtkIdTypeArray,vtkIdType>(elt4,zeSizeOfOutCellArr,nbOfCellsOfInput,myOffsetsPtr,nbPtsPerCellArrPtr));
701 ret->GetCellData()->AddArray(arr);
705 for(std::vector<vtkDataArray *>::const_iterator it=arrsOnCells.begin();it!=arrsOnCells.end();it++)
707 vtkDataArray *elt(*it);
708 vtkDoubleArray *elt2(vtkDoubleArray::SafeDownCast(elt));
709 vtkIntArray *elt3(vtkIntArray::SafeDownCast(elt));
710 vtkIdTypeArray *elt4(vtkIdTypeArray::SafeDownCast(elt));
713 vtkSmartPointer<vtkDoubleArray> arr(ExtractCellFieldArr<vtkDoubleArray,double>(elt2,zeSizeOfOutCellArr,nbOfCellsOfInput,ids->begin(),nbPtsPerCellArrPtr));
714 ret->GetCellData()->AddArray(arr);
719 vtkSmartPointer<vtkIntArray> arr(ExtractCellFieldArr<vtkIntArray,int>(elt3,zeSizeOfOutCellArr,nbOfCellsOfInput,ids->begin(),nbPtsPerCellArrPtr));
720 ret->GetCellData()->AddArray(arr);
725 vtkSmartPointer<vtkIdTypeArray> arr(ExtractCellFieldArr<vtkIdTypeArray,vtkIdType>(elt4,zeSizeOfOutCellArr,nbOfCellsOfInput,ids->begin(),nbPtsPerCellArrPtr));
726 ret->GetCellData()->AddArray(arr);
733 vtkSmartPointer<vtkUnstructuredGrid> ComputeVoroGauss(vtkUnstructuredGrid *usgIn, const std::vector<double>& GaussAdvData)
735 OffsetKeeper zeOffsets;
736 std::string zeArrOffset;
737 std::vector<std::string> cellFieldNamesToDestroy;
739 int nArrays(usgIn->GetFieldData()->GetNumberOfArrays());
740 for(int i=0;i<nArrays;i++)
742 vtkDataArray *array(usgIn->GetFieldData()->GetArray(i));
745 const char* arrayOffsetName(array->GetInformation()->Get(vtkQuadratureSchemeDefinition::QUADRATURE_OFFSET_ARRAY_NAME()));
748 std::string arrOffsetNameCpp(arrayOffsetName);
749 cellFieldNamesToDestroy.push_back(arrOffsetNameCpp);
750 if(arrOffsetNameCpp.find("ELGA@")==std::string::npos)
752 if(zeArrOffset.empty())
753 zeArrOffset=arrOffsetNameCpp;
755 if(zeArrOffset!=arrOffsetNameCpp)
757 throw INTERP_KERNEL::Exception("ComputeVoroGauss : error in QUADRATURE_OFFSET_ARRAY_NAME for Gauss fields array !");
759 zeOffsets.pushBack(array);
761 if(zeArrOffset.empty())
762 throw INTERP_KERNEL::Exception("ComputeVoroGauss : no Gauss points fields in DataSet !");
764 std::vector< MCAuto<MEDCouplingUMesh> > ms;
765 std::vector< MCAuto<DataArrayIdType> > ids;
766 ConvertFromUnstructuredGrid(usgIn,ms,ids);
768 vtkDataArray *offTmp(usgIn->GetCellData()->GetArray(zeArrOffset.c_str()));
771 std::ostringstream oss; oss << "ComputeVoroGauss : cell field " << zeArrOffset << " not found !";
772 throw INTERP_KERNEL::Exception(oss.str());
774 vtkIdTypeArray *offsets(vtkIdTypeArray::SafeDownCast(offTmp));
777 std::ostringstream oss; oss << "ComputeVoroGauss : cell field " << zeArrOffset << " exists but not with the right type of data !";
778 throw INTERP_KERNEL::Exception(oss.str());
781 zeOffsets.setVTKArray(offsets);
784 std::vector<vtkDataArray *> arrsOnCells;
786 int nArrays(usgIn->GetCellData()->GetNumberOfArrays());
787 for(int i=0;i<nArrays;i++)
789 vtkDataArray *array(usgIn->GetCellData()->GetArray(i));
792 std::string name(array->GetName());
793 if(std::find(cellFieldNamesToDestroy.begin(),cellFieldNamesToDestroy.end(),name)==cellFieldNamesToDestroy.end())
795 arrsOnCells.push_back(array);
799 vtkDataArray *doNotKeepThis(zeOffsets.getVTKOffsets());
800 std::vector<vtkDataArray *>::iterator it2(std::find(arrsOnCells.begin(),arrsOnCells.end(),doNotKeepThis));
801 if(it2!=arrsOnCells.end())
802 arrsOnCells.erase(it2);
806 std::size_t sz(ms.size());
807 std::vector< vtkSmartPointer<vtkUnstructuredGrid> > res;
808 for(std::size_t i=0;i<sz;i++)
810 MCAuto<MEDCouplingUMesh> mmc(ms[i]);
811 MCAuto<DataArrayIdType> myIds(ids[i]);
812 vtkSmartPointer<vtkUnstructuredGrid> vor(Voronize(mmc,myIds,zeOffsets.getVTKOffsets(),zeOffsets.getOffsets(),zeOffsets.getArrayGauss(),GaussAdvData,arrsOnCells));
816 throw INTERP_KERNEL::Exception("Dataset is empty !");
817 vtkSmartPointer<vtkMultiBlockDataGroupFilter> mb(vtkSmartPointer<vtkMultiBlockDataGroupFilter>::New());
818 vtkSmartPointer<vtkCompositeDataToUnstructuredGridFilter> cd(vtkSmartPointer<vtkCompositeDataToUnstructuredGridFilter>::New());
819 for(std::vector< vtkSmartPointer<vtkUnstructuredGrid> >::const_iterator it=res.begin();it!=res.end();it++)
820 mb->AddInputData(*it);
821 cd->SetInputConnection(mb->GetOutputPort());
822 cd->SetMergePoints(0);
824 vtkSmartPointer<vtkUnstructuredGrid> ret;
831 vtkVoroGauss::vtkVoroGauss()
833 this->SetNumberOfInputPorts(1);
834 this->SetNumberOfOutputPorts(1);
837 vtkVoroGauss::~vtkVoroGauss()
841 int vtkVoroGauss::RequestInformation(vtkInformation * /*request*/, vtkInformationVector **inputVector, vtkInformationVector * /*outputVector*/)
843 //std::cerr << "########################################## vtkVoroGauss::RequestInformation ##########################################" << std::endl;
846 vtkUnstructuredGrid *usgIn(0);
847 ExtractInfo(inputVector[0],usgIn);
849 catch(INTERP_KERNEL::Exception& e)
851 std::ostringstream oss;
852 oss << "Exception has been thrown in vtkVoroGauss::RequestInformation : " << e.what() << std::endl;
853 if(this->HasObserver("ErrorEvent") )
854 this->InvokeEvent("ErrorEvent",const_cast<char *>(oss.str().c_str()));
856 vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
857 vtkObject::BreakOnError();
863 int vtkVoroGauss::RequestData(vtkInformation * /*request*/, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
865 //std::cerr << "########################################## vtkVoroGauss::RequestData ##########################################" << std::endl;
869 std::vector<double> GaussAdvData;
870 bool isOK(IsInformationOK(inputVector[0]->GetInformationObject(0),GaussAdvData));
872 throw INTERP_KERNEL::Exception("Sorry but no advanced gauss info found ! Expect to be called right after a MEDReader containing Gauss Points !");
873 vtkUnstructuredGrid *usgIn(0);
874 ExtractInfo(inputVector[0],usgIn);
876 vtkSmartPointer<vtkUnstructuredGrid> ret(ComputeVoroGauss(usgIn,GaussAdvData));
877 //vtkInformation *inInfo(inputVector[0]->GetInformationObject(0)); // todo: unused
878 vtkInformation *outInfo(outputVector->GetInformationObject(0));
879 vtkUnstructuredGrid *output(vtkUnstructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())));
880 output->ShallowCopy(ret);
882 catch(INTERP_KERNEL::Exception& e)
884 std::ostringstream oss;
885 oss << "Exception has been thrown in vtkVoroGauss::RequestData : " << e.what() << std::endl;
886 if(this->HasObserver("ErrorEvent") )
887 this->InvokeEvent("ErrorEvent",const_cast<char *>(oss.str().c_str()));
889 vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
890 vtkObject::BreakOnError();
896 void vtkVoroGauss::PrintSelf(ostream& os, vtkIndent indent)
898 this->Superclass::PrintSelf(os, indent);