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 }