src/ElectromagnetismRotation/plugin/ElectromagnetismRotationIO/vtkElectromagnetismRotation.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
  20
  21 #include "vtkElectromagnetismRotation.h"
  22
  23 #include "MEDCouplingMemArray.hxx"
  24
  25 #include "ElectromagnetismRotationHelper.h"
  26 #include "VTKMEDTraits.hxx"
  27
  28 #include "vtkAdjacentVertexIterator.h"
  29 #include "vtkDataArrayTemplate.h"
  30 #include "vtkIntArray.h"
  31 #include "vtkLongArray.h"
  32 #ifdef WIN32
  33 #include "vtkLongLongArray.h"
  34 #endif
  35 #include "vtkCellData.h"
  36 #include "vtkPointData.h"
  37
  38 #include "vtkStreamingDemandDrivenPipeline.h"
  39 #include "vtkUnstructuredGrid.h"
  40 #include  "vtkMultiBlockDataSet.h"
  41
  42 #include "vtkInformationStringKey.h"
  43 #include "vtkAlgorithmOutput.h"
  44 #include "vtkObjectFactory.h"
  45 #include "vtkMutableDirectedGraph.h"
  46 #include "vtkMultiBlockDataSet.h"
  47 #include "vtkDataSet.h"
  48 #include "vtkInformationVector.h"
  49 #include "vtkInformation.h"
  50 #include "vtkDataArraySelection.h"
  51 #include "vtkTimeStamp.h"
  52 #include "vtkInEdgeIterator.h"
  53 #include "vtkInformationDataObjectKey.h"
  54 #include "vtkExecutive.h"
  55 #include "vtkVariantArray.h"
  56 #include "vtkStringArray.h"
  57 #include "vtkDoubleArray.h"
  58 #include "vtkCharArray.h"
  59 #include "vtkUnsignedCharArray.h"
  60 #include "vtkDataSetAttributes.h"
  61 #include "vtkDemandDrivenPipeline.h"
  62 #include "vtkDataObjectTreeIterator.h"
  63 #include "vtkThreshold.h"
  64 #include "vtkMultiBlockDataGroupFilter.h"
  65 #include "vtkCompositeDataToUnstructuredGridFilter.h"
  66 #include "vtkInformationDataObjectMetaDataKey.h"
  67 #include "vtkTransform.h"
  68 #include "vtkFunctionParser.h"
  69
  70 #include <map>
  71 #include <deque>
  72 #include <cstring>
  73 #include <memory>
  74
  75 const char ZE_SEP[]="@@][@@";
  76
  77 const char TS_STR[]="TS";
  78
  79 const char COM_SUP_STR[]="ComSup";
  80
  81 const char FAMILY_ID_CELL_NAME[]="FamilyIdCell";
  82
  83 const char NUM_ID_CELL_NAME[]="NumIdCell";
  84
  85 const char FAMILY_ID_NODE_NAME[]="FamilyIdNode";
  86
  87 const char NUM_ID_NODE_NAME[]="NumIdNode";
  88
  89 const char GLOBAL_NODE_ID_NAME[]="GlobalNodeIds";
  90
  91 vtkStandardNewMacro(vtkElectromagnetismRotation);
  92
  93 vtkInformationDataObjectMetaDataKey* GetMEDReaderMetaDataIfAny()
  94 {
  95   static const char ZE_KEY[] = "vtkMEDReader::META_DATA";
  96   MEDCoupling::GlobalDict* gd(MEDCoupling::GlobalDict::GetInstance());
  97   if (!gd->hasKey(ZE_KEY))
  98     return 0;
  99   std::string ptSt(gd->value(ZE_KEY));
 100   void* pt(0);
 101   std::istringstream iss(ptSt);
 102   iss >> pt;
 103   return reinterpret_cast<vtkInformationDataObjectMetaDataKey*>(pt);
 104 }
 105
 106 bool IsInformationOK(vtkInformation* info)
 107 {
 108   vtkInformationDataObjectMetaDataKey* key(GetMEDReaderMetaDataIfAny());
 109   if (!key)
 110     return false;
 111   // Check the information contain meta data key
 112   if (!info->Has(key))
 113     return false;
 114   // Recover Meta Data
 115   vtkMutableDirectedGraph* sil(vtkMutableDirectedGraph::SafeDownCast(info->Get(key)));
 116   if (!sil)
 117     return false;
 118   int idNames(0);
 119   vtkAbstractArray* verticesNames(sil->GetVertexData()->GetAbstractArray("Names", idNames));
 120   vtkStringArray* verticesNames2(vtkStringArray::SafeDownCast(verticesNames));
 121   if (!verticesNames2)
 122     return false;
 123   for (int i = 0; i < verticesNames2->GetNumberOfValues(); i++)
 124   {
 125     vtkStdString& st(verticesNames2->GetValue(i));
 126     if (st == "MeshesFamsGrps")
 127       return true;
 128   }
 129   return false;
 130 }
 131
 132 class vtkElectromagnetismRotation::vtkElectromagnetismRotationInternal : public ElectromagnetismRotationInternal
 133 {
 134 };
 135
 136 ////////////////////
 137
 138 vtkElectromagnetismRotation::vtkElectromagnetismRotation():SIL(NULL),Internal(new vtkElectromagnetismRotationInternal),InsideOut(0),Axis(2),RotationRotor(1e300)
 139 {
 140 }
 141
 142 vtkElectromagnetismRotation::~vtkElectromagnetismRotation()
 143 {
 144   delete this->Internal;
 145 }
 146
 147 void vtkElectromagnetismRotation::SetInsideOut(int val)
 148 {
 149   if(this->InsideOut!=val)
 150     {
 151       this->InsideOut=val;
 152       this->Modified();
 153     }
 154 }
 155
 156 void vtkElectromagnetismRotation::SetAxis(int axis)
 157 {
 158   if(this->Axis!=axis)
 159     {
 160       this->Axis=axis;
 161       this->Modified();
 162     }
 163 }
 164
 165 void vtkElectromagnetismRotation::SetAngularStep(char *angStep)
 166 {
 167   this->AngularStep = angStep;
 168   this->Modified();
 169 }
 170
 171 int vtkElectromagnetismRotation::RequestInformation(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
 172 {
 173 //  vtkUnstructuredGridAlgorithm::RequestInformation(request,inputVector,outputVector);
 174   try
 175     {
 176 //      std::cerr << "########################################## vtkElectromagnetismRotation::RequestInformation ##########################################" << std::endl;
 177 //      request->Print(cout);
 178       vtkInformation *outInfo(outputVector->GetInformationObject(0));
 179       vtkInformation *inputInfo(inputVector[0]->GetInformationObject(0));
 180       if(!ElectromagnetismRotationInternal::IndependantIsInformationOK(GetMEDReaderMetaDataIfAny(),inputInfo))
 181         {
 182         vtkErrorMacro("No SIL Data available ! The source of this filter must be MEDReader !");
 183         return 0;
 184         }
 185
 186       this->SetSIL(vtkMutableDirectedGraph::SafeDownCast(inputInfo->Get(GetMEDReaderMetaDataIfAny())));
 187       this->Internal->loadFrom(this->SIL);
 188       //this->Internal->printMySelf(std::cerr);
 189     }
 190   catch(INTERP_KERNEL::Exception& e)
 191     {
 192       std::cerr << "Exception has been thrown in vtkElectromagnetismRotation::RequestInformation : " << e.what() << std::endl;
 193       return 0;
 194     }
 195   return 1;
 196 }
 197
 198 /*!
 199  * Do not use vtkCxxSetObjectMacro macro because input mdg comes from an already managed in the pipeline just a ref on it.
 200  */
 201 void vtkElectromagnetismRotation::SetSIL(vtkMutableDirectedGraph *mdg)
 202 {
 203   if(this->SIL==mdg)
 204     return ;
 205   this->SIL=mdg;
 206 }
 207
 208 template<class CellPointExtractor>
 209 vtkDataSet *FilterFamilies(vtkThreshold *thres,
 210                            vtkDataSet *input, const std::set<int>& idsToKeep, bool insideOut, const char *arrNameOfFamilyField,
 211                            const char *associationForThreshold, bool& catchAll, bool& catchSmth)
 212 {
 213   const int VTK_DATA_ARRAY_DELETE=vtkDataArrayTemplate<double>::VTK_DATA_ARRAY_DELETE;
 214   const char ZE_SELECTION_ARR_NAME[]="@@ZeSelection@@";
 215   vtkDataSet *output(input->NewInstance());
 216   output->ShallowCopy(input);
 217   thres->SetInputData(output);
 218   vtkDataSetAttributes *dscIn(input->GetCellData()),*dscIn2(input->GetPointData());
 219   vtkDataSetAttributes *dscOut(output->GetCellData()),*dscOut2(output->GetPointData());
 220   //
 221   double vMin(insideOut==0?1.:0.),vMax(insideOut==0?2.:1.);
 222   thres->ThresholdBetween(vMin,vMax);
 223   // OK for the output
 224   //
 225   CellPointExtractor cpe2(input);
 226   vtkDataArray *da(cpe2.Get()->GetScalars(arrNameOfFamilyField));
 227   if(!da)
 228     return 0;
 229   std::string daName(da->GetName());
 230   typedef MEDFileVTKTraits<mcIdType>::VtkType vtkMCIdTypeArray;
 231   vtkMCIdTypeArray *dai(vtkMCIdTypeArray::SafeDownCast(da));
 232   if(daName!=arrNameOfFamilyField || !dai)
 233     return 0;
 234   //
 235   int nbOfTuples(dai->GetNumberOfTuples());
 236   vtkCharArray *zeSelection(vtkCharArray::New());
 237   zeSelection->SetName(ZE_SELECTION_ARR_NAME);
 238   zeSelection->SetNumberOfComponents(1);
 239   char *pt(new char[nbOfTuples]);
 240   zeSelection->SetArray(pt,nbOfTuples,0,VTK_DATA_ARRAY_DELETE);
 241   const mcIdType *inPtr(dai->GetPointer(0));
 242   std::fill(pt,pt+nbOfTuples,0);
 243   catchAll=true; catchSmth=false;
 244   std::vector<bool> pt2(nbOfTuples,false);
 245   for(std::set<int>::const_iterator it=idsToKeep.begin();it!=idsToKeep.end();it++)
 246     {
 247       bool catchFid(false);
 248       for(int i=0;i<nbOfTuples;i++)
 249         if(inPtr[i]==*it)
 250           { pt2[i]=true; catchFid=true; }
 251       if(!catchFid)
 252         catchAll=false;
 253       else
 254         catchSmth=true;
 255     }
 256   for(int ii=0;ii<nbOfTuples;ii++)
 257     if(pt2[ii])
 258       pt[ii]=2;
 259   CellPointExtractor cpe3(output);
 260   int idx(cpe3.Get()->AddArray(zeSelection));
 261   cpe3.Get()->SetActiveAttribute(idx,vtkDataSetAttributes::SCALARS);
 262   cpe3.Get()->CopyScalarsOff();
 263   zeSelection->Delete();
 264   //
 265   thres->SetInputArrayToProcess(idx,0,0,associationForThreshold,ZE_SELECTION_ARR_NAME);
 266   thres->Update();
 267   vtkUnstructuredGrid *zeComputedOutput(thres->GetOutput());
 268   CellPointExtractor cpe(zeComputedOutput);
 269   cpe.Get()->RemoveArray(idx);
 270   output->Delete();
 271   zeComputedOutput->Register(0);
 272   return zeComputedOutput;
 273 }
 274
 275 class CellExtractor
 276 {
 277 public:
 278   CellExtractor(vtkDataSet *ds):_ds(ds) { }
 279   vtkDataSetAttributes *Get() { return _ds->GetCellData(); }
 280 private:
 281   vtkDataSet *_ds;
 282 };
 283
 284 class PointExtractor
 285 {
 286 public:
 287   PointExtractor(vtkDataSet *ds):_ds(ds) { }
 288   vtkDataSetAttributes *Get() { return _ds->GetPointData(); }
 289 private:
 290   vtkDataSet *_ds;
 291 };
 292
 293 int vtkElectromagnetismRotation::RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
 294 {
 295   try
 296     {
 297       // std::cerr << "########################################## vtkElectromagnetismRotation::RequestData        ##########################################" << std::endl;
 298       // request->Print(cout);
 299       vtkInformation* inputInfo=inputVector[0]->GetInformationObject(0);
 300       vtkMultiBlockDataSet *inputMB(vtkMultiBlockDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
 301       if(inputMB->GetNumberOfBlocks()!=1)
 302         {
 303           vtkErrorMacro(<< "vtkElectromagnetismRotation::RequestData : input has not the right number of parts ! Expected 1 !" ) ;
 304           return 0;
 305         }
 306       vtkDataSet *input(vtkDataSet::SafeDownCast(inputMB->GetBlock(0)));
 307       vtkInformation *info(input->GetInformation());
 308       vtkInformation *outInfo(outputVector->GetInformationObject(0));
 309       vtkMultiBlockDataSet *output(vtkMultiBlockDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())));
 310       output->SetNumberOfBlocks(2);
 311       std::set<int> idsToKeep(this->Internal->getIdsToKeep());
 312       this->Internal->clearSelection();
 313       // first shrink the input
 314       bool catchAll,catchSmth;
 315       vtkNew<vtkThreshold> thres1,thres2;
 316       vtkSmartPointer<vtkDataSet> rotor(FilterFamilies<CellExtractor>(thres1,input,idsToKeep,0,FAMILY_ID_CELL_NAME,"vtkDataObject::FIELD_ASSOCIATION_CELLS",catchAll,catchSmth));
 317       vtkSmartPointer<vtkDataSet> stator(FilterFamilies<CellExtractor>(thres2,input,idsToKeep,1,FAMILY_ID_CELL_NAME,"vtkDataObject::FIELD_ASSOCIATION_CELLS",catchAll,catchSmth));
 318       //
 319       double reqTS(0.);
 320       int nbOfSteps(-1);
 321       std::unique_ptr<double[]> timeSteps;
 322
 323       if(outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
 324         reqTS=outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
 325       if(outInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()))
 326       {
 327         nbOfSteps = outInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
 328         timeSteps.reset(new double[ nbOfSteps ]);
 329         outInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),timeSteps.get());
 330         //std::cerr << "nb : " << nbOfSteps << std::endl;
 331         //std::for_each(timeSteps.get(),timeSteps.get()+nbOfSteps,[](double v) { std::cerr << v << std::endl; });
 332       }
 333       if(nbOfSteps<2 || !timeSteps.get())
 334       {
 335         vtkErrorMacro(<< "vtkElectromagnetismRotation::RequestData : A temporal dataset is expected ! Here < 2 time steps found !" ) ;
 336         return 0;
 337       }
 338       // Calcul de l angle effectif de rotation
 339       double minTime(timeSteps[0]),maxTime(timeSteps[nbOfSteps-1]);
 340       if(minTime == maxTime)
 341       {
 342         vtkErrorMacro(<< "vtkElectromagnetismRotation::RequestData : minTime == maxTime !" ) ;
 343         return 0;
 344       }
 345       double angularStepD = 0;
 346       {
 347         vtkNew<vtkFunctionParser> fp;
 348         fp->SetFunction(this->AngularStep.c_str());
 349         angularStepD = fp->GetScalarResult();
 350       }
 351       double effectiveTime(reqTS); effectiveTime = std::max(effectiveTime,minTime); effectiveTime = std::min(effectiveTime,maxTime);
 352       this->RotationRotor = (angularStepD*((effectiveTime-minTime)/(maxTime-minTime)));
 353       //std::cout << "*** " << effectiveTime << " " << minTime << " " << maxTime << " " << angleDegree << std::endl << std::flush;
 354       //std::cout << "*** " << this->RotationRotor << std::endl << std::flush;
 355       //
 356       if(rotor)
 357         {
 358           if(catchAll)
 359             {
 360               vtkNew<vtkTransform> transformR;
 361               switch(this->Axis)
 362               {
 363                 case 0:
 364                 {
 365                   transformR->RotateX(this->RotationRotor);
 366                   break;
 367                 }
 368                 case 1:
 369                 {
 370                   transformR->RotateY(this->RotationRotor);
 371                   break;
 372                 }
 373                 case 2:
 374                 {
 375                   transformR->RotateZ(this->RotationRotor);
 376                   break;
 377                 }
 378                 default:
 379                 {
 380                   vtkErrorMacro(<< "vtkElectromagnetismRotation::RequestData : not recognized axis !" ) ;
 381                   return 0;
 382                 }
 383               }
 384               vtkNew<vtkPoints> newCoords;
 385               transformR->TransformPoints(vtkPointSet::SafeDownCast(rotor)->GetPoints(),newCoords);
 386               vtkPointSet::SafeDownCast(rotor)->SetPoints(newCoords);
 387               output->SetBlock(0,rotor);
 388               output->SetBlock(1,stator);
 389               return 1;
 390             }
 391           else
 392             return 0;
 393         }
 394     }
 395   catch(INTERP_KERNEL::Exception& e)
 396     {
 397       vtkErrorMacro(<< "Exception has been thrown in vtkElectromagnetismRotation::RequestData : " << e.what());
 398       return 0;
 399     }
 400 }
 401
 402 int vtkElectromagnetismRotation::GetSILUpdateStamp()
 403 {
 404   return (int)this->SILTime;
 405 }
 406
 407 const char* vtkElectromagnetismRotation::GetGrpStart()
 408 {
 409   return ElectromagnetismRotationGrp::START;
 410 }
 411
 412 const char* vtkElectromagnetismRotation::GetFamStart()
 413 {
 414   return ElectromagnetismRotationFam::START;
 415 }
 416
 417 void vtkElectromagnetismRotation::PrintSelf(ostream& os, vtkIndent indent)
 418 {
 419   this->Superclass::PrintSelf(os, indent);
 420 }
 421
 422 int vtkElectromagnetismRotation::GetNumberOfGroupsFlagsArrays()
 423 {
 424   int ret(this->Internal->getNumberOfEntries());
 425   //std::cerr << "vtkElectromagnetismRotation::GetNumberOfFieldsTreeArrays() -> " << ret << std::endl;
 426   return ret;
 427 }
 428
 429 const char *vtkElectromagnetismRotation::GetGroupsFlagsArrayName(int index)
 430 {
 431   const char *ret(this->Internal->getKeyOfEntry(index));
 432 //  std::cerr << "vtkElectromagnetismRotation::GetFieldsTreeArrayName(" << index << ") -> " << ret << std::endl;
 433   return ret;
 434 }
 435
 436 int vtkElectromagnetismRotation::GetGroupsFlagsArrayStatus(const char *name)
 437 {
 438   int ret((int)this->Internal->getStatusOfEntryStr(name));
 439 //  std::cerr << "vtkElectromagnetismRotation::GetGroupsFlagsArrayStatus(" << name << ") -> " << ret << std::endl;
 440   return ret;
 441 }
 442
 443 void vtkElectromagnetismRotation::SetGroupsFlagsStatus(const char *name, int status)
 444 {
 445   //std::cerr << "vtkElectromagnetismRotation::SetFieldsStatus(" << name << "," << status << ")" << std::endl;
 446   this->Internal->setStatusOfEntryStr(name,(bool)status);
 447   this->Modified();
 448   //this->Internal->printMySelf(std::cerr);
 449 }
 450
 451 const char *vtkElectromagnetismRotation::GetMeshName()
 452 {
 453   return this->Internal->getMeshName();
 454 }