src/Plugins/MEDReader/plugin/MEDReaderIO/vtkELNOFilter.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 "vtkELNOFilter.h"
  21 #include "vtkInformation.h"
  22 #include "vtkInformationVector.h"
  23 #include "vtkObjectFactory.h"
  24 #include "vtkPolyDataAlgorithm.h"
  25 #include "vtkPolyData.h"
  26 #include "vtkIdTypeArray.h"
  27 #include "vtkFieldData.h"
  28 #include "vtkCellData.h"
  29 #include "vtkPointData.h"
  30 #include "vtkCell.h"
  31 #include "vtkInformationQuadratureSchemeDefinitionVectorKey.h"
  32 #include "vtkQuadratureSchemeDefinition.h"
  33 #include "vtkUnstructuredGrid.h"
  34
  35 #include "MEDUtilities.hxx"
  36 #include "InterpKernelAutoPtr.hxx"
  37
  38 //vtkCxxRevisionMacro(vtkELNOFilter, "$Revision: 1.2.2.2 $");
  39 vtkStandardNewMacro(vtkELNOFilter)
  40
  41 vtkELNOFilter::vtkELNOFilter()
  42 {
  43   this->ShrinkFactor = 0.5;
  44 }
  45
  46 vtkELNOFilter::~vtkELNOFilter()
  47 {
  48 }
  49
  50 int vtkELNOFilter::RequestData(vtkInformation *request, vtkInformationVector **input, vtkInformationVector *output)
  51 {
  52   vtkUnstructuredGrid *usgIn(vtkUnstructuredGrid::SafeDownCast( input[0]->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT())));
  53   vtkPolyData *pdOut(vtkPolyData::SafeDownCast(output->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT())));
  54
  55   //vtkDataArray *array(this->GetInputArrayToProcess(0, input)); // todo: unused
  56   vtkIdTypeArray* offsets(vtkIdTypeArray::SafeDownCast(this->GetInputArrayToProcess(0, input)));
  57
  58   if(usgIn == NULL || offsets == NULL || pdOut == NULL)
  59     {
  60       vtkDebugMacro("vtkELNOFilter no correctly configured : offsets = " << offsets);
  61       return 1;
  62     }
  63
  64   vtkInformation *info(offsets->GetInformation());
  65   vtkInformationQuadratureSchemeDefinitionVectorKey *key(vtkQuadratureSchemeDefinition::DICTIONARY());
  66   if(!key->Has(info))
  67     {
  68       vtkDebugMacro("Dictionary is not present in array " << offsets->GetName() << " " << offsets << " Aborting." );
  69       return 1;
  70     }
  71
  72   int res(this->Superclass::RequestData(request, input, output));
  73   if(res == 0)
  74     return 0;
  75
  76   int dictSize(key->Size(info));
  77   vtkQuadratureSchemeDefinition **dict = new vtkQuadratureSchemeDefinition *[dictSize];
  78   key->GetRange(info, dict, 0, 0, dictSize);
  79
  80   vtkIdType ncell(usgIn->GetNumberOfCells());
  81   vtkPoints *points(pdOut->GetPoints());
  82   vtkIdType start(0);
  83   for(vtkIdType cellId = 0; cellId < ncell; cellId++)
  84     {
  85       //vtkIdType offset(offsets->GetValue(cellId)); // todo: unused
  86       int cellType(usgIn->GetCellType(cellId));
  87       // a simple check to see if a scheme really exists for this cell type.
  88       // should not happen if the cell type has not been modified.
  89       if(dict[cellType] == NULL)
  90         continue;
  91       int np = dict[cellType]->GetNumberOfQuadraturePoints();
  92       double center[3] = {0, 0, 0};
  93       for(int id = start; id < start + np; id++)
  94         {
  95           double *position = points->GetPoint(id);
  96           center[0] += position[0];
  97           center[1] += position[1];
  98           center[2] += position[2];
  99         }
 100       center[0] /= np;
 101       center[1] /= np;
 102       center[2] /= np;
 103       for(int id = start; id < start + np; id++)
 104         {
 105           double *position = points->GetPoint(id);
 106           double newpos[3];
 107           newpos[0] = position[0] * this->ShrinkFactor + center[0] * (1 - this->ShrinkFactor);
 108           newpos[1] = position[1] * this->ShrinkFactor + center[1] * (1 - this->ShrinkFactor);
 109           newpos[2] = position[2] * this->ShrinkFactor + center[2] * (1 - this->ShrinkFactor);
 110           points->SetPoint(id, newpos);
 111         }
 112       start += np;
 113     }
 114   //// bug EDF 8407 PARAVIS - mantis 22610
 115   vtkFieldData *fielddata(usgIn->GetFieldData());
 116   for(int index=0;index<fielddata->GetNumberOfArrays();index++)
 117     {
 118       vtkDataArray *data(fielddata->GetArray(index));
 119       vtkInformation *info(data->GetInformation());
 120       const char *arrayOffsetName(info->Get(vtkQuadratureSchemeDefinition::QUADRATURE_OFFSET_ARRAY_NAME()));
 121       vtkIdTypeArray *offData(0);
 122       bool isELNO(false);
 123       if(arrayOffsetName)
 124         {
 125           vtkCellData *cellData(usgIn->GetCellData());
 126           vtkDataArray *offDataTmp(cellData->GetArray(arrayOffsetName));
 127           isELNO=offDataTmp->GetInformation()->Get(MEDUtilities::ELNO())==1;
 128           offData=dynamic_cast<vtkIdTypeArray *>(offDataTmp);
 129         }
 130       if(isELNO && offData)
 131         {
 132           vtkIdType nbCellsInput(usgIn->GetNumberOfCells());
 133           if(nbCellsInput==0)
 134             continue ;//no cells -> no fields
 135           // First trying to detected if we are in the special case where data can be used directly. To detect that look at offData. If offData.front()==0 && offData->back()+NbOfNodesInLastCell==data->GetNumberOfTuples() OK.
 136           vtkCell *cell(usgIn->GetCell(nbCellsInput-1));
 137           bool statement0(offData->GetTuple1(0)==0);
 138           bool statement1(offData->GetTuple1(nbCellsInput-1)+cell->GetNumberOfPoints()==data->GetNumberOfTuples());
 139           if(statement0 && statement1)
 140             pdOut->GetPointData()->AddArray(data);//We are lucky ! No extraction needed.
 141           else
 142             {//not lucky ! Extract values from data. A previous threshold has been done... Bug EDF8662
 143               vtkDataArray *newArray(data->NewInstance());
 144               newArray->SetName(data->GetName());
 145               pdOut->GetPointData()->AddArray(newArray);
 146               newArray->SetNumberOfComponents(data->GetNumberOfComponents());
 147               newArray->SetNumberOfTuples(pdOut->GetNumberOfPoints());
 148               newArray->CopyComponentNames(data);
 149               newArray->Delete();
 150               vtkIdType *offsetPtr(offData->GetPointer(0));
 151               vtkIdType zeId(0);
 152               for(vtkIdType cellId=0;cellId<nbCellsInput;cellId++)
 153                 {
 154                   vtkCell *cell(usgIn->GetCell(cellId));
 155                   vtkIdType nbPoints(cell->GetNumberOfPoints())/*,offset(offsetPtr[cellId])*/; // todo: offset is unused
 156                   for(vtkIdType j=0;j<nbPoints;j++,zeId++)
 157                     newArray->SetTuple(zeId,offsetPtr[cellId]+j,data);
 158                 }
 159             }
 160         }
 161     }
 162   AttachCellFieldsOn(usgIn,pdOut->GetCellData(),pdOut->GetNumberOfCells());
 163   return 1;
 164 }
 165
 166 void vtkELNOFilter::PrintSelf(ostream& os, vtkIndent indent)
 167 {
 168   this->Superclass::PrintSelf(os, indent);
 169   os << indent << "ShrinkFactor : " << this->ShrinkFactor << endl;
 170 }
 171
 172 /*!
 173  * This method attach fields on cell of \a inGrid and add it as a point data in \a outData.
 174  */
 175 void vtkELNOFilter::AttachCellFieldsOn(vtkUnstructuredGrid *inGrid, vtkCellData *outData, int nbCellsOut)
 176 {
 177   vtkCellData *cd(inGrid->GetCellData());
 178   int nbOfArrays(cd->GetNumberOfArrays());
 179   vtkIdType nbCells(inGrid->GetNumberOfCells());
 180   if(nbOfArrays==0)
 181     return ;
 182   INTERP_KERNEL::AutoPtr<vtkIdType> tmpPtr(new vtkIdType[nbCells]);
 183   for(vtkIdType cellId=0;cellId<nbCells;cellId++)
 184     {
 185       vtkCell *cell(inGrid->GetCell(cellId));
 186       tmpPtr[cellId]=cell->GetNumberOfPoints();
 187     }
 188   for(int index=0;index<nbOfArrays;index++)
 189     {
 190       vtkDataArray *data(cd->GetArray(index));
 191       vtkDataArray *newArray(data->NewInstance());
 192       newArray->SetName(data->GetName());
 193       outData->AddArray(newArray);
 194       newArray->SetNumberOfComponents(data->GetNumberOfComponents());
 195       newArray->SetNumberOfTuples(nbCellsOut);
 196       newArray->CopyComponentNames(data);
 197       newArray->Delete();
 198       vtkIdType offset(0);
 199       for(vtkIdType cellId=0;cellId<nbCells;cellId++)
 200         {
 201           for(vtkIdType j=0;j<tmpPtr[cellId];j++,offset++)
 202             newArray->SetTuple(offset,cellId,data);
 203         }
 204     }
 205 }