src/Plugins/SimpleMode/IO/vtkSimpleMode.cxx

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