#include "vtkFieldData.h"
#include "vtkCellData.h"
#include "vtkPointData.h"
+#include "vtkCell.h"
#include "vtkInformationQuadratureSchemeDefinitionVectorKey.h"
#include "vtkQuadratureSchemeDefinition.h"
#include "vtkUnstructuredGrid.h"
{
}
-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;index<fielddata->GetNumberOfArrays();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<vtkIdTypeArray *>(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;cellId<nbCellsInput;cellId++)
+ {
+ vtkCell *cell(usgIn->GetCell(cellId));
+ vtkIdType nbPoints(cell->GetNumberOfPoints()),offset(offsetPtr[cellId]);
+ for(vtkIdType j=0;j<nbPoints;j++,zeId++)
+ newArray->SetTuple(zeId,offsetPtr[cellId]+j,data);
+ }
+ }
}
}
return 1;
void vtkELNOFilter::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os, indent);
-
os << indent << "ShrinkFactor : " << this->ShrinkFactor << endl;
}
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()))
{
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<vtkIdTypeArray *>(offDataTmp);
}
if(arrayOffsetName == NULL || isELGA)
}
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;cellId<ncell;cellId++)
{
- 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())
+ 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;cellId<ncell;cellId++)
+ {
+ vtkCell *cell(shrinked->GetCell(cellId));
+ vtkIdType nbPoints(cell->GetNumberOfPoints()),offset(offsetPtr[cellId]);
+ for(vtkIdType j=0;j<nbPoints;j++,zeId++)
+ newArray->SetTuple(zeId,offsetPtr[cellId]+j,data);
}
- vtkCell *cell(usgIn->GetCell(cellId));
- vtkIdType nbptsInCell(cell->GetNumberOfPoints());
- offset+=nbptsInCell;
}
- ids->FastDelete();
+ else
+ continue ;
}
delete [] dict;
}
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}")
--- /dev/null
+# -*- 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
