src/Plugins/MEDReader/plugin/MEDReaderIO/vtkELNOMeshFilter.cxx

   1 // Copyright (C) 2010-2022  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
  20 #include "vtkELNOMeshFilter.h"
  21 #include "vtkInformation.h"
  22 #include "vtkInformationVector.h"
  23 #include "vtkInformationIntegerKey.h"
  24 #include "vtkObjectFactory.h"
  25 #include "vtkPolyDataAlgorithm.h"
  26 #include "vtkPolyData.h"
  27 #include "vtkIdTypeArray.h"
  28 #include "vtkQuadratureSchemeDefinition.h"
  29 #include "vtkInformationQuadratureSchemeDefinitionVectorKey.h"
  30 #include "vtkUnstructuredGrid.h"
  31 #include "vtkShrinkFilter.h"
  32 #include "vtkSmartPointer.h"
  33 #include "vtkPointData.h"
  34 #include "vtkCellData.h"
  35 #include "vtkIdList.h"
  36 #include "vtkCell.h"
  37
  38 #include "MEDUtilities.hxx"
  39
  40 #include <map>
  41
  42 vtkStandardNewMacro(vtkELNOMeshFilter)
  43
  44 vtkELNOMeshFilter::vtkELNOMeshFilter():ShrinkFactor(0.9999)
  45 {
  46 }
  47
  48 vtkELNOMeshFilter::~vtkELNOMeshFilter()
  49 {
  50 }
  51
  52 int vtkELNOMeshFilter::RequestData(vtkInformation * /*request*/,
  53     vtkInformationVector **input, vtkInformationVector *output)
  54 {
  55   vtkUnstructuredGrid *usgIn = vtkUnstructuredGrid::SafeDownCast(
  56       input[0]->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT()));
  57
  58   vtkUnstructuredGrid *usgOut = vtkUnstructuredGrid::SafeDownCast(
  59       output->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT()));
  60
  61   if(usgIn == NULL || usgOut == NULL)
  62     {
  63       vtkDebugMacro("vtkELNOMeshFilter not correctly configured : Invalid input or output !");
  64       return 0;
  65     }
  66
  67   // creates offsets array
  68
  69   // first shrink the input
  70   vtkUnstructuredGrid* usgInClone = usgIn->NewInstance();
  71   usgInClone->ShallowCopy(usgIn);
  72   vtkSmartPointer<vtkShrinkFilter> shrink(vtkSmartPointer<vtkShrinkFilter>::New());
  73   shrink->SetInputData(usgInClone);
  74   shrink->SetShrinkFactor(this->ShrinkFactor);
  75   shrink->Update();
  76   vtkUnstructuredGrid *shrinked(shrink->GetOutput());
  77   usgInClone->Delete();
  78   usgOut->ShallowCopy(shrinked);
  79   // OK for the output
  80
  81   // now copy ELNO data. Start by verifying if it is possible to
  82   // shallow copy the array.
  83   //vtkInformation *info(usgIn->GetInformation()); // todo: unused
  84   //
  85   vtkIdType nVerts(shrinked->GetNumberOfPoints()),ncell(usgIn->GetNumberOfCells());
  86   // first loop through all cells to check if a shallow copy is possible
  87   bool shallowok(true);// Anthony : checks that shrink works well. Really necessary ?
  88   vtkIdType previous(0),offset(0);
  89
  90   for(vtkIdType cellId = 0; cellId < ncell; cellId++)
  91     {
  92       if(offset != previous)
  93         {
  94           shallowok = false;
  95           break;
  96         }
  97       vtkCell *cell(usgIn->GetCell(cellId));
  98       vtkIdType nbptsInCell(cell->GetNumberOfPoints());
  99       previous = offset + nbptsInCell;
 100       //
 101       offset += nbptsInCell ;
 102     }
 103   //
 104   if(shallowok)
 105     shallowok = (previous == nVerts);
 106
 107   vtkFieldData *fielddata(usgIn->GetFieldData());
 108   for(int index = 0; index < fielddata->GetNumberOfArrays(); index++)
 109     {
 110       vtkDataArray *data(fielddata->GetArray(index));
 111       vtkQuadratureSchemeDefinition **dict = 0;
 112       vtkInformationQuadratureSchemeDefinitionVectorKey *key(vtkQuadratureSchemeDefinition::DICTIONARY());
 113       if(key->Has(data->GetInformation()))
 114         {
 115           int dictSize(key->Size(data->GetInformation()));
 116           dict=new vtkQuadratureSchemeDefinition *[dictSize];
 117           key->GetRange(data->GetInformation(),dict,0,0,dictSize);
 118         }
 119       if(data == NULL)
 120         {
 121           delete [] dict;
 122           continue;
 123         }
 124
 125       vtkInformation *info(data->GetInformation());
 126       const char *arrayOffsetName(info->Get(vtkQuadratureSchemeDefinition::QUADRATURE_OFFSET_ARRAY_NAME()));
 127       vtkIdTypeArray *offData(0);
 128       bool isELGA(false),isELNO(false);
 129
 130       if(arrayOffsetName)
 131         {
 132           vtkFieldData *cellData(usgIn->GetCellData());
 133           vtkDataArray *offDataTmp(cellData->GetArray(arrayOffsetName));
 134           isELGA=offDataTmp->GetInformation()->Get(MEDUtilities::ELGA())==1;
 135           isELNO=offDataTmp->GetInformation()->Get(MEDUtilities::ELNO())==1;
 136           offData=dynamic_cast<vtkIdTypeArray *>(offDataTmp);
 137         }
 138
 139       if(arrayOffsetName == NULL || isELGA)
 140         {
 141           if(shallowok && data->GetNumberOfTuples()==nVerts )// Anthony : is it not a little confusing to assign a FieldData on Points because the number of tuples fits the number of nodes of shrinked mesh ?
 142             usgOut->GetPointData()->AddArray(data);
 143           else
 144             shrinked->GetFieldData()->AddArray(data);
 145           delete [] dict;
 146           continue;
 147         }
 148       else
 149         {
 150           vtkDataArray *newArray(data->NewInstance());
 151           newArray->SetName(data->GetName());
 152           usgOut->GetPointData()->AddArray(newArray);
 153           newArray->SetNumberOfComponents(data->GetNumberOfComponents());
 154           newArray->SetNumberOfTuples(usgOut->GetNumberOfPoints());
 155           newArray->CopyComponentNames(data);
 156           newArray->Delete();
 157           if(isELGA)
 158             {
 159               vtkIdList *ids(vtkIdList::New());
 160               vtkIdType offset(0);
 161               for(vtkIdType cellId=0;cellId<ncell;cellId++)
 162                 {
 163                   int cellType(shrinked->GetCellType(cellId));
 164                   shrinked->GetCellPoints(cellId, ids);
 165                   for(int id = 0; id < dict[cellType]->GetNumberOfQuadraturePoints(); id++)
 166                     {
 167                       const double * w = dict[cellType]->GetShapeFunctionWeights(id);
 168                       int j;
 169                       for(j = 0; j < dict[cellType]->GetNumberOfNodes(); j++)
 170                         {
 171                           if(w[j] == 1.0)
 172                             break;
 173                         }
 174                       if(j == dict[cellType]->GetNumberOfNodes())
 175                         {
 176                           j = id;
 177                         }
 178                       newArray->SetTuple(ids->GetId(id), offset + j, data);
 179                     }
 180                   vtkCell *cell(usgIn->GetCell(cellId));
 181                   vtkIdType nbptsInCell(cell->GetNumberOfPoints());
 182                   offset+=nbptsInCell;
 183                 }
 184               ids->FastDelete();
 185             }
 186           else if(offData && isELNO)
 187             {
 188               vtkIdType *offsetPtr(offData->GetPointer(0));
 189               vtkIdType zeId(0);
 190               for(vtkIdType cellId=0;cellId<ncell;cellId++)
 191                 {
 192                   vtkCell *cell(shrinked->GetCell(cellId));
 193                   vtkIdType nbPoints(cell->GetNumberOfPoints())/*,offset(offsetPtr[cellId])*/; // todo: offset is unused
 194                   for(vtkIdType j=0;j<nbPoints;j++,zeId++)
 195                     newArray->SetTuple(zeId,offsetPtr[cellId]+j,data);
 196                 }
 197             }
 198           else
 199             {
 200               delete [] dict;
 201               continue ;
 202             }
 203         }
 204       delete [] dict;
 205     }
 206   return 1;
 207 }
 208
 209 void vtkELNOMeshFilter::PrintSelf(ostream& os, vtkIndent indent)
 210 {
 211   this->Superclass::PrintSelf(os, indent);
 212 }