]> SALOME platform Git repositories - modules/paravis.git/blob - src/Plugins/GaussToCell/vtkGaussToCell.cxx
Salome HOME
Add missing MEDFile dependency in DevelopedSurface filter CMakeLists.txt
[modules/paravis.git] / src / Plugins / GaussToCell / vtkGaussToCell.cxx
1 // Copyright (C) 2018  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 "vtkGaussToCell.h"
22
23 #include "vtkAdjacentVertexIterator.h"
24 #include "vtkIntArray.h"
25 #include "vtkCellData.h"
26 #include "vtkPointData.h"
27 #include "vtkCellType.h"
28 #include "vtkCell.h"
29 #include "vtkCellArray.h"
30 #include "vtkIdTypeArray.h"
31
32 #include "vtkStreamingDemandDrivenPipeline.h"
33 #include "vtkUnstructuredGrid.h"
34 #include  "vtkMultiBlockDataSet.h"
35
36 #include "vtkInformationStringKey.h"
37 #include "vtkAlgorithmOutput.h"
38 #include "vtkObjectFactory.h"
39 #include "vtkMutableDirectedGraph.h"
40 #include "vtkMultiBlockDataSet.h"
41 #include "vtkDataSet.h"
42 #include "vtkInformationVector.h"
43 #include "vtkInformation.h"
44 #include "vtkDataArraySelection.h"
45 #include "vtkTimeStamp.h"
46 #include "vtkInEdgeIterator.h"
47 #include "vtkInformationDataObjectKey.h"
48 #include "vtkInformationDataObjectMetaDataKey.h"
49 #include "vtkInformationDoubleVectorKey.h"
50 #include "vtkExecutive.h"
51 #include "vtkVariantArray.h"
52 #include "vtkStringArray.h"
53 #include "vtkDoubleArray.h"
54 #include "vtkFloatArray.h"
55 #include "vtkCharArray.h"
56 #include "vtkLongArray.h"
57 #include "vtkUnsignedCharArray.h"
58 #include "vtkDataSetAttributes.h"
59 #include "vtkDemandDrivenPipeline.h"
60 #include "vtkDataObjectTreeIterator.h"
61 #include "vtkWarpScalar.h"
62 #include "vtkQuadratureSchemeDefinition.h"
63 #include "vtkInformationQuadratureSchemeDefinitionVectorKey.h"
64 #include "vtkCompositeDataToUnstructuredGridFilter.h"
65 #include "vtkMultiBlockDataGroupFilter.h"
66
67 #include "MEDCouplingMemArray.hxx"
68 #include "MEDCouplingUMesh.hxx"
69 #include "MEDCouplingFieldDouble.hxx"
70 #include "InterpKernelAutoPtr.hxx"
71 #include "InterpKernelGaussCoords.hxx"
72
73 #include <map>
74 #include <set>
75 #include <deque>
76 #include <sstream>
77
78 using MEDCoupling::DataArray;
79 using MEDCoupling::DataArrayInt;
80 using MEDCoupling::DataArrayDouble;
81 using MEDCoupling::MEDCouplingMesh;
82 using MEDCoupling::MEDCouplingUMesh;
83 using MEDCoupling::DynamicCastSafe;
84 using MEDCoupling::MEDCouplingFieldDouble;
85 using MEDCoupling::ON_GAUSS_PT;
86 using MEDCoupling::MCAuto;
87
88 vtkStandardNewMacro(vtkGaussToCell);
89
90 vtkInformationDoubleVectorKey *GetMEDReaderMetaDataIfAny()
91 {
92   static const char ZE_KEY[]="vtkMEDReader::GAUSS_DATA";
93   MEDCoupling::GlobalDict *gd(MEDCoupling::GlobalDict::GetInstance());
94   if(!gd->hasKey(ZE_KEY))
95     return 0;
96   std::string ptSt(gd->value(ZE_KEY));
97   void *pt(0);
98   std::istringstream iss(ptSt); iss >> pt;
99   return reinterpret_cast<vtkInformationDoubleVectorKey *>(pt);
100 }
101
102 bool IsInformationOK(vtkInformation *info, std::vector<double>& data)
103 {
104   vtkInformationDoubleVectorKey *key(GetMEDReaderMetaDataIfAny());
105   if(!key)
106     return false;
107   // Check the information contain meta data key
108   if(!info->Has(key))
109     return false;
110   int lgth(key->Length(info));
111   const double *data2(info->Get(key));
112   data.insert(data.end(),data2,data2+lgth);
113   return true;
114 }
115
116 void ExtractInfo(vtkInformationVector *inputVector, vtkUnstructuredGrid *& usgIn)
117 {
118   vtkInformation *inputInfo(inputVector->GetInformationObject(0));
119   vtkDataSet *input(0);
120   vtkDataSet *input0(vtkDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
121   vtkMultiBlockDataSet *input1(vtkMultiBlockDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
122   if(input0)
123     input=input0;
124   else
125     {
126       if(!input1)
127         throw INTERP_KERNEL::Exception("Input dataSet must be a DataSet or single elt multi block dataset expected !");
128       if(input1->GetNumberOfBlocks()!=1)
129         throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with not exactly one block ! Use MergeBlocks or ExtractBlocks filter before calling this filter !");
130       vtkDataObject *input2(input1->GetBlock(0));
131       if(!input2)
132         throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with exactly one block but this single element is NULL !");
133       vtkDataSet *input2c(vtkDataSet::SafeDownCast(input2));
134       if(!input2c)
135         throw INTERP_KERNEL::Exception("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 !");
136       input=input2c;
137     }
138   if(!input)
139     throw INTERP_KERNEL::Exception("Input data set is NULL !");
140   usgIn=vtkUnstructuredGrid::SafeDownCast(input);
141   if(!usgIn)
142     throw INTERP_KERNEL::Exception("Input data set is not an unstructured mesh ! This filter works only on unstructured meshes !");
143 }
144
145 vtkGaussToCell::vtkGaussToCell():avgStatus(true),maxStatus(false),minStatus(false)
146 {
147   this->SetNumberOfInputPorts(1);
148   this->SetNumberOfOutputPorts(1);
149 }
150
151 vtkGaussToCell::~vtkGaussToCell()
152 {
153 }
154
155 void vtkGaussToCell::SetAvgFlag(bool avgStatus)
156 {
157   if(this->avgStatus!=avgStatus)
158     {
159       this->avgStatus=avgStatus;
160       this->Modified();
161     }
162 }
163
164 void vtkGaussToCell::SetMaxFlag(bool maxStatus)
165 {
166   if(this->maxStatus!=maxStatus)
167     {
168       this->maxStatus=maxStatus;
169       this->Modified();
170     }
171 }
172
173 void vtkGaussToCell::SetMinFlag(bool minStatus)
174 {
175   if(this->minStatus!=minStatus)
176     {
177       this->minStatus=minStatus;
178       this->Modified();
179     }
180 }
181
182 int vtkGaussToCell::RequestInformation(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
183
184   //std::cerr << "########################################## vtkGaussToCell::RequestInformation ##########################################" << std::endl;
185   try
186     {
187       vtkUnstructuredGrid *usgIn(0);
188       ExtractInfo(inputVector[0],usgIn);
189     }
190   catch(INTERP_KERNEL::Exception& e)
191     {
192       std::ostringstream oss;
193       oss << "Exception has been thrown in vtkGaussToCell::RequestInformation : " << e.what() << std::endl;
194       if(this->HasObserver("ErrorEvent") )
195         this->InvokeEvent("ErrorEvent",const_cast<char *>(oss.str().c_str()));
196       else
197         vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
198       vtkObject::BreakOnError();
199       return 0;
200     }
201   return 1;
202 }
203
204 typedef void (*DataComputer) (const double *inData, const vtkIdType *offData, const std::vector<int> *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData);
205
206 void ComputeAvg(const double *inData, const vtkIdType *offData, const std::vector<int> *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData)
207 {
208   std::fill(outData,outData+outNbCells*zeNbCompo,0.);
209   for(auto i=0;i<outNbCells;i++,outData+=zeNbCompo)
210     {
211       auto NbGaussPt((*nbgPerCell)[i]);
212       for(auto j=0;j<NbGaussPt;j++)
213         std::transform(inData+(offData[i]+j)*zeNbCompo,inData+(offData[i]+(j+1))*zeNbCompo,outData,outData,[](const double& v, const double& w) { return v+w; });
214       double deno(1./(double)NbGaussPt);
215       std::transform(outData,outData+zeNbCompo,outData,[deno](const double& v) { return v*deno; });
216     }
217 }
218
219 void ComputeMax(const double *inData, const vtkIdType *offData, const std::vector<int> *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData)
220 {
221   std::fill(outData,outData+outNbCells*zeNbCompo,-std::numeric_limits<double>::max());
222   for(auto i=0;i<outNbCells;i++,outData+=zeNbCompo)
223     {
224       auto NbGaussPt((*nbgPerCell)[i]);
225       for(auto j=0;j<NbGaussPt;j++)
226         std::transform(inData+(offData[i]+j)*zeNbCompo,inData+(offData[i]+(j+1))*zeNbCompo,outData,outData,[](const double& v, const double& w) { return std::max(v,w); });
227     }
228 }
229
230 void ComputeMin(const double *inData, const vtkIdType *offData, const std::vector<int> *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData)
231 {
232   std::fill(outData,outData+outNbCells*zeNbCompo,std::numeric_limits<double>::max());
233   for(auto i=0;i<outNbCells;i++,outData+=zeNbCompo)
234     {
235       auto NbGaussPt((*nbgPerCell)[i]);
236       for(auto j=0;j<NbGaussPt;j++)
237         std::transform(inData+(offData[i]+j)*zeNbCompo,inData+(offData[i]+(j+1))*zeNbCompo,outData,outData,[](const double& v, const double& w) { return std::min(v,w); });
238     }
239 }
240
241 void DealWith(const char *postName, vtkDoubleArray *zearray, vtkIdTypeArray *offsets, const std::vector<int> *nbgPerCell, vtkIdType outNbCells, vtkDoubleArray *zeOutArray, DataComputer dc)
242 {
243   std::ostringstream oss;
244   oss << zearray->GetName() << '_' << postName;
245   int zeNbCompo(zearray->GetNumberOfComponents());
246   {
247     std::string st(oss.str());
248     zeOutArray->SetName(st.c_str());
249   }
250   zeOutArray->SetNumberOfComponents(zeNbCompo);
251   zeOutArray->SetNumberOfTuples(outNbCells);
252   double *outData(zeOutArray->GetPointer(0));
253   const double *inData(zearray->GetPointer(0));
254   const auto *offData(offsets->GetPointer(0));
255   for(auto i=0;i<zeNbCompo;i++)
256     {
257       const char *comp(zearray->GetComponentName(i));
258       if(comp)
259         zeOutArray->SetComponentName(i,comp);
260     }
261   dc(inData,offData,nbgPerCell,zeNbCompo,outNbCells,outData);
262 }
263
264 int vtkGaussToCell::RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
265 {
266   //std::cerr << "########################################## vtkGaussToCell::RequestData        ##########################################" << std::endl;
267   try
268     {
269       std::vector<double> GaussAdvData;
270       bool isOK(IsInformationOK(inputVector[0]->GetInformationObject(0),GaussAdvData));
271       if(!isOK)
272         throw INTERP_KERNEL::Exception("Sorry but no advanced gauss info found ! Expect to be called right after a MEDReader containing Gauss Points !");
273       vtkUnstructuredGrid *usgIn(0);
274       ExtractInfo(inputVector[0],usgIn);
275       vtkInformation *outInfo(outputVector->GetInformationObject(0));
276       vtkUnstructuredGrid *output(vtkUnstructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())));
277       output->ShallowCopy(usgIn);
278       //
279       std::string zeArrOffset;
280       int nArrays(usgIn->GetFieldData()->GetNumberOfArrays());
281       std::map<vtkIdTypeArray *,std::vector<int> > offsetKeyMap;//Map storing for each offsets array the corresponding nb of Gauss Points per cell
282       for(int i=0;i<nArrays;i++)
283         {
284           vtkDataArray *array(usgIn->GetFieldData()->GetArray(i));
285           if(!array)
286             continue;
287           const char* arrayOffsetName(array->GetInformation()->Get(vtkQuadratureSchemeDefinition::QUADRATURE_OFFSET_ARRAY_NAME()));
288           if(!arrayOffsetName)
289             continue;
290           std::string arrOffsetNameCpp(arrayOffsetName);
291           if(arrOffsetNameCpp.find("ELGA@")==std::string::npos)
292             continue;
293           if(zeArrOffset.empty())
294             zeArrOffset=arrOffsetNameCpp;
295           else
296             if(zeArrOffset!=arrOffsetNameCpp)
297               {
298                 throw INTERP_KERNEL::Exception("ComputeGaussToCell : error in QUADRATURE_OFFSET_ARRAY_NAME for Gauss fields array !");
299               }
300           vtkDataArray *offTmp(usgIn->GetCellData()->GetArray(zeArrOffset.c_str()));
301           if(!offTmp)
302             {
303               std::ostringstream oss; oss << "ComputeGaussToCell : cell field " << zeArrOffset << " not found !";
304               throw INTERP_KERNEL::Exception(oss.str());
305             }
306           vtkIdTypeArray *offsets(vtkIdTypeArray::SafeDownCast(offTmp));
307           if(!offsets)
308             {
309               std::ostringstream oss; oss << "ComputeGaussToCell : cell field " << zeArrOffset << " exists but not with the right type of data !";
310               throw INTERP_KERNEL::Exception(oss.str());
311             }
312           vtkDoubleArray *zearray(vtkDoubleArray::SafeDownCast(array));
313           if(!zearray)
314             continue ;
315           //
316           std::map<vtkIdTypeArray *,std::vector<int> >::iterator nbgPerCellPt(offsetKeyMap.find(offsets));
317           const std::vector<int> *nbgPerCell(nullptr);
318           if(nbgPerCellPt==offsetKeyMap.end())
319             {
320               // fini la parlote
321               vtkInformation *info(offsets->GetInformation());
322               if(!info)
323                 throw INTERP_KERNEL::Exception("info is null ! Internal error ! Looks bad !");
324               vtkInformationQuadratureSchemeDefinitionVectorKey *key(vtkQuadratureSchemeDefinition::DICTIONARY());
325               if(!key->Has(info))
326                 throw INTERP_KERNEL::Exception("No quadrature key in info included in offets array ! Internal error ! Looks bad !");
327               int dictSize(key->Size(info));
328               INTERP_KERNEL::AutoPtr<vtkQuadratureSchemeDefinition *> dict(new vtkQuadratureSchemeDefinition *[dictSize]);
329               key->GetRange(info,dict,0,0,dictSize);
330               auto nbOfCells(output->GetNumberOfCells());
331               std::vector<int> nbg(nbOfCells);
332               for(auto cellId=0;cellId<nbOfCells;cellId++)
333                 {
334                   int ct(output->GetCellType(cellId));
335                   vtkQuadratureSchemeDefinition *gaussLoc(dict[ct]);
336                   if(!gaussLoc)
337                     {
338                       std::ostringstream oss; oss << "For cell " << cellId << " no Gauss info attached !";
339                       throw INTERP_KERNEL::Exception(oss.str());
340                     }
341                   int np(gaussLoc->GetNumberOfQuadraturePoints());
342                   nbg[cellId]=np;
343                 }
344               nbgPerCell=&((*(offsetKeyMap.emplace(offsets,std::move(nbg)).first)).second);
345             }
346           else
347             {
348               nbgPerCell=&((*nbgPerCellPt).second);
349             }
350           auto outNbCells(nbgPerCell->size());
351           if(this->avgStatus)
352             {
353               vtkNew<vtkDoubleArray> zeOutArray;
354               DealWith("avg",zearray,offsets,nbgPerCell,outNbCells,zeOutArray,ComputeAvg);
355               output->GetCellData()->AddArray(zeOutArray);
356             }
357           if(this->maxStatus)
358             {
359               vtkNew<vtkDoubleArray> zeOutArray;
360               DealWith("max",zearray,offsets,nbgPerCell,outNbCells,zeOutArray,ComputeMax);
361               output->GetCellData()->AddArray(zeOutArray);
362             }
363           if(this->minStatus)
364             {
365               vtkNew<vtkDoubleArray> zeOutArray;
366               DealWith("min",zearray,offsets,nbgPerCell,outNbCells,zeOutArray,ComputeMin);
367               output->GetCellData()->AddArray(zeOutArray);
368             }
369         }
370     }
371   catch(INTERP_KERNEL::Exception& e)
372     {
373       std::ostringstream oss;
374       oss << "Exception has been thrown in vtkGaussToCell::RequestData : " << e.what() << std::endl;
375       if(this->HasObserver("ErrorEvent") )
376         this->InvokeEvent("ErrorEvent",const_cast<char *>(oss.str().c_str()));
377       else
378         vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
379       vtkObject::BreakOnError();
380       return 0;
381     }
382   return 1;
383 }
384
385 void vtkGaussToCell::PrintSelf(ostream& os, vtkIndent indent)
386 {
387   this->Superclass::PrintSelf(os, indent);
388 }