src/TorseurCIH/plugin/TorseurCIHModule/vtkTorseurCIH.cxx

   1 // Copyright (C) 2021  CEA/DEN, EDF R&D
   2 //
   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.
   7 //
   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.
  12 //
  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
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (EDF R&D)
  20
  21 #include "vtkTorseurCIH.h"
  22
  23 #include <vtkCell.h>
  24 #include <vtkCellArray.h>
  25 #include <vtkCellData.h>
  26 #include <vtkCellType.h>
  27 #include <vtkIdTypeArray.h>
  28 #include <vtkIntArray.h>
  29 #include <vtkPointData.h>
  30 #include <vtkMultiBlockDataSet.h>
  31 #include <vtkStreamingDemandDrivenPipeline.h>
  32 #include <vtkUnstructuredGrid.h>
  33 #include <vtkAlgorithmOutput.h>
  34 #include <vtkCompositeDataToUnstructuredGridFilter.h>
  35 #include <vtkDataArraySelection.h>
  36 #include <vtkDataObjectTreeIterator.h>
  37 #include <vtkDataSet.h>
  38 #include <vtkDataSetAttributes.h>
  39 #include <vtkDemandDrivenPipeline.h>
  40 #include <vtkDoubleArray.h>
  41 #include <vtkExecutive.h>
  42 #include <vtkFloatArray.h>
  43 #include <vtkInformation.h>
  44 #include <vtkInformationDataObjectKey.h>
  45 #include <vtkInformationDataObjectMetaDataKey.h>
  46 #include <vtkInformationDoubleVectorKey.h>
  47 #include <vtkInformationQuadratureSchemeDefinitionVectorKey.h>
  48 #include <vtkInformationStringKey.h>
  49 #include <vtkInformationVector.h>
  50 #include <vtkLongArray.h>
  51 #include <vtkMultiBlockDataGroupFilter.h>
  52 #include <vtkMultiBlockDataSet.h>
  53 #include <vtkMutableDirectedGraph.h>
  54 #include <vtkObjectFactory.h>
  55 #include <vtkQuadratureSchemeDefinition.h>
  56 #include <vtkStringArray.h>
  57 #include <vtkTable.h>
  58 #include <vtkUnsignedCharArray.h>
  59
  60 #include "InterpKernelAutoPtr.hxx"
  61 #include "InterpKernelGaussCoords.hxx"
  62 #include "MEDCouplingFieldDouble.hxx"
  63 #include "MEDCouplingMemArray.hxx"
  64 #include "MEDCouplingUMesh.hxx"
  65
  66 #include <deque>
  67 #include <map>
  68 #include <set>
  69 #include <sstream>
  70
  71 using MEDCoupling::DataArray;
  72 using MEDCoupling::DataArrayDouble;
  73 using MEDCoupling::DataArrayInt;
  74 using MEDCoupling::DataArrayInt64;
  75 using MEDCoupling::DynamicCastSafe;
  76 using MEDCoupling::MCAuto;
  77 using MEDCoupling::MEDCouplingFieldDouble;
  78 using MEDCoupling::MEDCouplingMesh;
  79 using MEDCoupling::MEDCouplingUMesh;
  80 using MEDCoupling::ON_GAUSS_PT;
  81
  82 vtkStandardNewMacro(vtkTorseurCIH);
  83 ///////////////////
  84
  85 std::map<int, int> ComputeMapOfType()
  86 {
  87   std::map<int, int> ret;
  88   int nbOfTypesInMC(sizeof(MEDCOUPLING2VTKTYPETRADUCER) / sizeof(int));
  89   for (int i = 0; i < nbOfTypesInMC; i++)
  90   {
  91     int vtkId(MEDCOUPLING2VTKTYPETRADUCER[i]);
  92     if (vtkId != -1)
  93       ret[vtkId] = i;
  94   }
  95   return ret;
  96 }
  97
  98 std::map<int, int> ComputeRevMapOfType()
  99 {
 100   std::map<int, int> ret;
 101   int nbOfTypesInMC(sizeof(MEDCOUPLING2VTKTYPETRADUCER) / sizeof(int));
 102   for (int i = 0; i < nbOfTypesInMC; i++)
 103   {
 104     int vtkId(MEDCOUPLING2VTKTYPETRADUCER[i]);
 105     if (vtkId != -1)
 106       ret[i] = vtkId;
 107   }
 108   return ret;
 109 }
 110
 111 ///////////////////
 112
 113 void ExtractInfo(vtkInformationVector* inputVector, vtkSmartPointer<vtkUnstructuredGrid>& usgIn)
 114 {
 115   vtkInformation* inputInfo(inputVector->GetInformationObject(0));
 116   vtkDataSet* input = nullptr;
 117   vtkDataSet* input0(vtkDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
 118   vtkMultiBlockDataSet* input1(
 119     vtkMultiBlockDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
 120   if (input0)
 121     input = input0;
 122   else
 123   {
 124     if (!input1)
 125       throw INTERP_KERNEL::Exception(
 126         "Input dataSet must be a DataSet or single elt multi block dataset expected !");
 127     if (input1->GetNumberOfBlocks() != 1)
 128       throw INTERP_KERNEL::Exception(
 129         "Input dataSet is a multiblock dataset with not exactly one block ! Use MergeBlocks or "
 130         "ExtractBlocks filter before calling this filter !");
 131     vtkDataObject* input2(input1->GetBlock(0));
 132     if (!input2)
 133       throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with exactly one block "
 134                                      "but this single element is NULL !");
 135     vtkDataSet* input2c(vtkDataSet::SafeDownCast(input2));
 136     if (!input2c)
 137       throw INTERP_KERNEL::Exception(
 138         "Input dataSet is a multiblock dataset with exactly one block but this single element is "
 139         "not a dataset ! Use MergeBlocks or ExtractBlocks filter before calling this filter !");
 140     input = input2c;
 141   }
 142   if (!input)
 143     throw INTERP_KERNEL::Exception("Input data set is NULL !");
 144   usgIn.TakeReference(vtkUnstructuredGrid::SafeDownCast(input));
 145   if (!usgIn.Get())
 146   {
 147     if (!input1)
 148     {
 149       vtkNew<vtkMultiBlockDataGroupFilter> mb;
 150       vtkNew<vtkCompositeDataToUnstructuredGridFilter> cd;
 151       mb->AddInputData(input);
 152       cd->SetInputConnection(mb->GetOutputPort());
 153       cd->SetMergePoints(0);
 154       cd->Update();
 155       usgIn = cd->GetOutput();
 156     }
 157     else
 158     {
 159       vtkNew<vtkCompositeDataToUnstructuredGridFilter> filter;
 160       filter->SetMergePoints(0);
 161       filter->SetInputData(input1);
 162       filter->Update();
 163       vtkUnstructuredGrid* res(filter->GetOutput());
 164       usgIn.TakeReference(res);
 165       if (res)
 166         res->Register(nullptr);
 167     }
 168   }
 169   else
 170     usgIn->Register(nullptr);
 171 }
 172
 173 DataArrayInt* ConvertVTKArrayToMCArrayInt(vtkDataArray* data)
 174 {
 175   if (!data)
 176     throw INTERP_KERNEL::Exception("ConvertVTKArrayToMCArrayInt : internal error !");
 177   int nbTuples(data->GetNumberOfTuples()), nbComp(data->GetNumberOfComponents());
 178   std::size_t nbElts(nbTuples * nbComp);
 179   MCAuto<DataArrayInt> ret(DataArrayInt::New());
 180   ret->alloc(nbTuples, nbComp);
 181   for (int i = 0; i < nbComp; i++)
 182   {
 183     const char* comp(data->GetComponentName(i));
 184     if (comp)
 185       ret->setInfoOnComponent(i, comp);
 186   }
 187   int* ptOut(ret->getPointer());
 188   vtkIntArray* d0(vtkIntArray::SafeDownCast(data));
 189   if (d0)
 190   {
 191     const int* pt(d0->GetPointer(0));
 192     std::copy(pt, pt + nbElts, ptOut);
 193     return ret.retn();
 194   }
 195   vtkLongArray* d1(vtkLongArray::SafeDownCast(data));
 196   if (d1)
 197   {
 198     const long* pt(d1->GetPointer(0));
 199     std::copy(pt, pt + nbElts, ptOut);
 200     return ret.retn();
 201   }
 202   vtkIdTypeArray* d2(vtkIdTypeArray::SafeDownCast(data));
 203   if (d2)
 204   {
 205     const vtkIdType* pt(d2->GetPointer(0));
 206     std::copy(pt, pt + nbElts, ptOut);
 207     return ret.retn();
 208   }
 209   std::ostringstream oss;
 210   oss << "ConvertVTKArrayToMCArrayInt : unrecognized array \"" << typeid(*data).name()
 211       << "\" type !";
 212   throw INTERP_KERNEL::Exception(oss.str());
 213 }
 214
 215 DataArrayDouble* ConvertVTKArrayToMCArrayDouble(vtkDataArray* data)
 216 {
 217   if (!data)
 218     throw INTERP_KERNEL::Exception("ConvertVTKArrayToMCArrayDouble : internal error !");
 219   int nbTuples(data->GetNumberOfTuples()), nbComp(data->GetNumberOfComponents());
 220   std::size_t nbElts(nbTuples * nbComp);
 221   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
 222   ret->alloc(nbTuples, nbComp);
 223   for (int i = 0; i < nbComp; i++)
 224   {
 225     const char* comp(data->GetComponentName(i));
 226     if (comp)
 227       ret->setInfoOnComponent(i, comp);
 228   }
 229   double* ptOut(ret->getPointer());
 230   vtkFloatArray* d0(vtkFloatArray::SafeDownCast(data));
 231   if (d0)
 232   {
 233     const float* pt(d0->GetPointer(0));
 234     for (std::size_t i = 0; i < nbElts; i++)
 235       ptOut[i] = pt[i];
 236     return ret.retn();
 237   }
 238   vtkDoubleArray* d1(vtkDoubleArray::SafeDownCast(data));
 239   if (d1)
 240   {
 241     const double* pt(d1->GetPointer(0));
 242     std::copy(pt, pt + nbElts, ptOut);
 243     return ret.retn();
 244   }
 245   std::ostringstream oss;
 246   oss << "ConvertVTKArrayToMCArrayDouble : unrecognized array \"" << typeid(*data).name()
 247       << "\" type !";
 248   throw INTERP_KERNEL::Exception(oss.str());
 249 }
 250
 251 DataArray* ConvertVTKArrayToMCArray(vtkDataArray* data)
 252 {
 253   if (!data)
 254     throw INTERP_KERNEL::Exception("ConvertVTKArrayToMCArray : internal error !");
 255   vtkFloatArray* d0(vtkFloatArray::SafeDownCast(data));
 256   vtkDoubleArray* d1(vtkDoubleArray::SafeDownCast(data));
 257   if (d0 || d1)
 258     return ConvertVTKArrayToMCArrayDouble(data);
 259   vtkIntArray* d2(vtkIntArray::SafeDownCast(data));
 260   vtkLongArray* d3(vtkLongArray::SafeDownCast(data));
 261   if (d2 || d3)
 262     return ConvertVTKArrayToMCArrayInt(data);
 263   std::ostringstream oss;
 264   oss << "ConvertVTKArrayToMCArray : unrecognized array \"" << typeid(*data).name() << "\" type !";
 265   throw INTERP_KERNEL::Exception(oss.str());
 266 }
 267
 268 DataArrayDouble* BuildCoordsFrom(vtkPointSet* ds)
 269 {
 270   if (!ds)
 271     throw INTERP_KERNEL::Exception("BuildCoordsFrom : internal error !");
 272   vtkPoints* pts(ds->GetPoints());
 273   if (!pts)
 274     throw INTERP_KERNEL::Exception("BuildCoordsFrom : internal error 2 !");
 275   vtkDataArray* data(pts->GetData());
 276   if (!data)
 277     throw INTERP_KERNEL::Exception("BuildCoordsFrom : internal error 3 !");
 278   MCAuto<DataArrayDouble> coords(ConvertVTKArrayToMCArrayDouble(data));
 279   return coords.retn();
 280 }
 281
 282 void ConvertFromUnstructuredGrid(vtkUnstructuredGrid* ds,
 283   std::vector<MCAuto<MEDCouplingUMesh> >& ms, std::vector<MCAuto<DataArrayIdType> >& ids)
 284 {
 285   MCAuto<DataArrayDouble> coords(BuildCoordsFrom(ds));
 286   vtkIdType nbCells(ds->GetNumberOfCells());
 287   vtkUnsignedCharArray* ct(ds->GetCellTypesArray());
 288   if (!ct)
 289     throw INTERP_KERNEL::Exception("ConvertFromUnstructuredGrid : internal error");
 290   const unsigned char* ctPtr(ct->GetPointer(0));
 291   std::map<int, int> m(ComputeMapOfType());
 292   MCAuto<DataArrayIdType> lev(DataArrayIdType::New());
 293   lev->alloc(nbCells, 1);
 294   mcIdType* levPtr(lev->getPointer());
 295   for (vtkIdType i = 0; i < nbCells; i++)
 296   {
 297     std::map<int, int>::iterator it(m.find(ctPtr[i]));
 298     if (it != m.end())
 299     {
 300       const INTERP_KERNEL::CellModel& cm(
 301         INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)(*it).second));
 302       levPtr[i] = cm.getDimension();
 303     }
 304     else
 305     {
 306       std::ostringstream oss;
 307       oss << "ConvertFromUnstructuredGrid : at pos #" << i
 308           << " unrecognized VTK cell with type =" << ctPtr[i];
 309       throw INTERP_KERNEL::Exception(oss.str());
 310     }
 311   }
 312   MCAuto<DataArrayIdType> levs(lev->getDifferentValues());
 313   vtkIdTypeArray *faces(ds->GetFaces()), *faceLoc(ds->GetFaceLocations());
 314   for (const mcIdType* curLev = levs->begin(); curLev != levs->end(); curLev++)
 315   {
 316     MCAuto<MEDCouplingUMesh> m0(MEDCouplingUMesh::New("", *curLev));
 317     m0->setCoords(coords);
 318     m0->allocateCells();
 319     MCAuto<DataArrayIdType> cellIdsCurLev(lev->findIdsEqual(*curLev));
 320     for (const mcIdType* cellId = cellIdsCurLev->begin(); cellId != cellIdsCurLev->end(); cellId++)
 321     {
 322       std::map<int, int>::iterator it(m.find(ctPtr[*cellId]));
 323       vtkIdType sz;
 324       vtkIdType const* pts;
 325       ds->GetCellPoints(*cellId, sz, pts);
 326       INTERP_KERNEL::NormalizedCellType ct((INTERP_KERNEL::NormalizedCellType)(*it).second);
 327       if (ct != INTERP_KERNEL::NORM_POLYHED)
 328       {
 329         std::vector<mcIdType> conn2(sz);
 330         for (int kk = 0; kk < sz; kk++)
 331           conn2[kk] = pts[kk];
 332         m0->insertNextCell(ct, sz, &conn2[0]);
 333       }
 334       else
 335       {
 336         if (!faces || !faceLoc)
 337           throw INTERP_KERNEL::Exception(
 338             "ConvertFromUnstructuredGrid : faces are expected when there are polyhedra !");
 339         const vtkIdType *facPtr(faces->GetPointer(0)), *facLocPtr(faceLoc->GetPointer(0));
 340         std::vector<mcIdType> conn;
 341         int off0(facLocPtr[*cellId]);
 342         int nbOfFaces(facPtr[off0++]);
 343         for (int k = 0; k < nbOfFaces; k++)
 344         {
 345           int nbOfNodesInFace(facPtr[off0++]);
 346           std::copy(facPtr + off0, facPtr + off0 + nbOfNodesInFace, std::back_inserter(conn));
 347           off0 += nbOfNodesInFace;
 348           if (k < nbOfFaces - 1)
 349             conn.push_back(-1);
 350         }
 351         m0->insertNextCell(ct, conn.size(), conn.data());
 352       }
 353     }
 354     ms.push_back(m0);
 355     ids.push_back(cellIdsCurLev);
 356   }
 357 }
 358
 359 vtkSmartPointer<vtkUnstructuredGrid> ConvertUMeshFromMCToVTK(const MEDCouplingUMesh* mVor)
 360 {
 361   std::map<int, int> zeMapRev(ComputeRevMapOfType());
 362   int nbCells(mVor->getNumberOfCells());
 363   vtkSmartPointer<vtkUnstructuredGrid> ret(vtkSmartPointer<vtkUnstructuredGrid>::New());
 364   ret->Initialize();
 365   ret->Allocate();
 366   vtkSmartPointer<vtkPoints> points(vtkSmartPointer<vtkPoints>::New());
 367   {
 368     const DataArrayDouble* vorCoords(mVor->getCoords());
 369     vtkSmartPointer<vtkDoubleArray> da(vtkSmartPointer<vtkDoubleArray>::New());
 370     da->SetNumberOfComponents(vorCoords->getNumberOfComponents());
 371     da->SetNumberOfTuples(vorCoords->getNumberOfTuples());
 372     std::copy(vorCoords->begin(), vorCoords->end(), da->GetPointer(0));
 373     points->SetData(da);
 374   }
 375   mVor->checkConsistencyLight();
 376   switch (mVor->getMeshDimension())
 377   {
 378     case 3:
 379     {
 380       vtkIdType *cPtr(nullptr), *dPtr(nullptr);
 381       unsigned char* aPtr(nullptr);
 382       vtkSmartPointer<vtkUnsignedCharArray> cellTypes(vtkSmartPointer<vtkUnsignedCharArray>::New());
 383       {
 384         cellTypes->SetNumberOfComponents(1);
 385         cellTypes->SetNumberOfTuples(nbCells);
 386         aPtr = cellTypes->GetPointer(0);
 387       }
 388       vtkSmartPointer<vtkIdTypeArray> cellLocations(vtkSmartPointer<vtkIdTypeArray>::New());
 389       {
 390         cellLocations->SetNumberOfComponents(1);
 391         cellLocations->SetNumberOfTuples(nbCells);
 392         cPtr = cellLocations->GetPointer(0);
 393       }
 394       vtkSmartPointer<vtkIdTypeArray> cells(vtkSmartPointer<vtkIdTypeArray>::New());
 395       {
 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);
 400       }
 401       const mcIdType *connPtr(mVor->getNodalConnectivity()->begin()),
 402         *connIPtr(mVor->getNodalConnectivityIndex()->begin());
 403       int k(0), kk(0);
 404       std::vector<vtkIdType> ee, ff;
 405       for (int i = 0; i < nbCells; i++, connIPtr++)
 406       {
 407         INTERP_KERNEL::NormalizedCellType ct(
 408           static_cast<INTERP_KERNEL::NormalizedCellType>(connPtr[connIPtr[0]]));
 409         *aPtr++ = zeMapRev[connPtr[connIPtr[0]]];
 410         if (ct != INTERP_KERNEL::NORM_POLYHED)
 411         {
 412           int sz(connIPtr[1] - connIPtr[0] - 1);
 413           *dPtr++ = sz;
 414           dPtr = std::copy(connPtr + connIPtr[0] + 1, connPtr + connIPtr[1], dPtr);
 415           *cPtr++ = k;
 416           k += sz + 1;
 417           ee.push_back(kk);
 418         }
 419         else
 420         {
 421           std::set<mcIdType> s(connPtr + connIPtr[0] + 1, connPtr + connIPtr[1]);
 422           s.erase(-1);
 423           int nbFace(std::count(connPtr + connIPtr[0] + 1, connPtr + connIPtr[1], -1) + 1);
 424           ff.push_back(nbFace);
 425           const mcIdType* work(connPtr + connIPtr[0] + 1);
 426           for (int j = 0; j < nbFace; j++)
 427           {
 428             const mcIdType* work2 = std::find(work, connPtr + connIPtr[1], -1);
 429             ff.push_back(std::distance(work, work2));
 430             ff.insert(ff.end(), work, work2);
 431             work = work2 + 1;
 432           }
 433           *dPtr++ = (int)s.size();
 434           dPtr = std::copy(s.begin(), s.end(), dPtr);
 435           *cPtr++ = k;
 436           k += (int)s.size() + 1;
 437           ee.push_back(kk);
 438           kk += connIPtr[1] - connIPtr[0] + 1;
 439         }
 440       }
 441       //
 442       vtkSmartPointer<vtkIdTypeArray> faceLocations(vtkSmartPointer<vtkIdTypeArray>::New());
 443       {
 444         faceLocations->SetNumberOfComponents(1);
 445         faceLocations->SetNumberOfTuples(ee.size());
 446         std::copy(ee.begin(), ee.end(), faceLocations->GetPointer(0));
 447       }
 448       vtkSmartPointer<vtkIdTypeArray> faces(vtkSmartPointer<vtkIdTypeArray>::New());
 449       {
 450         faces->SetNumberOfComponents(1);
 451         faces->SetNumberOfTuples(ff.size());
 452         std::copy(ff.begin(), ff.end(), faces->GetPointer(0));
 453       }
 454       vtkSmartPointer<vtkCellArray> cells2(vtkSmartPointer<vtkCellArray>::New());
 455       cells2->SetCells(nbCells, cells);
 456       ret->SetCells(cellTypes, cellLocations, cells2, faceLocations, faces);
 457       break;
 458     }
 459     case 2:
 460     {
 461       vtkSmartPointer<vtkUnsignedCharArray> cellTypes(vtkSmartPointer<vtkUnsignedCharArray>::New());
 462       {
 463         cellTypes->SetNumberOfComponents(1);
 464         cellTypes->SetNumberOfTuples(nbCells);
 465         unsigned char* ptr(cellTypes->GetPointer(0));
 466         std::fill(ptr, ptr + nbCells, zeMapRev[(int)INTERP_KERNEL::NORM_POLYGON]);
 467       }
 468       vtkIdType *cPtr(nullptr), *dPtr(nullptr);
 469       vtkSmartPointer<vtkIdTypeArray> cellLocations(vtkSmartPointer<vtkIdTypeArray>::New());
 470       {
 471         cellLocations->SetNumberOfComponents(1);
 472         cellLocations->SetNumberOfTuples(nbCells);
 473         cPtr = cellLocations->GetPointer(0);
 474       }
 475       vtkSmartPointer<vtkIdTypeArray> cells(vtkSmartPointer<vtkIdTypeArray>::New());
 476       {
 477         cells->SetNumberOfComponents(1);
 478         cells->SetNumberOfTuples(mVor->getNodalConnectivity()->getNumberOfTuples());
 479         dPtr = cells->GetPointer(0);
 480       }
 481       const mcIdType *connPtr(mVor->getNodalConnectivity()->begin()),
 482         *connIPtr(mVor->getNodalConnectivityIndex()->begin());
 483       mcIdType k(0);
 484       for (mcIdType i = 0; i < nbCells; i++, connIPtr++)
 485       {
 486         *dPtr++ = connIPtr[1] - connIPtr[0] - 1;
 487         dPtr = std::copy(connPtr + connIPtr[0] + 1, connPtr + connIPtr[1], dPtr);
 488         *cPtr++ = k;
 489         k += connIPtr[1] - connIPtr[0];
 490       }
 491       vtkSmartPointer<vtkCellArray> cells2(vtkSmartPointer<vtkCellArray>::New());
 492       cells2->SetCells(nbCells, cells);
 493       ret->SetCells(cellTypes, cellLocations, cells2);
 494       break;
 495     }
 496     case 1:
 497     {
 498       vtkSmartPointer<vtkUnsignedCharArray> cellTypes(vtkSmartPointer<vtkUnsignedCharArray>::New());
 499       {
 500         cellTypes->SetNumberOfComponents(1);
 501         cellTypes->SetNumberOfTuples(nbCells);
 502         unsigned char* ptr(cellTypes->GetPointer(0));
 503         std::fill(ptr, ptr + nbCells, zeMapRev[(int)INTERP_KERNEL::NORM_SEG2]);
 504       }
 505       vtkIdType *cPtr(nullptr), *dPtr(nullptr);
 506       vtkSmartPointer<vtkIdTypeArray> cellLocations(vtkSmartPointer<vtkIdTypeArray>::New());
 507       {
 508         cellLocations->SetNumberOfComponents(1);
 509         cellLocations->SetNumberOfTuples(nbCells);
 510         cPtr = cellLocations->GetPointer(0);
 511       }
 512       vtkSmartPointer<vtkIdTypeArray> cells(vtkSmartPointer<vtkIdTypeArray>::New());
 513       {
 514         cells->SetNumberOfComponents(1);
 515         cells->SetNumberOfTuples(mVor->getNodalConnectivity()->getNumberOfTuples());
 516         dPtr = cells->GetPointer(0);
 517       }
 518       const mcIdType *connPtr(mVor->getNodalConnectivity()->begin()),
 519         *connIPtr(mVor->getNodalConnectivityIndex()->begin());
 520       for (int i = 0; i < nbCells; i++, connIPtr++)
 521       {
 522         *dPtr++ = 2;
 523         dPtr = std::copy(connPtr + connIPtr[0] + 1, connPtr + connIPtr[1], dPtr);
 524         *cPtr++ = 3 * i;
 525       }
 526       vtkSmartPointer<vtkCellArray> cells2(vtkSmartPointer<vtkCellArray>::New());
 527       cells2->SetCells(nbCells, cells);
 528       ret->SetCells(cellTypes, cellLocations, cells2);
 529       break;
 530     }
 531     default:
 532       throw INTERP_KERNEL::Exception("Not implemented yet !");
 533   }
 534   ret->SetPoints(points);
 535   return ret;
 536 }
 537
 538 MCAuto<DataArrayDouble> ForceBuilder(const std::vector<std::size_t>& TAB, const DataArrayDouble* matrix, const DataArrayDouble* eqn)
 539 {
 540   MCAuto<DataArrayDouble> tmp0, tmp1, ret;
 541   tmp0 = matrix->keepSelectedComponents({ TAB[0] });
 542   tmp1 = eqn->keepSelectedComponents({ 0 });
 543   MCAuto<DataArrayDouble> p0(DataArrayDouble::Multiply(tmp0, tmp1));
 544   tmp0 = matrix->keepSelectedComponents({ TAB[1] });
 545   tmp1 = eqn->keepSelectedComponents({ 1 });
 546   MCAuto<DataArrayDouble> p1(DataArrayDouble::Multiply(tmp0, tmp1));
 547   ret = DataArrayDouble::Add(p0, p1);
 548   tmp0 = matrix->keepSelectedComponents({ TAB[2] });
 549   tmp1 = eqn->keepSelectedComponents({ 2 });
 550   MCAuto<DataArrayDouble> p2(DataArrayDouble::Multiply(tmp0, tmp1));
 551   ret = DataArrayDouble::Add(ret, p2);
 552   return ret;
 553 }
 554
 555 double ReturnInertia(
 556   const double pOut[3], const DataArrayDouble* OM, const DataArrayDouble* area_field_ids)
 557 {
 558   MCAuto<DataArrayDouble> base_X(DataArrayDouble::New());
 559   base_X->alloc(OM->getNumberOfTuples(), 3);
 560   base_X->setPartOfValuesSimple1(pOut[0], 0, OM->getNumberOfTuples(), 1, 0, 1, 1);
 561   base_X->setPartOfValuesSimple1(pOut[1], 0, OM->getNumberOfTuples(), 1, 1, 2, 1);
 562   base_X->setPartOfValuesSimple1(pOut[2], 0, OM->getNumberOfTuples(), 1, 2, 3, 1);
 563   MCAuto<DataArrayDouble> dist_to_base_X;
 564   {
 565     MCAuto<DataArrayDouble> tmp(DataArrayDouble::Dot(OM, base_X));
 566     tmp = DataArrayDouble::Multiply(tmp, base_X);
 567     tmp = DataArrayDouble::Substract(OM, tmp);
 568     dist_to_base_X = tmp->magnitude();
 569   }
 570   MCAuto<DataArrayDouble> inertiaArr;
 571   {
 572     MCAuto<DataArrayDouble> tmp(DataArrayDouble::Multiply(dist_to_base_X, dist_to_base_X));
 573     inertiaArr = DataArrayDouble::Multiply(tmp, area_field_ids);
 574   }
 575   double inertiaTmp;
 576   inertiaArr->accumulate(&inertiaTmp);
 577   return inertiaTmp;
 578 }
 579
 580 void FindPrincipalAxeInternal(const double startVector[3], const double normalFace[3],
 581   const DataArrayDouble* OM, const DataArrayDouble* area_field_ids,
 582   const std::vector<double>& posToIterate, double& angleDegree, double outputAxis[3],
 583   double& inertia)
 584 {
 585   constexpr double CENTER[3] = { 0, 0, 0 };
 586   inertia = -std::numeric_limits<double>::max();
 587   for (auto zePos : posToIterate)
 588   {
 589     double p[3] = { startVector[0], startVector[1], startVector[2] }, pOut[3];
 590     DataArrayDouble::Rotate3DAlg(CENTER, normalFace, zePos / 180. * M_PI, 1, p, pOut);
 591     double inertiaTmp = ReturnInertia(pOut, OM, area_field_ids);
 592     if (inertiaTmp > inertia)
 593     {
 594       inertia = inertiaTmp;
 595       std::copy(pOut, pOut + 3, outputAxis);
 596       angleDegree = zePos;
 597     }
 598   }
 599 }
 600
 601 void FindPrincipalAxe(const double startVector[3], const double normalFace[3],
 602   const DataArrayDouble* OM, const DataArrayDouble* area_field_ids, double& angleDegree,
 603   double outputAxis[3], double& inertia)
 604 {
 605   std::vector<double> R(
 606     { 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170 });
 607   FindPrincipalAxeInternal(
 608     startVector, normalFace, OM, area_field_ids, R, angleDegree, outputAxis, inertia);
 609   for (int i = 0; i < 5; ++i)
 610   {
 611     std::vector<double> Q({ -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 });
 612     const double CST(std::pow((double)10., (double)-i));
 613     std::for_each(Q.begin(), Q.end(), [angleDegree, CST](double& v) { v = CST * v + angleDegree; });
 614     FindPrincipalAxeInternal(
 615       startVector, normalFace, OM, area_field_ids, Q, angleDegree, outputAxis, inertia);
 616   }
 617 }
 618
 619 vtkSmartPointer<vtkTable> ComputeTorseurCIH(vtkUnstructuredGrid* usgIn)
 620 {
 621   std::vector<MCAuto<MEDCouplingUMesh> > m;
 622   {
 623     std::vector<MCAuto<DataArrayIdType> > ids;
 624     ConvertFromUnstructuredGrid(usgIn, m, ids);
 625   }
 626   vtkDataArray* sief(nullptr);
 627   {
 628     int nArrays(usgIn->GetPointData()->GetNumberOfArrays());
 629     for (int i = 0; i < nArrays; i++)
 630     {
 631       vtkDataArray* array(usgIn->GetPointData()->GetArray(i));
 632       if (!array)
 633         continue;
 634       std::string name(array->GetName());
 635       if (name.find("SIEF") != std::string::npos)
 636         if (array->GetNumberOfComponents() == 6)
 637         {
 638           if (sief)
 639           {
 640             std::ostringstream oss;
 641             oss << "ComputeTorseurCIH : several candidates for SIEF field !";
 642             throw INTERP_KERNEL::Exception(oss.str());
 643           }
 644           sief = array;
 645         }
 646     }
 647   }
 648   if (!sief)
 649     throw INTERP_KERNEL::Exception("ComputeTorseurCIH : unable to find a field for SIEF!");
 650   MCAuto<MEDCouplingFieldDouble> area_field(m[0]->getMeasureField(true));
 651   double area;
 652   area_field->accumulate(&area); // 1
 653   MCAuto<DataArrayDouble> centerOfMassField(m[0]->computeCellCenterOfMass());
 654   double centerOfMass[3];
 655   {
 656     MCAuto<DataArrayDouble> tmp(
 657       DataArrayDouble::Multiply(centerOfMassField, area_field->getArray()));
 658     tmp->accumulate(centerOfMass);
 659     std::for_each(centerOfMass, centerOfMass + 3, [area](double& v) { v /= area; });
 660   } // 2
 661   m[0]->unPolyze();
 662   std::set<INTERP_KERNEL::NormalizedCellType> s(m[0]->getAllGeoTypes());
 663   if (s.size() != 1)
 664   {
 665     std::ostringstream oss;
 666     oss << "ComputeTorseurCIH : Only TRI3 are supported !";
 667     throw INTERP_KERNEL::Exception(oss.str());
 668   }
 669   if (*(s.begin()) != INTERP_KERNEL::NORM_TRI3)
 670   {
 671     std::ostringstream oss;
 672     oss << "ComputeTorseurCIH : Only TRI3 are supported !";
 673     throw INTERP_KERNEL::Exception(oss.str());
 674   }
 675   MCAuto<MEDCouplingFieldDouble> f(MEDCouplingFieldDouble::New(MEDCoupling::ON_NODES));
 676   {
 677     f->setMesh(m[0]);
 678     MCAuto<DataArrayDouble> tmp(ConvertVTKArrayToMCArrayDouble(sief));
 679     f->setArray(tmp);
 680   }
 681   MCAuto<MEDCouplingFieldDouble> fCell(f->nodeToCellDiscretization());
 682   MCAuto<DataArrayIdType> ids;
 683   {
 684     MCAuto<DataArrayIdType> tmp(area_field->getArray()->findIdsLowerThan(1e-7));
 685     ids = tmp->buildComplement(m[0]->getNumberOfCells());
 686   }
 687   MCAuto<MEDCouplingUMesh> m_ids(m[0]->buildPartOfMySelf(ids->begin(), ids->end()));
 688   MCAuto<DataArrayDouble> eqn;
 689   {
 690     MCAuto<DataArrayDouble> tmp(m_ids->computePlaneEquationOf3DFaces());
 691     eqn = tmp->keepSelectedComponents({ 0, 1, 2 });
 692     tmp = eqn->magnitude();
 693     eqn = DataArrayDouble::Divide(eqn, tmp);
 694   }
 695   MCAuto<DataArrayDouble> area_field_ids(
 696     area_field->getArray()->selectByTupleId(ids->begin(), ids->end()));
 697   MCAuto<DataArrayDouble> area_vector(DataArrayDouble::Multiply(eqn, area_field_ids));
 698   MCAuto<DataArrayDouble> matrix(fCell->getArray()->selectByTupleId(ids->begin(), ids->end()));
 699   MCAuto<DataArrayDouble> F_x, F_y, F_z;
 700   {
 701     F_x = ForceBuilder({ 0, 3, 4 }, matrix, eqn);
 702     F_y = ForceBuilder({ 3, 1, 5 }, matrix, eqn);
 703     F_z = ForceBuilder({ 4, 5, 2 }, matrix, eqn);
 704   }
 705   //
 706   //
 707   MCAuto<DataArrayDouble> F(DataArrayDouble::Meld({ F_x, F_y, F_z }));
 708   F->multiplyEqual(area_vector);
 709   double ZeForce[3], normalFace[3];
 710   F->accumulate(ZeForce);
 711   eqn->accumulate(normalFace);
 712   {
 713     double normalFaceNorm(sqrt(normalFace[0] * normalFace[0] + normalFace[1] * normalFace[1] +
 714       normalFace[2] * normalFace[2]));
 715     std::for_each(normalFace, normalFace + 3, [normalFaceNorm](double& v) { v /= normalFaceNorm; });
 716   }
 717   double ForceNormale[3]; // 3
 718   {
 719     double tmp(
 720       ZeForce[0] * normalFace[0] + ZeForce[1] * normalFace[1] + ZeForce[2] * normalFace[2]);
 721     ForceNormale[0] = tmp * normalFace[0];
 722     ForceNormale[1] = tmp * normalFace[1];
 723     ForceNormale[2] = tmp * normalFace[2];
 724   }
 725   double TangentForce[3] = { ZeForce[0] - ForceNormale[0], ZeForce[1] - ForceNormale[1],
 726     ZeForce[2] - ForceNormale[2] }; // 4
 727   MCAuto<DataArrayDouble> bary(m_ids->computeCellCenterOfMass());
 728   MCAuto<DataArrayDouble> OM;
 729   {
 730     MCAuto<DataArrayDouble> centerOfMass2(DataArrayDouble::New());
 731     centerOfMass2->alloc(1, 3);
 732     std::copy(centerOfMass, centerOfMass + 3, centerOfMass2->getPointer());
 733     OM = DataArrayDouble::Substract(bary, centerOfMass2);
 734   }
 735   double momentum[3]; // 5
 736   {
 737     MCAuto<DataArrayDouble> tmp(DataArrayDouble::CrossProduct(OM, F));
 738     tmp->accumulate(momentum);
 739   }
 740   double InertiaNormale(ReturnInertia(normalFace, OM, area_field_ids)); // 6
 741   double base[9];
 742   DataArrayDouble::GiveBaseForPlane(normalFace, base);
 743   double angleDegree, outputAxis[3], inertia;
 744   FindPrincipalAxe(base, normalFace, OM, area_field_ids, angleDegree, outputAxis, inertia);
 745   double tangentOther[3] = { normalFace[1] * outputAxis[2] - normalFace[2] * outputAxis[1],
 746     normalFace[2] * outputAxis[0] - normalFace[0] * outputAxis[2],
 747     normalFace[0] * outputAxis[1] - normalFace[1] * outputAxis[0] };
 748   double inertiaOther(ReturnInertia(tangentOther, OM, area_field_ids));
 749   vtkSmartPointer<vtkTable> ret(vtkSmartPointer<vtkTable>::New());
 750   vtkSmartPointer<vtkStringArray> col0(vtkSmartPointer<vtkStringArray>::New());
 751   constexpr int NB_ROWS = 11;
 752   col0->SetNumberOfComponents(1);
 753   col0->SetNumberOfTuples(NB_ROWS);
 754   col0->SetName("Grandeur");
 755   // scalaire
 756   col0->SetValue(0, strdup("Aire"));
 757   col0->SetValue(1, strdup("Inertie Normal"));
 758   col0->SetValue(2, strdup("Inertie Tangentielle principale"));
 759   col0->SetValue(3, strdup("Inertie Tangentielle secondaire"));
 760   // vectoriel
 761   col0->SetValue(4, strdup("Position du centre de gravite"));
 762   col0->SetValue(5, strdup("Effort Normal"));
 763   col0->SetValue(6, strdup("Effort Tangentiel"));
 764   col0->SetValue(7, strdup("Axe Normal"));
 765   col0->SetValue(8, strdup("Axe Tangentiel principal"));
 766   col0->SetValue(9, strdup("Axe Tangentiel secondaire"));
 767   col0->SetValue(10, strdup("Moment au centre de gravite"));
 768   ret->AddColumn(col0);
 769   //
 770   vtkSmartPointer<vtkDoubleArray> col1(vtkSmartPointer<vtkDoubleArray>::New());
 771   col1->SetName("X");
 772   col1->SetNumberOfComponents(1);
 773   col1->SetNumberOfTuples(NB_ROWS);
 774   col1->SetValue(0, area);
 775   col1->SetValue(1, InertiaNormale);
 776   col1->SetValue(2, inertia);
 777   col1->SetValue(3, inertiaOther);
 778   col1->SetValue(4, centerOfMass[0]);
 779   col1->SetValue(5, ForceNormale[0]);
 780   col1->SetValue(6, TangentForce[0]);
 781   col1->SetValue(7, normalFace[0]);
 782   col1->SetValue(8, outputAxis[0]);
 783   col1->SetValue(9, tangentOther[0]);
 784   col1->SetValue(10, momentum[0]);
 785   ret->AddColumn(col1);
 786   //
 787   vtkSmartPointer<vtkDoubleArray> col2(vtkSmartPointer<vtkDoubleArray>::New());
 788   col2->SetName("Y");
 789   col2->SetNumberOfComponents(1);
 790   col2->SetNumberOfTuples(NB_ROWS);
 791   col2->SetValue(0, 0.);
 792   col2->SetValue(1, 0.);
 793   col2->SetValue(2, 0.);
 794   col2->SetValue(3, 0.);
 795   col2->SetValue(4, centerOfMass[1]);
 796   col2->SetValue(5, ForceNormale[1]);
 797   col2->SetValue(6, TangentForce[1]);
 798   col2->SetValue(7, normalFace[1]);
 799   col2->SetValue(8, outputAxis[1]);
 800   col2->SetValue(9, tangentOther[1]);
 801   col2->SetValue(10, momentum[1]);
 802   ret->AddColumn(col2);
 803   //
 804   vtkSmartPointer<vtkDoubleArray> col3(vtkSmartPointer<vtkDoubleArray>::New());
 805   col3->SetName("Z");
 806   col3->SetNumberOfComponents(1);
 807   col3->SetNumberOfTuples(NB_ROWS);
 808   col3->SetValue(0, 0.);
 809   col3->SetValue(1, 0.);
 810   col3->SetValue(2, 0.);
 811   col3->SetValue(3, 0.);
 812   col3->SetValue(4, centerOfMass[2]);
 813   col3->SetValue(5, ForceNormale[2]);
 814   col3->SetValue(6, TangentForce[2]);
 815   col3->SetValue(7, normalFace[2]);
 816   col3->SetValue(8, outputAxis[2]);
 817   col3->SetValue(9, tangentOther[2]);
 818   col3->SetValue(10, momentum[2]);
 819   ret->AddColumn(col3);
 820   return ret;
 821 }
 822
 823 ////////////////////
 824
 825 int vtkTorseurCIH::FillOutputPortInformation(int vtkNotUsed(port), vtkInformation* info)
 826 {
 827   info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkTable");
 828   return 1;
 829 }
 830
 831 int vtkTorseurCIH::RequestInformation(
 832   vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
 833 {
 834   // std::cerr << "########################################## vtkTorseurCIH::RequestInformation
 835   // ##########################################" << std::endl;
 836   try
 837   {
 838     vtkSmartPointer<vtkUnstructuredGrid> usgIn;
 839     ExtractInfo(inputVector[0], usgIn);
 840   }
 841   catch (INTERP_KERNEL::Exception& e)
 842   {
 843     std::ostringstream oss;
 844     oss << "Exception has been thrown in vtkTorseurCIH::RequestInformation : " << e.what()
 845         << std::endl;
 846     if (this->HasObserver("ErrorEvent"))
 847       this->InvokeEvent("ErrorEvent", const_cast<char*>(oss.str().c_str()));
 848     else
 849       vtkOutputWindowDisplayErrorText(const_cast<char*>(oss.str().c_str()));
 850     vtkObject::BreakOnError();
 851     return 0;
 852   }
 853   return 1;
 854 }
 855
 856 int vtkTorseurCIH::RequestData(
 857   vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
 858 {
 859   // std::cerr << "########################################## vtkTorseurCIH::RequestData
 860   // ##########################################" << std::endl;
 861   try
 862   {
 863     vtkSmartPointer<vtkUnstructuredGrid> usgIn;
 864     ExtractInfo(inputVector[0], usgIn);
 865     //
 866     vtkSmartPointer<vtkTable> ret(ComputeTorseurCIH(usgIn));
 867     vtkInformation* inInfo(inputVector[0]->GetInformationObject(0));
 868     vtkInformation* outInfo(outputVector->GetInformationObject(0));
 869     vtkTable* output(vtkTable::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())));
 870     output->ShallowCopy(ret);
 871   }
 872   catch (INTERP_KERNEL::Exception& e)
 873   {
 874     std::ostringstream oss;
 875     oss << "Exception has been thrown in vtkTorseurCIH::RequestData : " << e.what() << std::endl;
 876     if (this->HasObserver("ErrorEvent"))
 877       this->InvokeEvent("ErrorEvent", const_cast<char*>(oss.str().c_str()));
 878     else
 879       vtkOutputWindowDisplayErrorText(const_cast<char*>(oss.str().c_str()));
 880     vtkObject::BreakOnError();
 881     return 0;
 882   }
 883   return 1;
 884 }
 885
 886 void vtkTorseurCIH::PrintSelf(ostream& os, vtkIndent indent)
 887 {
 888   this->Superclass::PrintSelf(os, indent);
 889 }