1 // Copyright (C) 2021 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (EDF R&D)
21 #include "vtkElectromagnetismFluxDisc.h"
23 #include "vtkAdjacentVertexIterator.h"
24 #include "vtkIntArray.h"
25 #include "vtkLongArray.h"
26 #include "vtkCellData.h"
27 #include "vtkPointData.h"
28 #include "vtkCylinder.h"
30 #include "vtkCutter.h"
31 #include "vtkTransform.h"
33 #include "vtkStreamingDemandDrivenPipeline.h"
34 #include "vtkUnstructuredGrid.h"
35 #include "vtkMultiBlockDataSet.h"
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"
71 #define _USE_MATH_DEFINES
80 vtkStandardNewMacro(vtkElectromagnetismFluxDisc);
82 class VTK_EXPORT MZCException : public std::exception
85 MZCException(const std::string& s):_reason(s) { }
86 virtual const char *what() const throw() { return _reason.c_str(); }
87 virtual ~MZCException() throw() { }
93 void ExtractInfo(vtkInformationVector *inputVector, vtkDataSet *& usgIn)
95 vtkInformation *inputInfo(inputVector->GetInformationObject(0));
97 vtkDataSet *input0(vtkDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
98 vtkMultiBlockDataSet *input1(vtkMultiBlockDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
104 throw MZCException("Input dataSet must be a DataSet or single elt multi block dataset expected !");
105 if(input1->GetNumberOfBlocks()!=1)
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 !");
110 vtkDataObject *input2(input1->GetBlock(0));
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));
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 !");
119 throw MZCException("Input data set is NULL !");
120 vtkPointData *att(input->GetPointData());
122 throw MZCException("Input dataset has no point data attribute ! Impossible to deduce a developed surface on it !");
126 class vtkElectromagnetismFluxDisc::vtkInternals
129 vtkNew<vtkCutter> Cutter;
134 vtkElectromagnetismFluxDisc::vtkElectromagnetismFluxDisc():_cyl(nullptr),Internal(new vtkInternals),RadialResolution(80),CircumferentialResolution(80)
138 vtkElectromagnetismFluxDisc::~vtkElectromagnetismFluxDisc()
140 delete this->Internal;
143 int vtkElectromagnetismFluxDisc::RequestInformation(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
145 //std::cerr << "########################################## vtkElectromagnetismFluxDisc::RequestInformation ##########################################" << std::endl;
148 vtkDataSet *usgIn(0);
149 ExtractInfo(inputVector[0],usgIn);
151 catch(MZCException& e)
153 vtkErrorMacro("Exception has been thrown in vtkElectromagnetismFluxDisc::RequestInformation : " << e.what());
159 double ComputeFlux(vtkIdType nbOfTuples, const double *area, const double *vector3Field, const double axis[3])
162 for( vtkIdType i = 0 ; i < nbOfTuples ; ++i )
163 ret += area[i] * vtkMath::Dot(vector3Field+i*3,axis);
168 void Rotate3DAlg(const double *center, const double *vect, double angle, vtkIdType nbNodes, const T *coordsIn, T *coordsOut)
170 double sina(sin(angle));
171 double cosa(cos(angle));
172 double vectorNorm[3];
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.
193 for(vtkIdType i=0; i<nbNodes; i++)
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];
202 int vtkElectromagnetismFluxDisc::RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
204 //std::cerr << "########################################## vtkElectromagnetismFluxDisc::RequestData ##########################################" << std::endl;
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());
213 _cyl->GetCenter(center);
214 _cyl->GetAxis(axis[0],axis[1],axis[2]);
215 radius=_cyl->GetRadius();
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; });
223 if(std::isnan(axis[0]) && std::isnan(axis[1]) && std::isnan(axis[2]))
224 { axis[0]=0.; axis[1]=-1.; axis[2]=0.; }
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);
231 vtkNew<vtkUnstructuredGrid> outputMesh;
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)
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)
244 coordsPtr[0] = c*(double(radI+1)/double(this->RadialResolution))*radius;
245 coordsPtr[1] = s*(double(radI+1)/double(this->RadialResolution))*radius;
250 vtkNew<vtkPoints> pts;
251 pts->SetData(coords);
252 outputMesh->SetPoints(pts);
254 vtkIdType nbOfCells(this->CircumferentialResolution*this->RadialResolution);
255 vtkNew<vtkCellArray> cells;
256 vtkNew<vtkIdTypeArray> cellsData;
257 vtkNew<vtkIdTypeArray> cellLocations;
258 vtkNew<vtkUnsignedCharArray> cellTypes;
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 )
268 vtkIdType zeDelta(iCirc*deltaPt);
269 for( int iRadial = 0 ; iRadial < this->RadialResolution ; ++iRadial )
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);
282 *ctPtr++ = VTK_TRIANGLE;
283 cdPtr[0] = 3 ; cdPtr[1] = 0 ; cdPtr[2] = zeDelta + 1; cdPtr[3] = (zeDelta + deltaPt +1);
289 vtkIdType zeDelta((this->CircumferentialResolution - 1)*deltaPt);
290 for( int iRadial = 0 ; iRadial < this->RadialResolution ; ++iRadial )
296 cdPtr[0] = 4 ; cdPtr[1] = zeDelta + iRadial ; cdPtr[2] = zeDelta + iRadial+1;
297 cdPtr[3] = iRadial+1; cdPtr[4] = iRadial;
303 *ctPtr++ = VTK_TRIANGLE;
304 cdPtr[0] = 3 ; cdPtr[1] = 0 ; cdPtr[2] = zeDelta + 1; cdPtr[3] = 1;
310 cells->SetCells(nbOfCells,cellsData);
311 outputMesh->SetCells(cellTypes,cellLocations,cells);
315 vtkIdType nbPoints(outputMesh->GetNumberOfPoints());
316 vtkMath::Cross(ZVec,axis,orthoAxis);
317 double normOrthoAxis( vtkMath::Norm(orthoAxis) );
318 if(normOrthoAxis > 1e-5)
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);
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]; }
342 vtkNew<vtkResampleWithDataSet> probeFilter;
343 probeFilter->SetInputData(outputMesh);
344 probeFilter->SetSourceData(usgIn);
346 vtkNew<vtkPointDataToCellData> pd2cd;
347 pd2cd->SetInputConnection(probeFilter->GetOutputPort());
350 vtkNew<vtkMeshQuality> mq;
351 mq->SetInputData(outputMesh);
352 mq->SetTriangleQualityMeasureToArea();
353 mq->SetQuadQualityMeasureToArea();
355 double *area( ( vtkDoubleArray::SafeDownCast(mq->GetOutput()->GetCellData()->GetArray("Quality")) )->GetPointer(0) );
357 vtkInformation *outInfo(outputVector->GetInformationObject(0));
358 vtkUnstructuredGrid *output(vtkUnstructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())));
359 output->ShallowCopy(outputMesh);
361 vtkIdType nbOfCells(output->GetNumberOfCells());
362 vtkFieldData* dsa(pd2cd->GetOutput()->GetCellData());
363 int nbOfArrays(dsa->GetNumberOfArrays());
364 for(int i = 0 ; i < nbOfArrays ; ++i )
366 vtkDoubleArray *arr( vtkDoubleArray::SafeDownCast(dsa->GetArray(i)) );
367 if( arr && arr->GetNumberOfComponents() == 3 )
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);
382 catch(MZCException& e)
384 vtkErrorMacro("Exception has been thrown in vtkElectromagnetismFluxDisc::RequestInformation : " << e.what());
390 void vtkElectromagnetismFluxDisc::PrintSelf(ostream& os, vtkIndent indent)
392 this->Superclass::PrintSelf(os, indent);
395 void vtkElectromagnetismFluxDisc::SetCutFunction(vtkImplicitFunction* func)
397 vtkCylinder *cyl(vtkCylinder::SafeDownCast(func));
405 vtkMTimeType vtkElectromagnetismFluxDisc::GetMTime()
407 vtkMTimeType maxMTime = this->Superclass::GetMTime(); // My MTime
410 maxMTime=std::max(maxMTime,_cyl->GetMTime());