1 // Copyright (C) 2018-2022 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 "vtkGaussToCell.h"
23 #include "vtkAdjacentVertexIterator.h"
24 #include "vtkIntArray.h"
25 #include "vtkCellData.h"
26 #include "vtkPointData.h"
27 #include "vtkCellType.h"
29 #include "vtkCellArray.h"
30 #include "vtkIdTypeArray.h"
32 #include "vtkStreamingDemandDrivenPipeline.h"
33 #include "vtkUnstructuredGrid.h"
34 #include "vtkMultiBlockDataSet.h"
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 "vtkMergeBlocks.h"
65 #include "vtkMultiBlockDataGroupFilter.h"
67 #include "MEDCouplingMemArray.hxx"
68 #include "MEDCouplingUMesh.hxx"
69 #include "MEDCouplingFieldDouble.hxx"
70 #include "InterpKernelAutoPtr.hxx"
71 #include "InterpKernelGaussCoords.hxx"
78 using MEDCoupling::DataArray;
79 using MEDCoupling::DataArrayInt32;
80 using MEDCoupling::DataArrayInt64;
81 using MEDCoupling::DataArrayDouble;
82 using MEDCoupling::MEDCouplingMesh;
83 using MEDCoupling::MEDCouplingUMesh;
84 using MEDCoupling::DynamicCastSafe;
85 using MEDCoupling::MEDCouplingFieldDouble;
86 using MEDCoupling::ON_GAUSS_PT;
87 using MEDCoupling::MCAuto;
89 vtkStandardNewMacro(vtkGaussToCell)
91 vtkInformationDoubleVectorKey *GetMEDReaderMetaDataIfAny()
93 static const char ZE_KEY[]="vtkFileSeriesGroupReader::GAUSS_DATA";
94 MEDCoupling::GlobalDict *gd(MEDCoupling::GlobalDict::GetInstance());
95 if(!gd->hasKey(ZE_KEY))
97 std::string ptSt(gd->value(ZE_KEY));
99 std::istringstream iss(ptSt); iss >> pt;
100 return reinterpret_cast<vtkInformationDoubleVectorKey *>(pt);
103 bool IsInformationOK(vtkInformation *info, std::vector<double>& data)
105 vtkInformationDoubleVectorKey *key(GetMEDReaderMetaDataIfAny());
108 // Check the information contain meta data key
111 int lgth(key->Length(info));
112 const double *data2(info->Get(key));
113 data.insert(data.end(),data2,data2+lgth);
117 void ExtractInfo(vtkInformationVector *inputVector, vtkUnstructuredGrid *& usgIn)
119 vtkInformation *inputInfo(inputVector->GetInformationObject(0));
120 vtkDataSet *input(0);
121 vtkDataSet *input0(vtkDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
122 vtkMultiBlockDataSet *input1(vtkMultiBlockDataSet::SafeDownCast(inputInfo->Get(vtkDataObject::DATA_OBJECT())));
128 throw INTERP_KERNEL::Exception("Input dataSet must be a DataSet or single elt multi block dataset expected !");
129 if(input1->GetNumberOfBlocks()!=1)
130 throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with not exactly one block ! Use MergeBlocks or ExtractBlocks filter before calling this filter !");
131 vtkDataObject *input2(input1->GetBlock(0));
133 throw INTERP_KERNEL::Exception("Input dataSet is a multiblock dataset with exactly one block but this single element is NULL !");
134 vtkDataSet *input2c(vtkDataSet::SafeDownCast(input2));
136 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 !");
140 throw INTERP_KERNEL::Exception("Input data set is NULL !");
141 usgIn=vtkUnstructuredGrid::SafeDownCast(input);
143 throw INTERP_KERNEL::Exception("Input data set is not an unstructured mesh ! This filter works only on unstructured meshes !");
146 vtkGaussToCell::vtkGaussToCell():avgStatus(true),maxStatus(false),minStatus(false)
148 this->SetNumberOfInputPorts(1);
149 this->SetNumberOfOutputPorts(1);
152 vtkGaussToCell::~vtkGaussToCell()
156 void vtkGaussToCell::SetAvgFlag(bool avgStatus)
158 if(this->avgStatus!=avgStatus)
160 this->avgStatus=avgStatus;
165 void vtkGaussToCell::SetMaxFlag(bool maxStatus)
167 if(this->maxStatus!=maxStatus)
169 this->maxStatus=maxStatus;
174 void vtkGaussToCell::SetMinFlag(bool minStatus)
176 if(this->minStatus!=minStatus)
178 this->minStatus=minStatus;
183 int vtkGaussToCell::RequestInformation(vtkInformation * /*request*/, vtkInformationVector **inputVector, vtkInformationVector * /*outputVector*/)
185 //std::cerr << "########################################## vtkGaussToCell::RequestInformation ##########################################" << std::endl;
188 vtkUnstructuredGrid *usgIn(0);
189 ExtractInfo(inputVector[0],usgIn);
191 catch(INTERP_KERNEL::Exception& e)
193 std::ostringstream oss;
194 oss << "Exception has been thrown in vtkGaussToCell::RequestInformation : " << e.what() << std::endl;
195 if(this->HasObserver("ErrorEvent") )
196 this->InvokeEvent("ErrorEvent",const_cast<char *>(oss.str().c_str()));
198 vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
199 vtkObject::BreakOnError();
205 typedef void (*DataComputer) (const double *inData, const vtkIdType *offData, const std::vector<int> *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData);
207 void ComputeAvg(const double *inData, const vtkIdType *offData, const std::vector<int> *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData)
209 std::fill(outData,outData+outNbCells*zeNbCompo,0.);
210 for(auto i=0;i<outNbCells;i++,outData+=zeNbCompo)
212 auto NbGaussPt((*nbgPerCell)[i]);
213 for(auto j=0;j<NbGaussPt;j++)
214 std::transform(inData+(offData[i]+j)*zeNbCompo,inData+(offData[i]+(j+1))*zeNbCompo,outData,outData,[](const double& v, const double& w) { return v+w; });
215 double deno(1./(double)NbGaussPt);
216 std::transform(outData,outData+zeNbCompo,outData,[deno](const double& v) { return v*deno; });
220 void ComputeMax(const double *inData, const vtkIdType *offData, const std::vector<int> *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData)
222 std::fill(outData,outData+outNbCells*zeNbCompo,-std::numeric_limits<double>::max());
223 for(auto i=0;i<outNbCells;i++,outData+=zeNbCompo)
225 auto NbGaussPt((*nbgPerCell)[i]);
226 for(auto j=0;j<NbGaussPt;j++)
227 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); });
231 void ComputeMin(const double *inData, const vtkIdType *offData, const std::vector<int> *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData)
233 std::fill(outData,outData+outNbCells*zeNbCompo,std::numeric_limits<double>::max());
234 for(auto i=0;i<outNbCells;i++,outData+=zeNbCompo)
236 auto NbGaussPt((*nbgPerCell)[i]);
237 for(auto j=0;j<NbGaussPt;j++)
238 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); });
242 void DealWith(const char *postName, vtkDoubleArray *zearray, vtkIdTypeArray *offsets, const std::vector<int> *nbgPerCell, vtkIdType outNbCells, vtkDoubleArray *zeOutArray, DataComputer dc)
244 std::ostringstream oss;
245 oss << zearray->GetName() << '_' << postName;
246 int zeNbCompo(zearray->GetNumberOfComponents());
248 std::string st(oss.str());
249 zeOutArray->SetName(st.c_str());
251 zeOutArray->SetNumberOfComponents(zeNbCompo);
252 zeOutArray->SetNumberOfTuples(outNbCells);
253 double *outData(zeOutArray->GetPointer(0));
254 const double *inData(zearray->GetPointer(0));
255 const auto *offData(offsets->GetPointer(0));
256 for(auto i=0;i<zeNbCompo;i++)
258 const char *comp(zearray->GetComponentName(i));
260 zeOutArray->SetComponentName(i,comp);
262 dc(inData,offData,nbgPerCell,zeNbCompo,outNbCells,outData);
265 int vtkGaussToCell::RequestData(vtkInformation * /*request*/, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
267 //std::cerr << "########################################## vtkGaussToCell::RequestData ##########################################" << std::endl;
270 std::vector<double> GaussAdvData;
271 bool isOK(IsInformationOK(inputVector[0]->GetInformationObject(0),GaussAdvData));
273 throw INTERP_KERNEL::Exception("Sorry but no advanced gauss info found ! Expect to be called right after a MEDReader containing Gauss Points !");
274 vtkUnstructuredGrid *usgIn(0);
275 ExtractInfo(inputVector[0],usgIn);
276 vtkInformation *outInfo(outputVector->GetInformationObject(0));
277 vtkUnstructuredGrid *output(vtkUnstructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())));
278 output->ShallowCopy(usgIn);
280 std::string zeArrOffset;
281 int nArrays(usgIn->GetFieldData()->GetNumberOfArrays());
282 std::map<vtkIdTypeArray *,std::vector<int> > offsetKeyMap;//Map storing for each offsets array the corresponding nb of Gauss Points per cell
283 for(int i=0;i<nArrays;i++)
285 vtkDataArray *array(usgIn->GetFieldData()->GetArray(i));
288 const char* arrayOffsetName(array->GetInformation()->Get(vtkQuadratureSchemeDefinition::QUADRATURE_OFFSET_ARRAY_NAME()));
291 std::string arrOffsetNameCpp(arrayOffsetName);
292 if(arrOffsetNameCpp.find("ELGA@")==std::string::npos)
294 if(zeArrOffset.empty())
295 zeArrOffset=arrOffsetNameCpp;
297 if(zeArrOffset!=arrOffsetNameCpp)
299 throw INTERP_KERNEL::Exception("ComputeGaussToCell : error in QUADRATURE_OFFSET_ARRAY_NAME for Gauss fields array !");
301 vtkDataArray *offTmp(usgIn->GetCellData()->GetArray(zeArrOffset.c_str()));
304 std::ostringstream oss; oss << "ComputeGaussToCell : cell field " << zeArrOffset << " not found !";
305 throw INTERP_KERNEL::Exception(oss.str());
307 vtkIdTypeArray *offsets(vtkIdTypeArray::SafeDownCast(offTmp));
310 std::ostringstream oss; oss << "ComputeGaussToCell : cell field " << zeArrOffset << " exists but not with the right type of data !";
311 throw INTERP_KERNEL::Exception(oss.str());
313 vtkDoubleArray *zearray(vtkDoubleArray::SafeDownCast(array));
317 std::map<vtkIdTypeArray *,std::vector<int> >::iterator nbgPerCellPt(offsetKeyMap.find(offsets));
318 const std::vector<int> *nbgPerCell(nullptr);
319 if(nbgPerCellPt==offsetKeyMap.end())
322 vtkInformation *info(offsets->GetInformation());
324 throw INTERP_KERNEL::Exception("info is null ! Internal error ! Looks bad !");
325 vtkInformationQuadratureSchemeDefinitionVectorKey *key(vtkQuadratureSchemeDefinition::DICTIONARY());
327 throw INTERP_KERNEL::Exception("No quadrature key in info included in offets array ! Internal error ! Looks bad !");
328 int dictSize(key->Size(info));
329 INTERP_KERNEL::AutoPtr<vtkQuadratureSchemeDefinition *> dict(new vtkQuadratureSchemeDefinition *[dictSize]);
330 key->GetRange(info,dict,0,0,dictSize);
331 auto nbOfCells(output->GetNumberOfCells());
332 std::vector<int> nbg(nbOfCells);
333 for(auto cellId=0;cellId<nbOfCells;cellId++)
335 int ct(output->GetCellType(cellId));
336 vtkQuadratureSchemeDefinition *gaussLoc(dict[ct]);
339 std::ostringstream oss; oss << "For cell " << cellId << " no Gauss info attached !";
340 throw INTERP_KERNEL::Exception(oss.str());
342 int np(gaussLoc->GetNumberOfQuadraturePoints());
345 nbgPerCell=&((*(offsetKeyMap.emplace(offsets,std::move(nbg)).first)).second);
349 nbgPerCell=&((*nbgPerCellPt).second);
351 auto outNbCells(nbgPerCell->size());
354 vtkNew<vtkDoubleArray> zeOutArray;
355 DealWith("avg",zearray,offsets,nbgPerCell,outNbCells,zeOutArray,ComputeAvg);
356 output->GetCellData()->AddArray(zeOutArray);
360 vtkNew<vtkDoubleArray> zeOutArray;
361 DealWith("max",zearray,offsets,nbgPerCell,outNbCells,zeOutArray,ComputeMax);
362 output->GetCellData()->AddArray(zeOutArray);
366 vtkNew<vtkDoubleArray> zeOutArray;
367 DealWith("min",zearray,offsets,nbgPerCell,outNbCells,zeOutArray,ComputeMin);
368 output->GetCellData()->AddArray(zeOutArray);
372 catch(INTERP_KERNEL::Exception& e)
374 std::ostringstream oss;
375 oss << "Exception has been thrown in vtkGaussToCell::RequestData : " << e.what() << std::endl;
376 if(this->HasObserver("ErrorEvent") )
377 this->InvokeEvent("ErrorEvent",const_cast<char *>(oss.str().c_str()));
379 vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
380 vtkObject::BreakOnError();
386 void vtkGaussToCell::PrintSelf(ostream& os, vtkIndent indent)
388 this->Superclass::PrintSelf(os, indent);