src/Plugins/GaussToCell/plugin/GaussToCellModule/vtkGaussToCell.cxx

   1 // Copyright (C) 2018-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 "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::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;
  88
  89 vtkStandardNewMacro(vtkGaussToCell)
  90
  91 vtkInformationDoubleVectorKey *GetMEDReaderMetaDataIfAny()
  92 {
  93   static const char ZE_KEY[]="vtkMEDReader::GAUSS_DATA";
  94   MEDCoupling::GlobalDict *gd(MEDCoupling::GlobalDict::GetInstance());
  95   if(!gd->hasKey(ZE_KEY))
  96     return 0;
  97   std::string ptSt(gd->value(ZE_KEY));
  98   void *pt(0);
  99   std::istringstream iss(ptSt); iss >> pt;
 100   return reinterpret_cast<vtkInformationDoubleVectorKey *>(pt);
 101 }
 102
 103 bool IsInformationOK(vtkInformation *info, std::vector<double>& data)
 104 {
 105   vtkInformationDoubleVectorKey *key(GetMEDReaderMetaDataIfAny());
 106   if(!key)
 107     return false;
 108   // Check the information contain meta data key
 109   if(!info->Has(key))
 110     return false;
 111   int lgth(key->Length(info));
 112   const double *data2(info->Get(key));
 113   data.insert(data.end(),data2,data2+lgth);
 114   return true;
 115 }
 116
 117 void ExtractInfo(vtkInformationVector *inputVector, vtkUnstructuredGrid *& usgIn)
 118 {
 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())));
 123   if(input0)
 124     input=input0;
 125   else
 126     {
 127       if(!input1)
 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));
 132       if(!input2)
 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));
 135       if(!input2c)
 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 !");
 137       input=input2c;
 138     }
 139   if(!input)
 140     throw INTERP_KERNEL::Exception("Input data set is NULL !");
 141   usgIn=vtkUnstructuredGrid::SafeDownCast(input);
 142   if(!usgIn)
 143     throw INTERP_KERNEL::Exception("Input data set is not an unstructured mesh ! This filter works only on unstructured meshes !");
 144 }
 145
 146 vtkGaussToCell::vtkGaussToCell():avgStatus(true),maxStatus(false),minStatus(false)
 147 {
 148   this->SetNumberOfInputPorts(1);
 149   this->SetNumberOfOutputPorts(1);
 150 }
 151
 152 vtkGaussToCell::~vtkGaussToCell()
 153 {
 154 }
 155
 156 void vtkGaussToCell::SetAvgFlag(bool avgStatus)
 157 {
 158   if(this->avgStatus!=avgStatus)
 159     {
 160       this->avgStatus=avgStatus;
 161       this->Modified();
 162     }
 163 }
 164
 165 void vtkGaussToCell::SetMaxFlag(bool maxStatus)
 166 {
 167   if(this->maxStatus!=maxStatus)
 168     {
 169       this->maxStatus=maxStatus;
 170       this->Modified();
 171     }
 172 }
 173
 174 void vtkGaussToCell::SetMinFlag(bool minStatus)
 175 {
 176   if(this->minStatus!=minStatus)
 177     {
 178       this->minStatus=minStatus;
 179       this->Modified();
 180     }
 181 }
 182
 183 int vtkGaussToCell::RequestInformation(vtkInformation * /*request*/, vtkInformationVector **inputVector, vtkInformationVector * /*outputVector*/)
 184 {
 185   //std::cerr << "########################################## vtkGaussToCell::RequestInformation ##########################################" << std::endl;
 186   try
 187     {
 188       vtkUnstructuredGrid *usgIn(0);
 189       ExtractInfo(inputVector[0],usgIn);
 190     }
 191   catch(INTERP_KERNEL::Exception& e)
 192     {
 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()));
 197       else
 198         vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
 199       vtkObject::BreakOnError();
 200       return 0;
 201     }
 202   return 1;
 203 }
 204
 205 typedef void (*DataComputer) (const double *inData, const vtkIdType *offData, const std::vector<int> *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData);
 206
 207 void ComputeAvg(const double *inData, const vtkIdType *offData, const std::vector<int> *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData)
 208 {
 209   std::fill(outData,outData+outNbCells*zeNbCompo,0.);
 210   for(auto i=0;i<outNbCells;i++,outData+=zeNbCompo)
 211     {
 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; });
 217     }
 218 }
 219
 220 void ComputeMax(const double *inData, const vtkIdType *offData, const std::vector<int> *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData)
 221 {
 222   std::fill(outData,outData+outNbCells*zeNbCompo,-std::numeric_limits<double>::max());
 223   for(auto i=0;i<outNbCells;i++,outData+=zeNbCompo)
 224     {
 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); });
 228     }
 229 }
 230
 231 void ComputeMin(const double *inData, const vtkIdType *offData, const std::vector<int> *nbgPerCell, int zeNbCompo, vtkIdType outNbCells, double *outData)
 232 {
 233   std::fill(outData,outData+outNbCells*zeNbCompo,std::numeric_limits<double>::max());
 234   for(auto i=0;i<outNbCells;i++,outData+=zeNbCompo)
 235     {
 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); });
 239     }
 240 }
 241
 242 void DealWith(const char *postName, vtkDoubleArray *zearray, vtkIdTypeArray *offsets, const std::vector<int> *nbgPerCell, vtkIdType outNbCells, vtkDoubleArray *zeOutArray, DataComputer dc)
 243 {
 244   std::ostringstream oss;
 245   oss << zearray->GetName() << '_' << postName;
 246   int zeNbCompo(zearray->GetNumberOfComponents());
 247   {
 248     std::string st(oss.str());
 249     zeOutArray->SetName(st.c_str());
 250   }
 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++)
 257     {
 258       const char *comp(zearray->GetComponentName(i));
 259       if(comp)
 260         zeOutArray->SetComponentName(i,comp);
 261     }
 262   dc(inData,offData,nbgPerCell,zeNbCompo,outNbCells,outData);
 263 }
 264
 265 int vtkGaussToCell::RequestData(vtkInformation * /*request*/, vtkInformationVector **inputVector, vtkInformationVector *outputVector)
 266 {
 267   //std::cerr << "########################################## vtkGaussToCell::RequestData        ##########################################" << std::endl;
 268   try
 269     {
 270       std::vector<double> GaussAdvData;
 271       bool isOK(IsInformationOK(inputVector[0]->GetInformationObject(0),GaussAdvData));
 272       if(!isOK)
 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);
 279       //
 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++)
 284         {
 285           vtkDataArray *array(usgIn->GetFieldData()->GetArray(i));
 286           if(!array)
 287             continue;
 288           const char* arrayOffsetName(array->GetInformation()->Get(vtkQuadratureSchemeDefinition::QUADRATURE_OFFSET_ARRAY_NAME()));
 289           if(!arrayOffsetName)
 290             continue;
 291           std::string arrOffsetNameCpp(arrayOffsetName);
 292           if(arrOffsetNameCpp.find("ELGA@")==std::string::npos)
 293             continue;
 294           if(zeArrOffset.empty())
 295             zeArrOffset=arrOffsetNameCpp;
 296           else
 297             if(zeArrOffset!=arrOffsetNameCpp)
 298               {
 299                 throw INTERP_KERNEL::Exception("ComputeGaussToCell : error in QUADRATURE_OFFSET_ARRAY_NAME for Gauss fields array !");
 300               }
 301           vtkDataArray *offTmp(usgIn->GetCellData()->GetArray(zeArrOffset.c_str()));
 302           if(!offTmp)
 303             {
 304               std::ostringstream oss; oss << "ComputeGaussToCell : cell field " << zeArrOffset << " not found !";
 305               throw INTERP_KERNEL::Exception(oss.str());
 306             }
 307           vtkIdTypeArray *offsets(vtkIdTypeArray::SafeDownCast(offTmp));
 308           if(!offsets)
 309             {
 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());
 312             }
 313           vtkDoubleArray *zearray(vtkDoubleArray::SafeDownCast(array));
 314           if(!zearray)
 315             continue ;
 316           //
 317           std::map<vtkIdTypeArray *,std::vector<int> >::iterator nbgPerCellPt(offsetKeyMap.find(offsets));
 318           const std::vector<int> *nbgPerCell(nullptr);
 319           if(nbgPerCellPt==offsetKeyMap.end())
 320             {
 321               // fini la parlote
 322               vtkInformation *info(offsets->GetInformation());
 323               if(!info)
 324                 throw INTERP_KERNEL::Exception("info is null ! Internal error ! Looks bad !");
 325               vtkInformationQuadratureSchemeDefinitionVectorKey *key(vtkQuadratureSchemeDefinition::DICTIONARY());
 326               if(!key->Has(info))
 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++)
 334                 {
 335                   int ct(output->GetCellType(cellId));
 336                   vtkQuadratureSchemeDefinition *gaussLoc(dict[ct]);
 337                   if(!gaussLoc)
 338                     {
 339                       std::ostringstream oss; oss << "For cell " << cellId << " no Gauss info attached !";
 340                       throw INTERP_KERNEL::Exception(oss.str());
 341                     }
 342                   int np(gaussLoc->GetNumberOfQuadraturePoints());
 343                   nbg[cellId]=np;
 344                 }
 345               nbgPerCell=&((*(offsetKeyMap.emplace(offsets,std::move(nbg)).first)).second);
 346             }
 347           else
 348             {
 349               nbgPerCell=&((*nbgPerCellPt).second);
 350             }
 351           auto outNbCells(nbgPerCell->size());
 352           if(this->avgStatus)
 353             {
 354               vtkNew<vtkDoubleArray> zeOutArray;
 355               DealWith("avg",zearray,offsets,nbgPerCell,outNbCells,zeOutArray,ComputeAvg);
 356               output->GetCellData()->AddArray(zeOutArray);
 357             }
 358           if(this->maxStatus)
 359             {
 360               vtkNew<vtkDoubleArray> zeOutArray;
 361               DealWith("max",zearray,offsets,nbgPerCell,outNbCells,zeOutArray,ComputeMax);
 362               output->GetCellData()->AddArray(zeOutArray);
 363             }
 364           if(this->minStatus)
 365             {
 366               vtkNew<vtkDoubleArray> zeOutArray;
 367               DealWith("min",zearray,offsets,nbgPerCell,outNbCells,zeOutArray,ComputeMin);
 368               output->GetCellData()->AddArray(zeOutArray);
 369             }
 370         }
 371     }
 372   catch(INTERP_KERNEL::Exception& e)
 373     {
 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()));
 378       else
 379         vtkOutputWindowDisplayErrorText(const_cast<char *>(oss.str().c_str()));
 380       vtkObject::BreakOnError();
 381       return 0;
 382     }
 383   return 1;
 384 }
 385
 386 void vtkGaussToCell::PrintSelf(ostream& os, vtkIndent indent)
 387 {
 388   this->Superclass::PrintSelf(os, indent);
 389 }