src/Plugins/SimpleMode/plugin/SimpleModeModule/vtkSimpleMode.cxx

   1 // Copyright (C) 2017-2023  CEA, EDF
   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 "vtkSimpleMode.h"
  22
  23 #include <vtkAdjacentVertexIterator.h>
  24 #include <vtkAlgorithmOutput.h>
  25 #include <vtkCellData.h>
  26 #include <vtkCharArray.h>
  27 #include <vtkMergeBlocks.h>
  28 #include <vtkDataArraySelection.h>
  29 #include <vtkDataObjectTreeIterator.h>
  30 #include <vtkDataSet.h>
  31 #include <vtkDataSetAttributes.h>
  32 #include <vtkDataSetSurfaceFilter.h>
  33 #include <vtkDemandDrivenPipeline.h>
  34 #include <vtkDoubleArray.h>
  35 #include <vtkExecutive.h>
  36 #include <vtkFloatArray.h>
  37 #include <vtkInEdgeIterator.h>
  38 #include <vtkInformation.h>
  39 #include <vtkInformationDataObjectKey.h>
  40 #include <vtkInformationStringKey.h>
  41 #include <vtkInformationVector.h>
  42 #include <vtkIntArray.h>
  43 #include <vtkMultiBlockDataGroupFilter.h>
  44 #include <vtkMultiBlockDataSet.h>
  45 #include <vtkMutableDirectedGraph.h>
  46 #include <vtkObjectFactory.h>
  47 #include <vtkPointData.h>
  48 #include <vtkStreamingDemandDrivenPipeline.h>
  49 #include <vtkStringArray.h>
  50 #include <vtkTimeStamp.h>
  51 #include <vtkUnsignedCharArray.h>
  52 #include <vtkUnstructuredGrid.h>
  53 #include <vtkVariantArray.h>
  54 #include <vtkWarpScalar.h>
  55 #include <vtkWarpVector.h>
  56
  57 #ifdef WIN32
  58 #define _USE_MATH_DEFINES
  59 #include <functional>
  60 #include <math.h>
  61 #endif
  62
  63 #include <deque>
  64 #include <map>
  65 #include <sstream>
  66
  67 vtkStandardNewMacro(vtkSimpleMode)
  68
  69 static const char ZE_DISPLACEMENT_NAME1[] = "@@ForReal?@@";
  70
  71 static const char ZE_DISPLACEMENT_NAME2[] = "@@ForImag?@@";
  72
  73 static const char ZE_DISPLACEMENT_NAME3[] = "MagnitudeOfCpxDisp";
  74
  75 static const double EPS = 1e-12;
  76
  77 ///////////////////
  78
  79 class MZCException : public std::exception
  80 {
  81 public:
  82   MZCException(const std::string& s)
  83     : _reason(s)
  84   {
  85   }
  86   virtual const char* what() const noexcept { return _reason.c_str(); }
  87   virtual ~MZCException() noexcept {}
  88
  89 private:
  90   std::string _reason;
  91 };
  92
  93 // ValueTypeT
  94
  95 /*template<class T>
  96 struct ArrayTraits
  97 {
  98   typedef T EltType;
  99   };*/
 100
 101 template <class VTK_ARRAY_T>
 102 vtkSmartPointer<VTK_ARRAY_T> ForceTo3CompoImpl(VTK_ARRAY_T* arr)
 103 {
 104   using ELT_TYPE = typename VTK_ARRAY_T::ValueType;
 105   if (!arr)
 106     return vtkSmartPointer<VTK_ARRAY_T>();
 107   vtkIdType nbCompo(arr->GetNumberOfComponents()), nbTuples(arr->GetNumberOfTuples());
 108   if (nbCompo == 3)
 109   {
 110     vtkSmartPointer<VTK_ARRAY_T> ret(arr);
 111     return ret;
 112   }
 113   if (nbCompo == 6)
 114   {
 115     vtkSmartPointer<VTK_ARRAY_T> ret(vtkSmartPointer<VTK_ARRAY_T>::New());
 116     ret->SetNumberOfComponents(3);
 117     ret->SetNumberOfTuples(nbTuples);
 118     const ELT_TYPE* srcPt(arr->Begin());
 119     ELT_TYPE* destPt(ret->Begin());
 120     for (vtkIdType i = 0; i < nbTuples; i++, destPt += 3, srcPt += 6)
 121       std::copy(srcPt, srcPt + 3, destPt);
 122     return ret;
 123   }
 124   throw MZCException("ForceTo3CompoImpl : internal error ! 6 or 3 compo arrays expected !");
 125 }
 126
 127 vtkSmartPointer<vtkDataArray> ForceTo3Compo(vtkDataArray* arr)
 128 {
 129   vtkDoubleArray* arr0(vtkDoubleArray::SafeDownCast(arr));
 130   if (arr0)
 131     return ForceTo3CompoImpl<vtkDoubleArray>(arr0);
 132   vtkFloatArray* arr1(vtkFloatArray::SafeDownCast(arr));
 133   if (arr1)
 134     return ForceTo3CompoImpl<vtkFloatArray>(arr1);
 135   throw MZCException("ForceTo3Compo : array is NEITHER float64 NOR float32 array !");
 136 }
 137
 138 template <class VTK_ARRAY_T>
 139 void FeedDataInternal(VTK_ARRAY_T* arrReal, double cst1, double* ptToFeed1)
 140 {
 141   vtkIdType nbTuples(arrReal->GetNumberOfTuples());
 142   using ELT_TYPE = typename VTK_ARRAY_T::ValueType;
 143   const ELT_TYPE* srcPt1(arrReal->Begin());
 144   std::for_each(srcPt1, srcPt1 + 3 * nbTuples, [&ptToFeed1, cst1](const ELT_TYPE& elt) {
 145     *ptToFeed1 = (double)elt * cst1;
 146     ptToFeed1++;
 147   });
 148 }
 149
 150 void FeedData(vtkDataArray* arr, double cst1, double* ptToFeed1)
 151 {
 152   vtkDoubleArray* arr0(vtkDoubleArray::SafeDownCast(arr));
 153   if (arr0)
 154     return FeedDataInternal<vtkDoubleArray>(arr0, cst1, ptToFeed1);
 155   vtkFloatArray* arr1(vtkFloatArray::SafeDownCast(arr));
 156   if (arr1)
 157     return FeedDataInternal<vtkFloatArray>(arr1, cst1, ptToFeed1);
 158   throw MZCException("FeedData : array is NEITHER float64 NOR float32 array !");
 159 }
 160
 161 std::vector<std::string> GetPossibleArrayNames(vtkDataSet* dataset)
 162 {
 163   if (!dataset)
 164     throw MZCException("The input dataset is null !");
 165   std::vector<std::string> ret;
 166   vtkPointData* att(dataset->GetPointData());
 167   for (int i = 0; i < att->GetNumberOfArrays(); i++)
 168   {
 169     vtkDataArray* locArr(att->GetArray(i));
 170     int nbComp(locArr->GetNumberOfComponents());
 171     if (nbComp != 3 && nbComp != 6)
 172       continue;
 173     std::string s(locArr->GetName());
 174     ret.push_back(s);
 175   }
 176   return ret;
 177 }
 178
 179 std::string FindTheBest(
 180   const std::vector<std::string>& arrNames, const std::string& key0, const std::string& key1)
 181 {
 182   std::string ret;
 183   char points(0);
 184   if (arrNames.empty())
 185     return ret;
 186   for (std::vector<std::string>::const_iterator it = arrNames.begin(); it != arrNames.end(); it++)
 187   {
 188     char curNbPts(1);
 189     if ((*it).find(key0, 0) != std::string::npos)
 190       curNbPts++;
 191     if ((*it).find(key1, 0) != std::string::npos)
 192       curNbPts++;
 193     if (curNbPts > points)
 194     {
 195       points = curNbPts;
 196       ret = *it;
 197     }
 198   }
 199   return ret;
 200 }
 201
 202 std::string FindBestRealAmong(const std::vector<std::string>& arrNames)
 203 {
 204   static const char KEY1[] = "DEPL";
 205   static const char KEY2[] = "REEL";
 206   return FindTheBest(arrNames, KEY1, KEY2);
 207 }
 208
 209 std::string FindBestImagAmong(const std::vector<std::string>& arrNames)
 210 {
 211   static const char KEY1[] = "DEPL";
 212   static const char KEY2[] = "IMAG";
 213   return FindTheBest(arrNames, KEY1, KEY2);
 214 }
 215
 216 vtkUnstructuredGrid* ExtractInfo1(vtkInformationVector* inputVector)
 217 {
 218   vtkInformation* inputInfo(inputVector->GetInformationObject(0));
 219   vtkDataSet* input(0);
 220   vtkDataSet* input0(vtkDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
 221   vtkMultiBlockDataSet* input1(
 222     vtkMultiBlockDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
 223   if (input0)
 224     input = input0;
 225   else
 226   {
 227     if (!input1)
 228       throw MZCException(
 229         "Input dataSet must be a DataSet or single elt multi block dataset expected !");
 230     if (input1->GetNumberOfBlocks() != 1)
 231       throw MZCException("Input dataSet is a multiblock dataset with not exactly one block ! Use "
 232                          "MergeBlocks or ExtractBlocks filter before calling this filter !");
 233     vtkDataObject* input2(input1->GetBlock(0));
 234     if (!input2)
 235       throw MZCException("Input dataSet is a multiblock dataset with exactly one block but this "
 236                          "single element is NULL !");
 237     vtkDataSet* input2c(vtkDataSet::SafeDownCast(input2));
 238     if (!input2c)
 239       throw MZCException(
 240         "Input dataSet is a multiblock dataset with exactly one block but this single element is "
 241         "not a dataset ! Use MergeBlocks or ExtractBlocks filter before calling this filter !");
 242     input = input2c;
 243   }
 244   if (!input)
 245     throw MZCException("Input data set is NULL !");
 246   vtkUnstructuredGrid* usgIn(vtkUnstructuredGrid::SafeDownCast(input));
 247   if (!usgIn)
 248     throw MZCException("Input data set is not an unstructured mesh ! This filter works only on "
 249                        "unstructured meshes !");
 250   return usgIn;
 251 }
 252
 253 void ExtractInfo3(vtkDataSet* ds, const std::string& arrName, vtkDataArray*& arr, int& idx)
 254 {
 255   vtkPointData* att(ds->GetPointData());
 256   if (!att)
 257     throw MZCException("Input dataset has no point data attribute ! Impossible to move mesh !");
 258   vtkDataArray* zeArr(0);
 259   int i(0);
 260   for (; i < att->GetNumberOfArrays(); i++)
 261   {
 262     vtkDataArray* locArr(att->GetArray(i));
 263     std::string s(locArr->GetName());
 264     if (s == arrName)
 265     {
 266       zeArr = locArr;
 267       break;
 268     }
 269   }
 270   if (!zeArr)
 271   {
 272     std::ostringstream oss;
 273     oss << "Impossible to locate the array called \"" << arrName << "\" used to move mesh !";
 274     throw MZCException(oss.str());
 275   }
 276   arr = zeArr;
 277   idx = i;
 278 }
 279
 280 void ExtractInfo2(vtkDataSet* ds, const std::string& arrName, vtkDataArray*& arr)
 281 {
 282   int dummy;
 283   ExtractInfo3(ds, arrName, arr, dummy);
 284   vtkDoubleArray* arr1(vtkDoubleArray::SafeDownCast(arr));
 285   vtkFloatArray* arr2(vtkFloatArray::SafeDownCast(arr));
 286   if (!arr1 && !arr2)
 287   {
 288     std::ostringstream oss;
 289     oss << "Array called \"" << arrName
 290         << "\" has been located but this is NEITHER float64 NOR float32 array !";
 291     throw MZCException(oss.str());
 292   }
 293   if (arr->GetNumberOfComponents() != 3 && arr->GetNumberOfComponents() != 6)
 294   {
 295     std::ostringstream oss;
 296     oss << "Float64 array called \"" << arrName
 297         << "\" has been located but this array has not exactly 3 or 6 components as it should !";
 298     throw MZCException(oss.str());
 299   }
 300   if (arr->GetNumberOfTuples() != ds->GetNumberOfPoints())
 301   {
 302     std::ostringstream oss;
 303     oss << "Float64-1 components array called \"" << arrName
 304         << "\" has been located but the number of tuples is invalid ! Should be "
 305         << ds->GetNumberOfPoints() << " instead of " << arr->GetNumberOfTuples() << " !";
 306     throw MZCException(oss.str());
 307   }
 308 }
 309
 310 ////////////////////
 311
 312 class vtkSimpleMode::vtkSimpleModeInternal
 313 {
 314 public:
 315   vtkSimpleModeInternal()
 316     : _surface(0)
 317   {
 318   }
 319   vtkPolyData* performConnection(vtkDataSet* ds);
 320   void setFieldForReal(const std::string& st) { _real = st; }
 321   std::string getFieldForReal() const { return _real; }
 322   ~vtkSimpleModeInternal();
 323
 324 private:
 325   vtkDataSetSurfaceFilter* _surface;
 326
 327 private:
 328   std::string _real;
 329 };
 330
 331 vtkPolyData* vtkSimpleMode::vtkSimpleModeInternal::performConnection(vtkDataSet* ds)
 332 {
 333   if (!_surface)
 334   {
 335     _surface = vtkDataSetSurfaceFilter::New();
 336     _surface->SetInputData(ds);
 337   }
 338   _surface->Update();
 339   return _surface->GetOutput();
 340 }
 341
 342 vtkSimpleMode::vtkSimpleModeInternal::~vtkSimpleModeInternal()
 343 {
 344   if (_surface)
 345     _surface->Delete();
 346 }
 347
 348 vtkSimpleMode::vtkSimpleMode()
 349   : Factor(1.)
 350   , AnimationTime(0.)
 351   , Internal(new vtkSimpleMode::vtkSimpleModeInternal)
 352 {
 353   // this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,vtkDataSetAttributes::VECTORS);
 354 }
 355
 356 vtkSimpleMode::~vtkSimpleMode()
 357 {
 358   delete this->Internal;
 359 }
 360
 361 void vtkSimpleMode::SetInputArrayToProcess(
 362   int idx, int port, int connection, int ff, const char* name)
 363 {
 364   if (idx == 0)
 365     this->Internal->setFieldForReal(name);
 366   vtkDataSetAlgorithm::SetInputArrayToProcess(idx, port, connection, ff, name);
 367 }
 368
 369 double GetOptimalRatioFrom(vtkUnstructuredGrid* dataset, vtkDoubleArray* array)
 370 {
 371   if (!dataset || !array)
 372     throw MZCException("The input dataset and or array is null !");
 373   vtkDataArray* coords(dataset->GetPoints()->GetData());
 374   vtkDoubleArray* coords2(vtkDoubleArray::SafeDownCast(coords));
 375   if (!coords2)
 376     throw MZCException("Input coordinates are not float64 !");
 377   int nbCompo(array->GetNumberOfComponents());
 378   if (coords2->GetNumberOfComponents() != 3 || (nbCompo != 3 && nbCompo != 6))
 379     throw MZCException("Input coordinates do not have 3 components as it should !");
 380   int nbPts(dataset->GetNumberOfPoints());
 381   const double* srcPt1(array->Begin());
 382   dataset->ComputeBounds();
 383   double* minmax1(dataset->GetBounds());
 384   double minmax2[3] = { 0., 0., 0. };
 385   for (int i = 0; i < nbPts; i++, srcPt1 += nbCompo)
 386   {
 387     minmax2[0] = std::max(fabs(srcPt1[0]), minmax2[0]);
 388     minmax2[1] = std::max(fabs(srcPt1[1]), minmax2[1]);
 389     minmax2[2] = std::max(fabs(srcPt1[2]), minmax2[2]);
 390   }
 391   double maxDispDelta(*std::max_element(minmax2, minmax2 + 3));
 392   if (maxDispDelta < EPS)
 393     maxDispDelta = 1.;
 394   for (int i = 0; i < 3; i++)
 395     minmax2[i] = minmax1[2 * i + 1] - minmax1[2 * i];
 396   double maxGeoDelta(*std::max_element(minmax2, minmax2 + 3));
 397   if (maxDispDelta < EPS)
 398     maxDispDelta = 1.;
 399   return maxGeoDelta / maxDispDelta;
 400 }
 401
 402 int vtkSimpleMode::FillOutputPortInformation(int vtkNotUsed(port), vtkInformation* info)
 403 {
 404   info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData");
 405   return 1;
 406 }
 407
 408 int vtkSimpleMode::RequestInformation(
 409       vtkInformation* /*request*/, vtkInformationVector** /*inputVector*/, vtkInformationVector* /*outputVector*/)
 410 {
 411   // std::cerr << "########################################## vtkSimpleMode::RequestInformation
 412   // ##########################################" << std::endl;
 413   try
 414   {
 415     if (this->Internal->getFieldForReal().empty())
 416       return 1;
 417     /*vtkUnstructuredGrid *usgIn(0);
 418     vtkDoubleArray *arr(0);
 419     ExtractInfo(inputVector[0],usgIn,this->Internal->getFieldForReal(),arr);
 420     std::vector<std::string> candidatesArrName(GetPossibleArrayNames(usgIn));
 421     //
 422     double ratio(GetOptimalRatioFrom(usgIn,arr));
 423     std::string optArrNameForReal(FindBestRealAmong(candidatesArrName));
 424     std::string optArrNameForImag(FindBestImagAmong(candidatesArrName));*/
 425     // std::cerr << ratio << std::endl;
 426     // std::cerr << optArrNameForReal << " * " << optArrNameForImag << std::endl;
 427   }
 428   catch (MZCException& e)
 429   {
 430     std::ostringstream oss;
 431     oss << "Exception has been thrown in vtkSimpleMode::RequestInformation : " << e.what()
 432         << std::endl;
 433     if (this->HasObserver("ErrorEvent"))
 434       this->InvokeEvent("ErrorEvent", const_cast<char*>(oss.str().c_str()));
 435     else
 436       vtkOutputWindowDisplayErrorText(const_cast<char*>(oss.str().c_str()));
 437     vtkObject::BreakOnError();
 438     return 0;
 439   }
 440   return 1;
 441 }
 442
 443 int vtkSimpleMode::RequestData(
 444     vtkInformation* /*request*/, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
 445 {
 446   // std::cerr << "########################################## vtkSimpleMode::RequestData
 447   // ##########################################" << std::endl;
 448   try
 449   {
 450     vtkUnstructuredGrid* usgIn(0);
 451     usgIn = ExtractInfo1(inputVector[0]);
 452     //
 453     int nbPts(usgIn->GetNumberOfPoints());
 454     vtkPolyData* outSurface(this->Internal->performConnection(usgIn));
 455     vtkDataArray* arrRealBase(nullptr);
 456     ExtractInfo2(outSurface, this->Internal->getFieldForReal(), arrRealBase);
 457     vtkSmartPointer<vtkDataArray> arrReal(ForceTo3Compo(arrRealBase));
 458     vtkSmartPointer<vtkDoubleArray> arr1(vtkSmartPointer<vtkDoubleArray>::New());
 459     arr1->SetName(ZE_DISPLACEMENT_NAME1);
 460     arr1->SetNumberOfComponents(3);
 461     arr1->SetNumberOfTuples(nbPts);
 462     double* ptToFeed1(arr1->Begin());
 463     double cst1(Factor * cos(AnimationTime * 2 * M_PI));
 464     FeedData(arrReal, cst1, ptToFeed1);
 465     int idx1(outSurface->GetPointData()->AddArray(arr1));
 466     outSurface->GetPointData()->SetActiveAttribute(idx1, vtkDataSetAttributes::VECTORS);
 467     //
 468     //
 469     vtkSmartPointer<vtkWarpVector> ws(vtkSmartPointer<vtkWarpVector>::New()); // vtkNew
 470     ws->SetInputData(outSurface);
 471     ws->SetScaleFactor(1.);
 472     ws->SetInputArrayToProcess(
 473       idx1, 0, 0, "vtkDataObject::FIELD_ASSOCIATION_POINTS", ZE_DISPLACEMENT_NAME1);
 474     ws->Update();
 475     vtkSmartPointer<vtkDataSet> ds(ws->GetOutput());
 476     ds->GetPointData()->RemoveArray(idx1);
 477     //
 478     int idx2(0);
 479     {
 480       vtkDataArray* dummy(0);
 481       ExtractInfo3(ds, this->Internal->getFieldForReal(), dummy, idx2);
 482     }
 483     //
 484     vtkInformation* outInfo(outputVector->GetInformationObject(0));
 485     vtkPolyData* output(vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())));
 486     output->ShallowCopy(ds);
 487     output->GetPointData()->DeepCopy(ds->GetPointData());
 488     //
 489     for (int i = 0; i < output->GetPointData()->GetNumberOfArrays(); i++)
 490     {
 491       vtkDataArray* arr(output->GetPointData()->GetArray(i));
 492       vtkDoubleArray* arr2(vtkDoubleArray::SafeDownCast(arr));
 493       if (!arr2)
 494         continue;
 495       vtkIdType nbCompo(arr2->GetNumberOfComponents()), nbTuples(arr2->GetNumberOfTuples());
 496       if (nbCompo != 3 && nbCompo != 2)
 497         continue;
 498       double* arrPtr(arr2->GetPointer(0));
 499       std::transform(arrPtr, arrPtr + nbCompo * nbTuples, arrPtr,
 500         std::bind(std::multiplies<double>(),std::placeholders::_1, cos(AnimationTime * 2 * M_PI)));
 501     }
 502     //
 503     vtkDataArray* array = output->GetPointData()->GetArray(idx2);
 504     vtkSmartPointer<vtkDataArray> result =
 505       vtkSmartPointer<vtkDataArray>::Take(vtkDataArray::CreateDataArray(array->GetDataType()));
 506     result->ShallowCopy(array);
 507     result->SetName("__NormalModesAnimation__");
 508     output->GetPointData()->SetScalars(result);
 509   }
 510   catch (MZCException& e)
 511   {
 512     std::ostringstream oss;
 513     oss << "Exception has been thrown in vtkSimpleMode::RequestInformation : " << e.what()
 514         << std::endl;
 515     if (this->HasObserver("ErrorEvent"))
 516       this->InvokeEvent("ErrorEvent", const_cast<char*>(oss.str().c_str()));
 517     else
 518       vtkOutputWindowDisplayErrorText(const_cast<char*>(oss.str().c_str()));
 519     vtkObject::BreakOnError();
 520     return 0;
 521   }
 522   return 1;
 523 }
 524
 525 void vtkSimpleMode::PrintSelf(ostream& os, vtkIndent indent)
 526 {
 527   this->Superclass::PrintSelf(os, indent);
 528 }