src/ElectromagnetismFluxDisc/plugin/ElectromagnetismFluxDiscModule/vtkElectromagnetismFluxDisc.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 "vtkElectromagnetismFluxDisc.h"
  22
  23 #include "vtkAdjacentVertexIterator.h"
  24 #include "vtkIntArray.h"
  25 #include "vtkLongArray.h"
  26 #include "vtkCellData.h"
  27 #include "vtkPointData.h"
  28 #include "vtkCylinder.h"
  29 #include "vtkNew.h"
  30 #include "vtkCutter.h"
  31 #include "vtkTransform.h"
  32
  33 #include "vtkStreamingDemandDrivenPipeline.h"
  34 #include "vtkUnstructuredGrid.h"
  35 #include  "vtkMultiBlockDataSet.h"
  36
  37 #include "vtkInformationStringKey.h"
  38 #include "vtkAlgorithmOutput.h"
  39 #include "vtkObjectFactory.h"
  40 #include "vtkMutableDirectedGraph.h"
  41 #include "vtkMultiBlockDataSet.h"
  42 #include "vtkDataSet.h"
  43 #include "vtkInformationVector.h"
  44 #include "vtkInformation.h"
  45 #include "vtkDataArraySelection.h"
  46 #include "vtkTimeStamp.h"
  47 #include "vtkInEdgeIterator.h"
  48 #include "vtkInformationDataObjectKey.h"
  49 #include "vtkExecutive.h"
  50 #include "vtkVariantArray.h"
  51 #include "vtkStringArray.h"
  52 #include "vtkDoubleArray.h"
  53 #include "vtkFloatArray.h"
  54 #include "vtkCharArray.h"
  55 #include "vtkUnsignedCharArray.h"
  56 #include "vtkDataSetAttributes.h"
  57 #include "vtkDemandDrivenPipeline.h"
  58 #include "vtkDataObjectTreeIterator.h"
  59 #include "vtkWarpScalar.h"
  60 #include "vtkDiskSource.h"
  61 #include "vtkTransform.h"
  62 #include "vtkTransformPolyDataFilter.h"
  63 #include "vtkResampleWithDataSet.h"
  64 #include "vtkPointDataToCellData.h"
  65 #include "vtkMeshQuality.h"
  66
  67 #include <cmath>
  68
  69
  70 #ifdef WIN32
  71 #define _USE_MATH_DEFINES
  72 #endif
  73 #include <math.h>
  74
  75 #include <map>
  76 #include <deque>
  77 #include <sstream>
  78 #include <algorithm>
  79
  80 vtkStandardNewMacro(vtkElectromagnetismFluxDisc);
  81
  82 class VTK_EXPORT MZCException : public std::exception
  83 {
  84 public:
  85   MZCException(const std::string& s):_reason(s) { }
  86   virtual const char *what() const throw() { return _reason.c_str(); }
  87   virtual ~MZCException() throw() { }
  88 private:
  89   std::string _reason;
  90 };
  91
  92
  93 void ExtractInfo(vtkInformationVector *inputVector, vtkDataSet *& usgIn)
  94 {
  95   vtkInformation *inputInfo(inputVector->GetInformationObject(0));
  96   vtkDataSet *input(0);
  97   vtkDataSet *input0(vtkDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
  98   vtkMultiBlockDataSet *input1(vtkMultiBlockDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
  99   if(input0)
 100     input=input0;
 101   else
 102     {
 103       if(!input1)
 104         throw MZCException("Input dataSet must be a DataSet or single elt multi block dataset expected !");
 105       if(input1->GetNumberOfBlocks()!=1)
 106       {
 107         std::cerr << "**** " << input1->GetNumberOfBlocks() << std::endl;
 108         throw MZCException("Input dataSet is a multiblock dataset with not exactly one block ! Use MergeBlocks or ExtractBlocks filter before calling this filter !");
 109       }
 110       vtkDataObject *input2(input1->GetBlock(0));
 111       if(!input2)
 112         throw MZCException("Input dataSet is a multiblock dataset with exactly one block but this single element is NULL !");
 113       vtkDataSet *input2c(vtkDataSet::SafeDownCast(input2));
 114       if(!input2c)
 115         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 !");
 116       input=input2c;
 117     }
 118   if(!input)
 119     throw MZCException("Input data set is NULL !");
 120   vtkPointData *att(input->GetPointData());
 121   if(!att)
 122     throw MZCException("Input dataset has no point data attribute ! Impossible to deduce a developed surface on it !");
 123   usgIn=input;
 124 }
 125
 126 class vtkElectromagnetismFluxDisc::vtkInternals
 127 {
 128 public:
 129   vtkNew<vtkCutter> Cutter;
 130 };
 131
 132 ////////////////////
 133
 134 vtkElectromagnetismFluxDisc::vtkElectromagnetismFluxDisc():_cyl(nullptr),Internal(new vtkInternals),RadialResolution(80),CircumferentialResolution(80)
 135 {
 136 }
 137
 138 vtkElectromagnetismFluxDisc::~vtkElectromagnetismFluxDisc()
 139 {
 140   delete this->Internal;
 141 }
 142
 143 int vtkElectromagnetismFluxDisc::RequestInformation(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
 144 {
 145   //std::cerr << "########################################## vtkElectromagnetismFluxDisc::RequestInformation ##########################################" << std::endl;
 146   try
 147     {
 148       vtkDataSet *usgIn(0);
 149       ExtractInfo(inputVector[0],usgIn);
 150     }
 151   catch(MZCException& e)
 152     {
 153       vtkErrorMacro("Exception has been thrown in vtkElectromagnetismFluxDisc::RequestInformation : " << e.what());
 154       return 0;
 155     }
 156   return 1;
 157 }
 158
 159 double ComputeFlux(vtkIdType nbOfTuples, const double *area, const double *vector3Field, const double axis[3])
 160 {
 161   double ret(0.0);
 162   for( vtkIdType i = 0 ; i < nbOfTuples ; ++i )
 163     ret += area[i] * vtkMath::Dot(vector3Field+i*3,axis);
 164   return ret;
 165 }
 166
 167 template<class T>
 168 void Rotate3DAlg(const double *center, const double *vect, double angle, vtkIdType nbNodes, const T *coordsIn, T *coordsOut)
 169 {
 170   double sina(sin(angle));
 171   double cosa(cos(angle));
 172   double vectorNorm[3];
 173   T matrix[9];
 174   T matrixTmp[9];
 175   double norm(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
 176   if(norm<std::numeric_limits<double>::min())
 177     throw MZCException("Rotate3DAlg : magnitude of input vector is too close of 0. !");
 178   std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
 179   //rotation matrix computation
 180   matrix[0]=cosa; matrix[1]=0.; matrix[2]=0.; matrix[3]=0.; matrix[4]=cosa; matrix[5]=0.; matrix[6]=0.; matrix[7]=0.; matrix[8]=cosa;
 181   matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
 182   matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
 183   matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
 184   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<T>(),1-cosa));
 185   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<T>());
 186   matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
 187   matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
 188   matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
 189   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<T>(),sina));
 190   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<T>());
 191   //rotation matrix computed.
 192   T tmp[3];
 193   for(vtkIdType i=0; i<nbNodes; i++)
 194     {
 195       std::transform(coordsIn+i*3,coordsIn+(i+1)*3,center,tmp,std::minus<T>());
 196       coordsOut[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+(T)center[0];
 197       coordsOut[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+(T)center[1];
 198       coordsOut[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+(T)center[2];
 199     }
 200 }
 201
 202 int vtkElectromagnetismFluxDisc::RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
 203 {
 204   //std::cerr << "########################################## vtkElectromagnetismFluxDisc::RequestData        ##########################################" << std::endl;
 205   try
 206     {
 207       if(!_cyl)
 208         throw MZCException("No cylinder object as cut function !");
 209       double center[3],axis[3],radius,orthoAxis[3];
 210       const double ZVec[3] = {0.,0.,1.};
 211       vtkAbstractTransform* trf(_cyl->GetTransform());
 212       {
 213         _cyl->GetCenter(center);
 214         _cyl->GetAxis(axis[0],axis[1],axis[2]);
 215         radius=_cyl->GetRadius();
 216       }
 217       if(trf)
 218         {
 219           double axis3[3]={center[0]+0.,center[1]+1.,center[2]+0.},axis4[3];
 220           trf->TransformPoint(axis3,axis4);
 221           std::transform(axis4,axis4+3,center,axis,[](double a, double b) { return b-a; });
 222           axis[1]=-axis[1];
 223           if(std::isnan(axis[0]) && std::isnan(axis[1]) && std::isnan(axis[2]))
 224             { axis[0]=0.; axis[1]=-1.; axis[2]=0.; }
 225         }
 226
 227       //std::cerr << trf << " jjj " << axis[0] << " " << axis[1] << " " << axis[2]  << " : " << center[0] <<  " " << center[1] << " " << center[2] <<  " " " " << " -> " << radius << std::endl;
 228       vtkDataSet *usgIn(0);
 229       ExtractInfo(inputVector[0],usgIn);
 230       //
 231       vtkNew<vtkUnstructuredGrid> outputMesh;
 232       {
 233         vtkIdType nbPoints(this->RadialResolution*this->CircumferentialResolution+1);
 234         vtkNew<vtkDoubleArray> coords;
 235         coords->SetNumberOfComponents(3); coords->SetNumberOfTuples(nbPoints);
 236         double *coordsPtr(coords->GetPointer(0));
 237         coordsPtr[0] = 0; coordsPtr[1] = 0; coordsPtr[2] = 0; coordsPtr+=3;
 238         for(int circI = 0 ; circI < this->CircumferentialResolution ; ++circI)
 239         {
 240           double a(2*M_PI*(double(circI)/double(this->CircumferentialResolution)));
 241           double c(cos(a)),s(sin(a));
 242           for(int radI = 0 ; radI < this->RadialResolution ; ++radI)
 243           {
 244             coordsPtr[0] = c*(double(radI+1)/double(this->RadialResolution))*radius;
 245             coordsPtr[1] = s*(double(radI+1)/double(this->RadialResolution))*radius;
 246             coordsPtr[2] = 0;
 247             coordsPtr+=3;
 248           }
 249         }
 250         vtkNew<vtkPoints> pts;
 251         pts->SetData(coords);
 252         outputMesh->SetPoints(pts);
 253         //
 254         vtkIdType nbOfCells(this->CircumferentialResolution*this->RadialResolution);
 255         vtkNew<vtkCellArray> cells;
 256         vtkNew<vtkIdTypeArray> cellsData;
 257         vtkNew<vtkIdTypeArray> cellLocations;
 258         vtkNew<vtkUnsignedCharArray> cellTypes;
 259         //
 260         cellTypes->SetNumberOfComponents(1); cellTypes->SetNumberOfTuples(nbOfCells);
 261         cellLocations->SetNumberOfComponents(1); cellLocations->SetNumberOfTuples(nbOfCells);
 262         cellsData->SetNumberOfComponents(1); cellsData->SetNumberOfTuples( ( 5*(this->RadialResolution)-1 ) * this->CircumferentialResolution );
 263         vtkIdType *clPtr(cellLocations->GetPointer(0)),*cdPtr(cellsData->GetPointer(0));
 264         vtkIdType offset(0),deltaPt(this->RadialResolution);
 265         unsigned char *ctPtr(cellTypes->GetPointer(0));
 266         for( int iCirc = 0 ; iCirc < this->CircumferentialResolution - 1; ++iCirc )
 267         {
 268           vtkIdType zeDelta(iCirc*deltaPt);
 269           for( int iRadial = 0 ; iRadial < this->RadialResolution ; ++iRadial )
 270           {
 271             *clPtr++ = offset;
 272             if(iRadial!=0)
 273             {
 274               *ctPtr++ = VTK_QUAD;
 275               cdPtr[0] = 4 ; cdPtr[1] = zeDelta + iRadial ; cdPtr[2] = zeDelta + iRadial+1;
 276               cdPtr[3] = (zeDelta + deltaPt + iRadial+1); cdPtr[4] = ( zeDelta + deltaPt + iRadial);
 277               cdPtr+=5;
 278               offset += 4;
 279             }
 280             else
 281             {
 282               *ctPtr++ = VTK_TRIANGLE;
 283               cdPtr[0] = 3 ; cdPtr[1] = 0 ; cdPtr[2] = zeDelta + 1; cdPtr[3] = (zeDelta + deltaPt +1);
 284               cdPtr += 4;
 285               offset += 3;
 286             }
 287           }
 288         }
 289         vtkIdType zeDelta((this->CircumferentialResolution - 1)*deltaPt);
 290         for( int iRadial = 0 ; iRadial < this->RadialResolution ; ++iRadial )
 291           {
 292             *clPtr++ = offset;
 293             if(iRadial!=0)
 294             {
 295               *ctPtr++ = VTK_QUAD;
 296               cdPtr[0] = 4 ; cdPtr[1] = zeDelta + iRadial ; cdPtr[2] = zeDelta + iRadial+1;
 297               cdPtr[3] = iRadial+1; cdPtr[4] = iRadial;
 298               cdPtr+=5;
 299               offset += 4;
 300             }
 301             else
 302             {
 303               *ctPtr++ = VTK_TRIANGLE;
 304               cdPtr[0] = 3 ; cdPtr[1] = 0 ; cdPtr[2] = zeDelta + 1; cdPtr[3] = 1;
 305               cdPtr += 4;
 306               offset += 3;
 307             }
 308           }
 309         //
 310         cells->SetCells(nbOfCells,cellsData);
 311         outputMesh->SetCells(cellTypes,cellLocations,cells);
 312       }
 313       // Rotation
 314       {
 315         vtkIdType nbPoints(outputMesh->GetNumberOfPoints());
 316         vtkMath::Cross(ZVec,axis,orthoAxis);
 317         double normOrthoAxis( vtkMath::Norm(orthoAxis) );
 318         if(normOrthoAxis > 1e-5)
 319         {
 320           //std::cerr << "ortho : " << normOrthoAxis << " X = " << orthoAxis[0] << " Y = " << orthoAxis[1] << " Z = " << orthoAxis[2] << std::endl;
 321           orthoAxis[0] *= normOrthoAxis; orthoAxis[1] *= normOrthoAxis; orthoAxis[2] *= normOrthoAxis;
 322           const double Center[3] = {0.,0.,0.};
 323           vtkNew<vtkDoubleArray> newArray;
 324           newArray->SetNumberOfComponents(3); newArray->SetNumberOfTuples(nbPoints);
 325           vtkDoubleArray *oldPts( vtkDoubleArray::SafeDownCast( outputMesh->GetPoints()->GetData() ) );
 326           double angle(asin(normOrthoAxis));
 327           if( vtkMath::Dot(ZVec,axis) < 0. )
 328             angle = M_PI - angle;
 329           Rotate3DAlg<double>(Center,orthoAxis,angle,nbPoints,oldPts->GetPointer(0),newArray->GetPointer(0));
 330           outputMesh->GetPoints()->SetData(newArray);
 331         }
 332       }
 333       // Translation
 334       {
 335         vtkDoubleArray *coords(vtkDoubleArray::SafeDownCast( outputMesh->GetPoints()->GetData()));
 336         vtkIdType nbPts(coords->GetNumberOfTuples());
 337         double *coordsPtr(coords->GetPointer(0));
 338         for(vtkIdType i = 0 ; i < nbPts ; ++i)
 339         { coordsPtr[3*i] += center[0]; coordsPtr[3*i+1] += center[1]; coordsPtr[3*i+2] += center[2]; }
 340       }
 341       //
 342       vtkNew<vtkResampleWithDataSet> probeFilter;
 343       probeFilter->SetInputData(outputMesh);
 344       probeFilter->SetSourceData(usgIn);
 345       //
 346       vtkNew<vtkPointDataToCellData> pd2cd;
 347       pd2cd->SetInputConnection(probeFilter->GetOutputPort());
 348       pd2cd->Update();
 349       //
 350       vtkNew<vtkMeshQuality> mq;
 351       mq->SetInputData(outputMesh);
 352       mq->SetTriangleQualityMeasureToArea();
 353       mq->SetQuadQualityMeasureToArea();
 354       mq->Update();
 355       double *area( ( vtkDoubleArray::SafeDownCast(mq->GetOutput()->GetCellData()->GetArray("Quality")) )->GetPointer(0) );
 356       //
 357       vtkInformation *outInfo(outputVector->GetInformationObject(0));
 358       vtkUnstructuredGrid *output(vtkUnstructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())));
 359       output->ShallowCopy(outputMesh);
 360       //
 361       vtkIdType nbOfCells(output->GetNumberOfCells());
 362       vtkFieldData* dsa(pd2cd->GetOutput()->GetCellData());
 363       int nbOfArrays(dsa->GetNumberOfArrays());
 364       for(int i = 0 ; i < nbOfArrays ; ++i )
 365       {
 366         vtkDoubleArray *arr( vtkDoubleArray::SafeDownCast(dsa->GetArray(i)) );
 367         if( arr && arr->GetNumberOfComponents() == 3 )
 368         {
 369           vtkNew<vtkDoubleArray> arr2;
 370           arr2->ShallowCopy(arr);
 371           output->GetCellData()->AddArray(arr2);
 372           double flux(ComputeFlux(nbOfCells,area,arr->GetPointer(0),axis));
 373           std::ostringstream oss; oss << dsa->GetArrayName(i) << "_flux";
 374           vtkNew<vtkDoubleArray> arrFlux;
 375           arrFlux->SetName(oss.str().c_str());
 376           arrFlux->SetNumberOfComponents(1); arrFlux->SetNumberOfTuples(nbOfCells);
 377           std::for_each(arrFlux->GetPointer(0),arrFlux->GetPointer(nbOfCells),[flux](double& elt) { elt = flux; });
 378           output->GetCellData()->AddArray(arrFlux);
 379         }
 380       }
 381     }
 382   catch(MZCException& e)
 383     {
 384       vtkErrorMacro("Exception has been thrown in vtkElectromagnetismFluxDisc::RequestInformation : " << e.what());
 385       return 0;
 386     }
 387   return 1;
 388 }
 389
 390 void vtkElectromagnetismFluxDisc::PrintSelf(ostream& os, vtkIndent indent)
 391 {
 392   this->Superclass::PrintSelf(os, indent);
 393 }
 394
 395 void vtkElectromagnetismFluxDisc::SetCutFunction(vtkImplicitFunction* func)
 396 {
 397   vtkCylinder *cyl(vtkCylinder::SafeDownCast(func));
 398   if(cyl)
 399     {
 400       _cyl=cyl;
 401       this->Modified();
 402     }
 403 }
 404
 405 vtkMTimeType vtkElectromagnetismFluxDisc::GetMTime()
 406 {
 407   vtkMTimeType maxMTime = this->Superclass::GetMTime(); // My MTime
 408   if(_cyl)
 409     {
 410       maxMTime=std::max(maxMTime,_cyl->GetMTime());
 411     }
 412   return maxMTime;
 413 }