From 31c4fb7eb65d19f8b21266c38a0460e37fc9d39a Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Thu, 18 Dec 2014 15:42:29 +0100 Subject: [PATCH] Correction of issue EDF8662 (ELNOPoint, ELNOMesh after ExtractGroup) --- src/Plugins/MEDReader/IO/vtkELNOFilter.cxx | 145 ++++++++++-------- .../MEDReader/IO/vtkELNOMeshFilter.cxx | 73 +++++---- src/Plugins/MEDReader/Test/CMakeLists.txt | 2 + src/Plugins/MEDReader/Test/testMEDReader16.py | 103 +++++++++++++ 4 files changed, 232 insertions(+), 91 deletions(-) create mode 100644 src/Plugins/MEDReader/Test/testMEDReader16.py diff --git a/src/Plugins/MEDReader/IO/vtkELNOFilter.cxx b/src/Plugins/MEDReader/IO/vtkELNOFilter.cxx index aee584bc..8af02927 100644 --- a/src/Plugins/MEDReader/IO/vtkELNOFilter.cxx +++ b/src/Plugins/MEDReader/IO/vtkELNOFilter.cxx @@ -27,6 +27,7 @@ #include "vtkFieldData.h" #include "vtkCellData.h" #include "vtkPointData.h" +#include "vtkCell.h" #include "vtkInformationQuadratureSchemeDefinitionVectorKey.h" #include "vtkQuadratureSchemeDefinition.h" #include "vtkUnstructuredGrid.h" @@ -45,99 +46,116 @@ vtkELNOFilter::~vtkELNOFilter() { } -int vtkELNOFilter::RequestData(vtkInformation *request, - vtkInformationVector **input, vtkInformationVector *output) +int vtkELNOFilter::RequestData(vtkInformation *request, vtkInformationVector **input, vtkInformationVector *output) { - vtkUnstructuredGrid *usgIn = vtkUnstructuredGrid::SafeDownCast( - input[0]->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT())); + vtkUnstructuredGrid *usgIn(vtkUnstructuredGrid::SafeDownCast( input[0]->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT()))); + vtkPolyData *pdOut(vtkPolyData::SafeDownCast(output->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT()))); - vtkPolyData *pdOut = vtkPolyData::SafeDownCast( - output->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT())); - - vtkDataArray* array = this->GetInputArrayToProcess(0, input); - vtkIdTypeArray* offsets = vtkIdTypeArray::SafeDownCast( - this->GetInputArrayToProcess(0, input)); + vtkDataArray *array(this->GetInputArrayToProcess(0, input)); + vtkIdTypeArray* offsets(vtkIdTypeArray::SafeDownCast(this->GetInputArrayToProcess(0, input))); if(usgIn == NULL || offsets == NULL || pdOut == NULL) { - vtkDebugMacro("vtkELNOFilter no correctly configured : offsets = " << offsets); - return 1; + vtkDebugMacro("vtkELNOFilter no correctly configured : offsets = " << offsets); + return 1; } - vtkInformation *info = offsets->GetInformation(); - vtkInformationQuadratureSchemeDefinitionVectorKey *key = - vtkQuadratureSchemeDefinition::DICTIONARY(); + vtkInformation *info(offsets->GetInformation()); + vtkInformationQuadratureSchemeDefinitionVectorKey *key(vtkQuadratureSchemeDefinition::DICTIONARY()); if(!key->Has(info)) { - vtkDebugMacro("Dictionary is not present in array " << offsets->GetName() << " " << offsets << " Aborting." ); - return 1; + vtkDebugMacro("Dictionary is not present in array " << offsets->GetName() << " " << offsets << " Aborting." ); + return 1; } - int res = this->Superclass::RequestData(request, input, output); + int res(this->Superclass::RequestData(request, input, output)); if(res == 0) - { return 0; - } - - int dictSize = key->Size(info); - vtkQuadratureSchemeDefinition **dict = - new vtkQuadratureSchemeDefinition *[dictSize]; + + int dictSize(key->Size(info)); + vtkQuadratureSchemeDefinition **dict(new vtkQuadratureSchemeDefinition *[dictSize]); key->GetRange(info, dict, 0, 0, dictSize); - vtkIdType ncell = usgIn->GetNumberOfCells(); - vtkPoints *points = pdOut->GetPoints(); - vtkIdType start = 0; + vtkIdType ncell(usgIn->GetNumberOfCells()); + vtkPoints *points(pdOut->GetPoints()); + vtkIdType start(0); for(vtkIdType cellId = 0; cellId < ncell; cellId++) { - vtkIdType offset = offsets->GetValue(cellId); - int cellType = usgIn->GetCellType(cellId); - // a simple check to see if a scheme really exists for this cell type. - // should not happen if the cell type has not been modified. - if(dict[cellType] == NULL) - continue; - int np = dict[cellType]->GetNumberOfQuadraturePoints(); - double center[3] = {0, 0, 0}; - for(int id = start; id < start + np; id++) - { - double *position = points->GetPoint(id); - center[0] += position[0]; - center[1] += position[1]; - center[2] += position[2]; - } - center[0] /= np; - center[1] /= np; - center[2] /= np; - for(int id = start; id < start + np; id++) - { - double *position = points->GetPoint(id); - double newpos[3]; - newpos[0] = position[0] * this->ShrinkFactor + center[0] * (1 - - this->ShrinkFactor); - newpos[1] = position[1] * this->ShrinkFactor + center[1] * (1 - - this->ShrinkFactor); - newpos[2] = position[2] * this->ShrinkFactor + center[2] * (1 - - this->ShrinkFactor); - points->SetPoint(id, newpos); - } - start += np; + vtkIdType offset(offsets->GetValue(cellId)); + int cellType(usgIn->GetCellType(cellId)); + // a simple check to see if a scheme really exists for this cell type. + // should not happen if the cell type has not been modified. + if(dict[cellType] == NULL) + continue; + int np = dict[cellType]->GetNumberOfQuadraturePoints(); + double center[3] = {0, 0, 0}; + for(int id = start; id < start + np; id++) + { + double *position = points->GetPoint(id); + center[0] += position[0]; + center[1] += position[1]; + center[2] += position[2]; + } + center[0] /= np; + center[1] /= np; + center[2] /= np; + for(int id = start; id < start + np; id++) + { + double *position = points->GetPoint(id); + double newpos[3]; + newpos[0] = position[0] * this->ShrinkFactor + center[0] * (1 - this->ShrinkFactor); + newpos[1] = position[1] * this->ShrinkFactor + center[1] * (1 - this->ShrinkFactor); + newpos[2] = position[2] * this->ShrinkFactor + center[2] * (1 - this->ShrinkFactor); + points->SetPoint(id, newpos); + } + start += np; } //// bug EDF 8407 PARAVIS - mantis 22610 vtkFieldData *fielddata(usgIn->GetFieldData()); - for(int index = 0; index < fielddata->GetNumberOfArrays(); index++) + for(int index=0;indexGetNumberOfArrays();index++) { vtkDataArray *data(fielddata->GetArray(index)); vtkInformation *info(data->GetInformation()); const char *arrayOffsetName(info->Get(vtkQuadratureSchemeDefinition::QUADRATURE_OFFSET_ARRAY_NAME())); + vtkIdTypeArray *offData(0); bool isELNO(false); if(arrayOffsetName) { vtkCellData *cellData(usgIn->GetCellData()); - vtkDataArray *offData(cellData->GetArray(arrayOffsetName)); - isELNO=offData->GetInformation()->Get(MEDUtilities::ELNO())==1; + vtkDataArray *offDataTmp(cellData->GetArray(arrayOffsetName)); + isELNO=offDataTmp->GetInformation()->Get(MEDUtilities::ELNO())==1; + offData=dynamic_cast(offDataTmp); } - if(isELNO) + if(isELNO && offData) { - pdOut->GetPointData()->AddArray(data); + vtkIdType nbCellsInput(usgIn->GetNumberOfCells()); + if(nbCellsInput==0) + continue ;//no cells -> no fields + // 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. + vtkCell *cell(usgIn->GetCell(nbCellsInput-1)); + bool statement0(offData->GetTuple1(0)==0); + bool statement1(offData->GetTuple1(nbCellsInput-1)+cell->GetNumberOfPoints()==data->GetNumberOfTuples()); + if(statement0 && statement1) + pdOut->GetPointData()->AddArray(data);//We are lucky ! No extraction needed. + else + {//not lucky ! Extract values from data. A previous threshold has been done... Bug EDF8662 + vtkDataArray *newArray(data->NewInstance()); + newArray->SetName(data->GetName()); + pdOut->GetPointData()->AddArray(newArray); + newArray->SetNumberOfComponents(data->GetNumberOfComponents()); + newArray->SetNumberOfTuples(pdOut->GetNumberOfPoints()); + newArray->CopyComponentNames(data); + newArray->Delete(); + vtkIdType *offsetPtr(offData->GetPointer(0)); + vtkIdType zeId(0); + for(vtkIdType cellId=0;cellIdGetCell(cellId)); + vtkIdType nbPoints(cell->GetNumberOfPoints()),offset(offsetPtr[cellId]); + for(vtkIdType j=0;jSetTuple(zeId,offsetPtr[cellId]+j,data); + } + } } } return 1; @@ -146,6 +164,5 @@ int vtkELNOFilter::RequestData(vtkInformation *request, void vtkELNOFilter::PrintSelf(ostream& os, vtkIndent indent) { this->Superclass::PrintSelf(os, indent); - os << indent << "ShrinkFactor : " << this->ShrinkFactor << endl; } diff --git a/src/Plugins/MEDReader/IO/vtkELNOMeshFilter.cxx b/src/Plugins/MEDReader/IO/vtkELNOMeshFilter.cxx index deffb070..ec410710 100644 --- a/src/Plugins/MEDReader/IO/vtkELNOMeshFilter.cxx +++ b/src/Plugins/MEDReader/IO/vtkELNOMeshFilter.cxx @@ -108,7 +108,7 @@ int vtkELNOMeshFilter::RequestData(vtkInformation *request, for(int index = 0; index < fielddata->GetNumberOfArrays(); index++) { vtkDataArray *data(fielddata->GetArray(index)); - vtkQuadratureSchemeDefinition **dict = 0; + vtkQuadratureSchemeDefinition **dict(0); vtkInformationQuadratureSchemeDefinitionVectorKey *key(vtkQuadratureSchemeDefinition::DICTIONARY()); if(key->Has(data->GetInformation())) { @@ -120,15 +120,17 @@ int vtkELNOMeshFilter::RequestData(vtkInformation *request, continue; vtkInformation *info(data->GetInformation()); - const char *arrayOffsetName = info->Get(vtkQuadratureSchemeDefinition::QUADRATURE_OFFSET_ARRAY_NAME()); - - bool isELGA(false); + const char *arrayOffsetName(info->Get(vtkQuadratureSchemeDefinition::QUADRATURE_OFFSET_ARRAY_NAME())); + vtkIdTypeArray *offData(0); + bool isELGA(false),isELNO(false); if(arrayOffsetName) { vtkFieldData *cellData(usgIn->GetCellData()); - vtkDataArray *offData(cellData->GetArray(arrayOffsetName)); - isELGA=offData->GetInformation()->Get(MEDUtilities::ELGA())==1; + vtkDataArray *offDataTmp(cellData->GetArray(arrayOffsetName)); + isELGA=offDataTmp->GetInformation()->Get(MEDUtilities::ELGA())==1; + isELNO=offDataTmp->GetInformation()->Get(MEDUtilities::ELNO())==1; + offData=dynamic_cast(offDataTmp); } if(arrayOffsetName == NULL || isELGA) @@ -141,39 +143,56 @@ int vtkELNOMeshFilter::RequestData(vtkInformation *request, } else { - vtkDataArray* newArray = data->NewInstance(); + vtkDataArray *newArray(data->NewInstance()); newArray->SetName(data->GetName()); usgOut->GetPointData()->AddArray(newArray); newArray->SetNumberOfComponents(data->GetNumberOfComponents()); newArray->SetNumberOfTuples(usgOut->GetNumberOfPoints()); newArray->CopyComponentNames(data); newArray->Delete(); - vtkIdList *ids(vtkIdList::New()); - vtkIdType offset(0); - for(vtkIdType cellId = 0; cellId < ncell; cellId++) + if(isELGA) { - int cellType = shrinked->GetCellType(cellId); - shrinked->GetCellPoints(cellId, ids); - for(int id = 0; id < dict[cellType]->GetNumberOfQuadraturePoints(); id++) + vtkIdList *ids(vtkIdList::New()); + vtkIdType offset(0); + for(vtkIdType cellId=0;cellIdGetShapeFunctionWeights(id); - int j; - for(j = 0; j < dict[cellType]->GetNumberOfNodes(); j++) - { - if(w[j] == 1.0) - break; - } - if(j == dict[cellType]->GetNumberOfNodes()) + int cellType(shrinked->GetCellType(cellId)); + shrinked->GetCellPoints(cellId, ids); + for(int id = 0; id < dict[cellType]->GetNumberOfQuadraturePoints(); id++) { - j = id; + const double * w = dict[cellType]->GetShapeFunctionWeights(id); + int j; + for(j = 0; j < dict[cellType]->GetNumberOfNodes(); j++) + { + if(w[j] == 1.0) + break; + } + if(j == dict[cellType]->GetNumberOfNodes()) + { + j = id; + } + newArray->SetTuple(ids->GetId(id), offset + j, data); } - newArray->SetTuple(ids->GetId(id), offset + j, data); + vtkCell *cell(usgIn->GetCell(cellId)); + vtkIdType nbptsInCell(cell->GetNumberOfPoints()); + offset+=nbptsInCell; + } + ids->FastDelete(); + } + else if(offData && isELNO) + { + vtkIdType *offsetPtr(offData->GetPointer(0)); + vtkIdType zeId(0); + for(vtkIdType cellId=0;cellIdGetCell(cellId)); + vtkIdType nbPoints(cell->GetNumberOfPoints()),offset(offsetPtr[cellId]); + for(vtkIdType j=0;jSetTuple(zeId,offsetPtr[cellId]+j,data); } - vtkCell *cell(usgIn->GetCell(cellId)); - vtkIdType nbptsInCell(cell->GetNumberOfPoints()); - offset+=nbptsInCell; } - ids->FastDelete(); + else + continue ; } delete [] dict; } diff --git a/src/Plugins/MEDReader/Test/CMakeLists.txt b/src/Plugins/MEDReader/Test/CMakeLists.txt index bac570f6..42b017da 100644 --- a/src/Plugins/MEDReader/Test/CMakeLists.txt +++ b/src/Plugins/MEDReader/Test/CMakeLists.txt @@ -50,3 +50,5 @@ ADD_TEST(testMEDReader14 ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/testME SET_TESTS_PROPERTIES(testMEDReader14 PROPERTIES ENVIRONMENT "${tests_env}") ADD_TEST(testMEDReader15 ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/testMEDReader15.py) SET_TESTS_PROPERTIES(testMEDReader15 PROPERTIES ENVIRONMENT "${tests_env}") +ADD_TEST(testMEDReader16 ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/testMEDReader16.py) +SET_TESTS_PROPERTIES(testMEDReader16 PROPERTIES ENVIRONMENT "${tests_env}") diff --git a/src/Plugins/MEDReader/Test/testMEDReader16.py b/src/Plugins/MEDReader/Test/testMEDReader16.py new file mode 100644 index 00000000..012c5110 --- /dev/null +++ b/src/Plugins/MEDReader/Test/testMEDReader16.py @@ -0,0 +1,103 @@ +# -*- coding: iso-8859-1 -*- +# Copyright (C) 2007-2014 CEA/DEN, EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Anthony Geay (EDF R&D) + +from MEDLoader import * + +""" This test is a non regression test of EDF8662 : This bug revealed that ELNOMesh and ELNOPoints do not behave correctly after the call of ExtractGroup""" + +fname="testMEDReader16.med" +outImgName="testMEDReader16.png" + +arr=DataArrayDouble([0,1,2]) +m=MEDCouplingCMesh() ; m.setCoords(arr,arr) ; m=m.buildUnstructured() ; m.setName("Mesh") +mm=MEDFileUMesh() ; mm.setMeshAtLevel(0,m) +grp0=DataArrayInt([0,1]) ; grp0.setName("grp0") +grp1=DataArrayInt([2,3]) ; grp1.setName("grp1") +grp2=DataArrayInt([1,2]) ; grp2.setName("grp2") +grp3=DataArrayInt([0,1,2,3]) ; grp3.setName("grp3") +mm.setGroupsAtLevel(0,[grp0,grp1,grp2,grp3]) +f=MEDCouplingFieldDouble(ON_GAUSS_NE) ; f.setMesh(m) ; f.setName("MyField") ; f.setTime(0.,0,0) +arr2=DataArrayDouble(4*4*2) ; arr2.iota() ; arr2.rearrange(2) ; arr2.setInfoOnComponents(["aa","bbb"]) +f.setArray(arr2) ; arr2+=0.1 ; f.checkCoherency() +mm.write(fname,2) +MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fname,f) +# +from paraview.simple import * +from paraview import servermanager + +paraview.simple._DisableFirstRenderCameraReset() +reader=MEDReader(FileName=fname) +ExpectedEntries=['TS0/Mesh/ComSup0/MyField@@][@@GSSNE','TS1/Mesh/ComSup0/Mesh@@][@@P0'] +assert(reader.GetProperty("FieldsTreeInfo")[::2]==ExpectedEntries) +reader.AllArrays=['TS0/Mesh/ComSup0/MyField@@][@@GSSNE'] +ExtractGroup1 = ExtractGroup(Input=reader) +ExtractGroup1.UpdatePipelineInformation() +ExtractGroup1.AllGroups=["GRP_grp1"] +ELNOMesh1=ELNOMesh(Input=ExtractGroup1) +ELNOPoints1=ELNOPoints(Input=ExtractGroup1) +ELNOPoints1.SelectSourceArray=['ELNO@MyField'] +for elt in [ELNOMesh1,ELNOPoints1]: + elnoMesh=servermanager.Fetch(ELNOPoints1,0) + vtkArrToTest=elnoMesh.GetBlock(0).GetPointData().GetArray("MyField") + assert(vtkArrToTest.GetNumberOfTuples()==8) + assert(vtkArrToTest.GetNumberOfComponents()==2) + assert(vtkArrToTest.GetComponentName(0)==arr2.getInfoOnComponent(0)) + assert(vtkArrToTest.GetComponentName(1)==arr2.getInfoOnComponent(1)) + vals=[vtkArrToTest.GetValue(i) for i in xrange(16)] + assert(arr2[8:].isEqualWithoutConsideringStr(DataArrayDouble(vals,8,2),1e-12)) + pass +# +ExtractGroup1.AllGroups=["GRP_grp2"] +for elt in [ELNOMesh1,ELNOPoints1]: + elnoMesh=servermanager.Fetch(ELNOMesh1) + vtkArrToTest=elnoMesh.GetBlock(0).GetPointData().GetArray("MyField") + assert(vtkArrToTest.GetNumberOfTuples()==8) + assert(vtkArrToTest.GetNumberOfComponents()==2) + assert(vtkArrToTest.GetComponentName(0)==arr2.getInfoOnComponent(0)) + assert(vtkArrToTest.GetComponentName(1)==arr2.getInfoOnComponent(1)) + vals=[vtkArrToTest.GetValue(i) for i in xrange(16)] + assert(arr2[4:12].isEqualWithoutConsideringStr(DataArrayDouble(vals,8,2),1e-12)) + pass +# important to check that if all the field is present that it is OK (check of the optimization) +ExtractGroup1.AllGroups=["GRP_grp3"] +for elt in [ELNOMesh1,ELNOPoints1]: + elnoMesh=servermanager.Fetch(ELNOMesh1) + vtkArrToTest=elnoMesh.GetBlock(0).GetPointData().GetArray("MyField") + assert(vtkArrToTest.GetNumberOfTuples()==16) + assert(vtkArrToTest.GetNumberOfComponents()==2) + assert(vtkArrToTest.GetComponentName(0)==arr2.getInfoOnComponent(0)) + assert(vtkArrToTest.GetComponentName(1)==arr2.getInfoOnComponent(1)) + vals=[vtkArrToTest.GetValue(i) for i in xrange(32)] + assert(arr2.isEqualWithoutConsideringStr(DataArrayDouble(vals,16,2),1e-12)) + pass +ELNOMesh1=ELNOMesh(Input=reader) +ELNOPoints1=ELNOPoints(Input=reader) +ELNOPoints1.SelectSourceArray=['ELNO@MyField'] +for elt in [ELNOMesh1,ELNOPoints1]: + elnoMesh=servermanager.Fetch(ELNOMesh1) + vtkArrToTest=elnoMesh.GetBlock(0).GetPointData().GetArray("MyField") + assert(vtkArrToTest.GetNumberOfTuples()==16) + assert(vtkArrToTest.GetNumberOfComponents()==2) + assert(vtkArrToTest.GetComponentName(0)==arr2.getInfoOnComponent(0)) + assert(vtkArrToTest.GetComponentName(1)==arr2.getInfoOnComponent(1)) + vals=[vtkArrToTest.GetValue(i) for i in xrange(32)] + assert(arr2.isEqualWithoutConsideringStr(DataArrayDouble(vals,16,2),1e-12)) + pass -- 2.39.2