From fe342008f22de4504754b8503f339445f8d1d759 Mon Sep 17 00:00:00 2001 From: ageay Date: Tue, 27 Jul 2010 05:53:21 +0000 Subject: [PATCH] Addition of utilities for profile RW in MEDLoader. --- src/MEDCoupling/MEDCouplingMemArray.cxx | 56 ++++++ src/MEDCoupling/MEDCouplingMemArray.hxx | 2 + src/MEDCoupling/MEDCouplingUMesh.cxx | 63 +++++- src/MEDCoupling/MEDCouplingUMesh.hxx | 2 + .../Test/MEDCouplingBasicsTest.hxx | 2 + .../Test/MEDCouplingBasicsTest1.cxx | 43 ++++ src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 42 ++++ src/MEDCoupling_Swig/MEDCouplingTypemaps.i | 28 +++ src/MEDCoupling_Swig/libMEDCoupling_Swig.i | 12 ++ src/MEDLoader/MEDLoader.cxx | 190 +++++++++++++++++- src/MEDLoader/MEDLoader.hxx | 13 +- src/MEDLoader/Test/MEDLoaderTest.cxx | 9 + 12 files changed, 447 insertions(+), 15 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index 430ca9bee..84ade1a95 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -101,6 +101,34 @@ DataArrayInt *DataArrayDouble::convertToIntArr() const return ret; } +/*! + * This methods has a similar behaviour than std::string::substr. This method returns a newly created DataArrayInt that is part of this with same number of components. + * The intervall is specified by [tupleIdBg,tupleIdEnd) except if tupleIdEnd ==-1 in this case the [tupleIdBg,this->end()) will be kept. + * This method check that interval is valid regarding this, if not an exception will be thrown. + */ +DataArrayDouble *DataArrayDouble::substr(int tupleIdBg, int tupleIdEnd) const throw(INTERP_KERNEL::Exception) +{ + int nbt=getNumberOfTuples(); + if(tupleIdBg<0) + throw INTERP_KERNEL::Exception("DataArrayInt::substr : The tupleIdBg parameter must be greater than 0 !"); + if(tupleIdBg>=nbt) + throw INTERP_KERNEL::Exception("DataArrayInt::substr : The tupleIdBg parameter is greater or equal than number of tuples !"); + int trueEnd=tupleIdEnd; + if(tupleIdEnd!=-1) + { + if(tupleIdEnd>nbt) + throw INTERP_KERNEL::Exception("DataArrayInt::substr : The tupleIdBg parameter is greater or equal than number of tuples !"); + } + else + trueEnd=nbt; + int nbComp=getNumberOfComponents(); + DataArrayDouble *ret=DataArrayDouble::New(); + ret->alloc(trueEnd-tupleIdBg,nbComp); + ret->copyStringInfoFrom(*this); + std::copy(getConstPointer()+tupleIdBg*nbComp,getConstPointer()+trueEnd*nbComp,ret->getPointer()); + return ret; +} + void DataArrayDouble::setArrayIn(DataArrayDouble *newArray, DataArrayDouble* &arrayToSet) { if(newArray!=arrayToSet) @@ -359,6 +387,34 @@ DataArrayDouble *DataArrayInt::convertToDblArr() const return ret; } +/*! + * This methods has a similar behaviour than std::string::substr. This method returns a newly created DataArrayInt that is part of this with same number of components. + * The intervall is specified by [tupleIdBg,tupleIdEnd) except if tupleIdEnd ==-1 in this case the [tupleIdBg,this->end()) will be kept. + * This method check that interval is valid regarding this, if not an exception will be thrown. + */ +DataArrayInt *DataArrayInt::substr(int tupleIdBg, int tupleIdEnd) const throw(INTERP_KERNEL::Exception) +{ + int nbt=getNumberOfTuples(); + if(tupleIdBg<0) + throw INTERP_KERNEL::Exception("DataArrayInt::substr : The tupleIdBg parameter must be greater than 0 !"); + if(tupleIdBg>=nbt) + throw INTERP_KERNEL::Exception("DataArrayInt::substr : The tupleIdBg parameter is greater or equal than number of tuples !"); + int trueEnd=tupleIdEnd; + if(tupleIdEnd!=-1) + { + if(tupleIdEnd>nbt) + throw INTERP_KERNEL::Exception("DataArrayInt::substr : The tupleIdBg parameter is greater or equal than number of tuples !"); + } + else + trueEnd=nbt; + int nbComp=getNumberOfComponents(); + DataArrayInt *ret=DataArrayInt::New(); + ret->alloc(trueEnd-tupleIdBg,nbComp); + ret->copyStringInfoFrom(*this); + std::copy(getConstPointer()+tupleIdBg*nbComp,getConstPointer()+trueEnd*nbComp,ret->getPointer()); + return ret; +} + void DataArrayInt::reAlloc(int nbOfTuples) { _mem.reAlloc(_info_on_compo.size()*nbOfTuples); diff --git a/src/MEDCoupling/MEDCouplingMemArray.hxx b/src/MEDCoupling/MEDCouplingMemArray.hxx index 9974aa563..eec90bf5e 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.hxx +++ b/src/MEDCoupling/MEDCouplingMemArray.hxx @@ -113,6 +113,7 @@ namespace ParaMEDMEM //!alloc or useArray should have been called before. MEDCOUPLING_EXPORT void reAlloc(int nbOfTuples); MEDCOUPLING_EXPORT DataArrayInt *convertToIntArr() const; + MEDCOUPLING_EXPORT DataArrayDouble *substr(int tupleIdBg, int tupleIdEnd=-1) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void getTuple(int tupleId, double *res) const { std::copy(_mem.getConstPointerLoc(tupleId*_info_on_compo.size()),_mem.getConstPointerLoc((tupleId+1)*_info_on_compo.size()),res); } MEDCOUPLING_EXPORT double getIJ(int tupleId, int compoId) const { return _mem[tupleId*_info_on_compo.size()+compoId]; } MEDCOUPLING_EXPORT void setIJ(int tupleId, int compoId, double newVal) { _mem[tupleId*_info_on_compo.size()+compoId]=newVal; } @@ -150,6 +151,7 @@ namespace ParaMEDMEM //!alloc or useArray should have been called before. MEDCOUPLING_EXPORT void reAlloc(int nbOfTuples); MEDCOUPLING_EXPORT DataArrayDouble *convertToDblArr() const; + MEDCOUPLING_EXPORT DataArrayInt *substr(int tupleIdBg, int tupleIdEnd=-1) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void getTuple(int tupleId, int *res) const { std::copy(_mem.getConstPointerLoc(tupleId*_info_on_compo.size()),_mem.getConstPointerLoc((tupleId+1)*_info_on_compo.size()),res); } MEDCOUPLING_EXPORT int getIJ(int tupleId, int compoId) const { return _mem[tupleId*_info_on_compo.size()+compoId]; } MEDCOUPLING_EXPORT void setIJ(int tupleId, int compoId, int newVal) { _mem[tupleId*_info_on_compo.size()+compoId]=newVal; } diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 342fc33a8..26f538942 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -1837,7 +1837,7 @@ DataArrayInt *MEDCouplingUMesh::rearrange2ConsecutiveCellTypes() /*! * This methods split this into as mush as untructured meshes that consecutive set of same type cells. * So this method has typically a sense if MEDCouplingUMesh::checkConsecutiveCellTypes has a sense. - * This method makes asumption (no check) that connectivity is correctly set before calling. + * This method makes asumption that connectivity is correctly set before calling. */ std::vector MEDCouplingUMesh::splitByType() const { @@ -1863,6 +1863,67 @@ std::vector MEDCouplingUMesh::splitByType() const return ret; } +/*! + * This method makes the assumption that da->getNumberOfTuples()getNumberOfCells(). This method makes the assumption that ids contained in 'da' + * are in [0:getNumberOfCells()) + */ +DataArrayInt *MEDCouplingUMesh::convertCellArrayPerGeoType(const DataArrayInt *da) const throw(INTERP_KERNEL::Exception) +{ + checkFullyDefined(); + const int *conn=_nodal_connec->getConstPointer(); + const int *connI=_nodal_connec_index->getConstPointer(); + int nbOfCells=getNumberOfCells(); + std::set types=getAllTypes(); + int *tmp=new int[nbOfCells]; + for(std::set::const_iterator iter=types.begin();iter!=types.end();iter++) + { + int j=0; + for(const int *i=connI;i!=connI+nbOfCells;i++) + if((INTERP_KERNEL::NormalizedCellType)conn[*i]==(*iter)) + tmp[std::distance(connI,i)]=j++; + } + DataArrayInt *ret=DataArrayInt::New(); + ret->alloc(da->getNumberOfTuples(),da->getNumberOfComponents()); + ret->copyStringInfoFrom(*da); + int *retPtr=ret->getPointer(); + const int *daPtr=da->getConstPointer(); + int nbOfElems=da->getNbOfElems(); + for(int k=0;k& idsPerGeoType) +{ + std::vector idsTokeep; + int nbOfCells=getNumberOfCells(); + int j=0; + for(int i=0;i(ret); + if(!ret2) + { + ret->decrRef(); + return 0; + } + ret2->setName(getName()); + return ret2; +} + /*! * Returns a newly created mesh (with ref count ==1) that contains merge of 'this' and 'other'. */ diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 3dde3bde2..f3fdf9a77 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -96,6 +96,8 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT bool checkConsecutiveCellTypes() const; MEDCOUPLING_EXPORT DataArrayInt *rearrange2ConsecutiveCellTypes(); MEDCOUPLING_EXPORT std::vector splitByType() const; + MEDCOUPLING_EXPORT DataArrayInt *convertCellArrayPerGeoType(const DataArrayInt *da) const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT MEDCouplingUMesh *keepSpecifiedCells(INTERP_KERNEL::NormalizedCellType type, const std::vector& idsPerGeoType); // MEDCOUPLING_EXPORT MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const; MEDCOUPLING_EXPORT DataArrayDouble *getBarycenterAndOwner() const; diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx index 7307e3834..503eba3c0 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx @@ -34,6 +34,7 @@ namespace ParaMEDMEM CPPUNIT_TEST_SUITE(MEDCouplingBasicsTest); CPPUNIT_TEST( testArray ); CPPUNIT_TEST( testArray2 ); + CPPUNIT_TEST( testArray3 ); CPPUNIT_TEST( testMesh ); CPPUNIT_TEST( testMeshPointsCloud ); CPPUNIT_TEST( testMeshM1D ); @@ -144,6 +145,7 @@ namespace ParaMEDMEM public: void testArray(); void testArray2(); + void testArray3(); void testMesh(); void testMeshPointsCloud(); void testMeshM1D(); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx index 46e79f9ec..3ce905d5e 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx @@ -78,6 +78,49 @@ void MEDCouplingBasicsTest::testArray2() arr->decrRef(); } +void MEDCouplingBasicsTest::testArray3() +{ + DataArrayInt *arr1=DataArrayInt::New(); + arr1->alloc(7,2); + int *tmp=arr1->getPointer(); + const int arr1Ref[14]={0,10,1,11,2,12,3,13,4,14,5,15,6,16}; + std::copy(arr1Ref,arr1Ref+14,tmp); + CPPUNIT_ASSERT_EQUAL(7,arr1->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(2,arr1->getNumberOfComponents()); + CPPUNIT_ASSERT(std::equal(arr1Ref,arr1Ref+14,arr1->getConstPointer())); + DataArrayInt *arr2=arr1->substr(3); + CPPUNIT_ASSERT_EQUAL(4,arr2->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(2,arr2->getNumberOfComponents()); + CPPUNIT_ASSERT(std::equal(arr1Ref+6,arr1Ref+14,arr2->getConstPointer())); + arr2->decrRef(); + DataArrayInt *arr3=arr1->substr(2,5); + CPPUNIT_ASSERT_EQUAL(3,arr3->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(2,arr3->getNumberOfComponents()); + CPPUNIT_ASSERT(std::equal(arr1Ref+4,arr1Ref+10,arr3->getConstPointer())); + arr1->decrRef(); + arr3->decrRef(); + // + DataArrayDouble *arr4=DataArrayDouble::New(); + arr4->alloc(7,2); + double *tmp2=arr4->getPointer(); + const int arr4Ref[14]={0.8,10.8,1.9,11.9,2.1,12.1,3.2,13.2,4.3,14.3,5.4,15.4,6.5,16.5}; + std::copy(arr4Ref,arr4Ref+14,tmp2); + CPPUNIT_ASSERT_EQUAL(7,arr4->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(2,arr4->getNumberOfComponents()); + CPPUNIT_ASSERT(std::equal(arr4Ref,arr4Ref+14,arr4->getConstPointer())); + DataArrayDouble *arr5=arr4->substr(3); + CPPUNIT_ASSERT_EQUAL(4,arr5->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(2,arr5->getNumberOfComponents()); + CPPUNIT_ASSERT(std::equal(arr4Ref+6,arr4Ref+14,arr5->getConstPointer())); + arr5->decrRef(); + DataArrayDouble *arr6=arr4->substr(2,5); + CPPUNIT_ASSERT_EQUAL(3,arr6->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(2,arr6->getNumberOfComponents()); + CPPUNIT_ASSERT(std::equal(arr4Ref+4,arr4Ref+10,arr6->getConstPointer())); + arr4->decrRef(); + arr6->decrRef(); +} + void MEDCouplingBasicsTest::testMesh() { const int nbOfCells=6; diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index f6254d4fd..8fef761ba 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -35,6 +35,48 @@ class MEDCouplingBasicsTest(unittest.TestCase): arr3=arr2.convertToDblArr(); self.failUnless(arr.isEqual(arr3,1e-14)) pass + + def testArray3(self): + arr1=DataArrayInt.New(); + arr1Ref=[0,10,1,11,2,12,3,13,4,14,5,15,6,16] + arr1.setValues(arr1Ref,7,2); + self.failUnlessEqual(7,arr1.getNumberOfTuples()); + self.failUnlessEqual(2,arr1.getNumberOfComponents()); + self.failUnlessEqual(arr1Ref,arr1.getValues()); + arr2=arr1.substr(3); + self.failUnlessEqual(4,arr2.getNumberOfTuples()); + self.failUnlessEqual(2,arr2.getNumberOfComponents()); + self.failUnlessEqual(arr1Ref[6:],arr2.getValues()); + arr3=arr1.substr(2,5); + self.failUnlessEqual(3,arr3.getNumberOfTuples()); + self.failUnlessEqual(2,arr3.getNumberOfComponents()); + self.failUnlessEqual(arr1Ref[4:10],arr3.getValues()); + # + arr4=DataArrayDouble.New(); + arr4Ref=[0.8,10.8,1.9,11.9,2.1,12.1,3.2,13.2,4.3,14.3,5.4,15.4,6.5,16.5] + arr4.setValues(arr4Ref,7,2); + self.failUnlessEqual(7,arr4.getNumberOfTuples()); + self.failUnlessEqual(2,arr4.getNumberOfComponents()); + tmp=arr4.getValues() + for i in xrange(14): + self.failUnless(abs(arr4Ref[i]-tmp[i])<1e-14); + pass + arr5=arr4.substr(3); + self.failUnlessEqual(4,arr5.getNumberOfTuples()); + self.failUnlessEqual(2,arr5.getNumberOfComponents()); + tmp=arr5.getValues() + for i in xrange(8): + self.failUnless(abs(arr4Ref[6+i]-tmp[i])<1e-14); + pass + arr6=arr4.substr(2,5); + self.failUnlessEqual(3,arr6.getNumberOfTuples()); + self.failUnlessEqual(2,arr6.getNumberOfComponents()); + tmp=arr6.getValues() + for i in xrange(6): + self.failUnless(abs(arr4Ref[4+i]-tmp[i])<1e-14); + pass + pass + def testMesh(self): tab4=[1, 2, 8, 7, 2, 3, 9, 8, 3, 4, 10, 9, 4, 5, 11, 10, 5, diff --git a/src/MEDCoupling_Swig/MEDCouplingTypemaps.i b/src/MEDCoupling_Swig/MEDCouplingTypemaps.i index 6583ae968..5bc13876d 100644 --- a/src/MEDCoupling_Swig/MEDCouplingTypemaps.i +++ b/src/MEDCoupling_Swig/MEDCouplingTypemaps.i @@ -76,6 +76,34 @@ static int *convertPyToNewIntArr2(PyObject *pyLi, int *size) } } +static void convertPyToNewIntArr3(PyObject *pyLi, std::vector& arr) +{ + if(PyList_Check(pyLi)) + { + int size=PyList_Size(pyLi); + arr.resize(size); + for(int i=0;idecrRef();" @@ -355,6 +358,7 @@ namespace ParaMEDMEM //tools bool checkConsecutiveCellTypes() const; DataArrayInt *rearrange2ConsecutiveCellTypes(); + DataArrayInt *convertCellArrayPerGeoType(const DataArrayInt *da) const throw(INTERP_KERNEL::Exception); DataArrayInt *zipConnectivityTraducer(int compType); void getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const; MEDCouplingUMesh *buildDescendingConnectivity(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const; @@ -416,6 +420,14 @@ namespace ParaMEDMEM return ret; } + PyObject *keepSpecifiedCells(INTERP_KERNEL::NormalizedCellType type, PyObject *ids) + { + std::vector idsPerGeoType; + convertPyToNewIntArr3(ids,idsPerGeoType); + MEDCouplingUMesh *ret=self->keepSpecifiedCells(type,idsPerGeoType); + return SWIG_NewPointerObj(SWIG_as_voidptr(ret),SWIGTYPE_p_ParaMEDMEM__MEDCouplingUMesh, SWIG_POINTER_OWN | 0 ); + } + PyObject *getCellsContainingPoints(PyObject *p, int nbOfPoints, double eps) const { int sz; diff --git a/src/MEDLoader/MEDLoader.cxx b/src/MEDLoader/MEDLoader.cxx index fa0724c10..be1aa373f 100644 --- a/src/MEDLoader/MEDLoader.cxx +++ b/src/MEDLoader/MEDLoader.cxx @@ -109,6 +109,10 @@ med_geometrie_element typmai3[32] = { MED_POINT1,//0 MED_POLYEDRE//31 }; +double MEDLoader::_EPS_FOR_NODE_COMP=1.e-12; + +int MEDLoader::_COMP_FOR_CELL=0; + using namespace ParaMEDMEM; namespace MEDLoaderNS @@ -161,8 +165,11 @@ namespace MEDLoaderNS INTERP_KERNEL::NormalizedCellType type, std::vector& conn4MEDFile, std::vector& fam4MEDFile); ParaMEDMEM::MEDCouplingFieldDouble *readFieldDoubleLev1(const char *fileName, const char *meshName, int meshDimRelToMax, const char *fieldName, int iteration, int order, ParaMEDMEM::TypeOfField typeOfOutField); + med_idt appendFieldSimpleAtt(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, med_int& numdt, med_int& numo, med_float& dt); void appendFieldDirectly(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f); - void prepareCellFieldDoubleForWriting(const ParaMEDMEM::MEDCouplingFieldDouble *f, std::list& split); + void appendNodeProfileField(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, const int *thisMeshNodeIds); + void appendCellProfileField(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, const int *thisMeshCellIds); + void prepareCellFieldDoubleForWriting(const ParaMEDMEM::MEDCouplingFieldDouble *f, const int *cellIds, std::list& split); void writeUMeshDirectly(const char *fileName, ParaMEDMEM::MEDCouplingUMesh *mesh, const DataArrayInt *families, bool forceFromScratch); void writeUMeshesDirectly(const char *fileName, const char *meshName, const std::vector& meshes, bool forceFromScratch); void writeFieldAndMeshDirectly(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, bool forceFromScratch); @@ -171,6 +178,22 @@ namespace MEDLoaderNS const char WHITE_SPACES[]=" \n"; +/*! + * This method set the epsilon value used for node comparison when trying to buid a profile for a field on node/cell on an already written mesh. + */ +void MEDLoader::setEpsilonForNodeComp(double val) +{ + _EPS_FOR_NODE_COMP=val; +} + +/*! + * This method set the policy comparison when trying to fit the already written mesh on a field. The semantic of the policy is specified in MEDCouplingUMesh::zipConnectivityTraducer. + */ +void MEDLoader::setCompPolicyForCell(int val) +{ + _COMP_FOR_CELL=val; +} + /*! * @param lgth is the size of fam tab. For classical types conn is size of 'lgth'*number_of_nodes_in_type. * @param index is optionnal only for polys. Set it to 0 if it is not the case. @@ -201,8 +224,10 @@ void MEDLoader::MEDConnOfOneElemType::releaseArray() delete [] _global; } -MEDLoader::MEDFieldDoublePerCellType::MEDFieldDoublePerCellType(INTERP_KERNEL::NormalizedCellType type, double *values, int ncomp, int ntuple):_ntuple(ntuple),_ncomp(ncomp),_values(values),_type(type) +MEDLoader::MEDFieldDoublePerCellType::MEDFieldDoublePerCellType(INTERP_KERNEL::NormalizedCellType type, double *values, int ncomp, int ntuple, const int *cellIdPerType):_ntuple(ntuple),_ncomp(ncomp),_values(values),_type(type) { + if(cellIdPerType) + _cell_id_per_type.insert(_cell_id_per_type.end(),cellIdPerType,cellIdPerType+ntuple); } void MEDLoader::MEDFieldDoublePerCellType::releaseArray() @@ -571,7 +596,14 @@ void MEDLoaderNS::readFieldDoubleDataInMedFile(const char *fileName, const char } MEDchampLire(fid,(char *)meshName,(char *)fieldName,(unsigned char*)valr,MED_FULL_INTERLACE,MED_ALL,locname, pflname,MED_COMPACT,tabEnt[typeOfOutField],tabType[typeOfOutField][j],iteration,order); - field.push_back(MEDLoader::MEDFieldDoublePerCellType(typmai2[j],valr,ncomp,nval)); + int *pfl=0; + if(pflname[0]!='\0') + { + pfl=new int[nval]; + MEDprofilLire(fid,pfl,pflname); + } + field.push_back(MEDLoader::MEDFieldDoublePerCellType(typmai2[j],valr,ncomp,nval,pfl)); + delete [] pfl; delete [] dt_unit; delete [] maa_ass; } @@ -1208,6 +1240,20 @@ ParaMEDMEM::MEDCouplingFieldDouble *MEDLoaderNS::readFieldDoubleLev1(const char ParaMEDMEM::MEDCouplingUMesh *mesh=readUMeshFromFileLev1(fileName,meshName,meshDimRelToMax,familiesToKeep,typesToKeep,meshDim); if(typeOfOutField==ON_CELLS) MEDLoaderNS::keepSpecifiedMeshDim(fieldPerCellType,meshDim); + //for profiles + for(std::list::const_iterator iter=fieldPerCellType.begin();iter!=fieldPerCellType.end();iter++) + { + const std::vector& cellIds=(*iter).getCellIdPerType(); + if(!cellIds.empty()) + { + std::vector ci(cellIds.size()); + std::transform(cellIds.begin(),cellIds.end(),ci.begin(),std::bind2nd(std::plus(),-1)); + ParaMEDMEM::MEDCouplingUMesh *mesh2=mesh->keepSpecifiedCells((*iter).getType(),ci); + mesh->decrRef(); + mesh=mesh2; + } + } + // ParaMEDMEM::MEDCouplingFieldDouble *ret=ParaMEDMEM::MEDCouplingFieldDouble::New(typeOfOutField,ONE_TIME); ret->setName(fieldName); ret->setTime(time,iteration,order); @@ -1409,15 +1455,68 @@ void MEDLoaderNS::writeUMeshesDirectly(const char *fileName, const char *meshNam m->decrRef(); } -void MEDLoaderNS::appendFieldDirectly(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f) +/*! + * This method makes the assumption that f->getMesh() nodes are fully included in already written mesh in 'fileName'. + * @param thisMeshNodeIds points to a tab of size f->getMesh()->getNumberOfNodes() that says for a node i in f->getMesh() that its id is thisMeshNodeIds[i] is already written mesh. + */ +void MEDLoaderNS::appendNodeProfileField(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, const int *thisMeshNodeIds) +{ + //not implemented yet. + med_int numdt,numo; + med_float dt; + int nbComp=f->getNumberOfComponents(); + med_idt fid=appendFieldSimpleAtt(fileName,f,numdt,numo,dt); + MEDfermer(fid); +} + +/*! + * This method makes the assumption that f->getMesh() cells are fully included in already written mesh in 'fileName'. + * @param thisMeshCellIdsPerType points to a tab of size f->getMesh()->getNumberOfCells() that says for a cell i in f->getMesh() that its id is thisMeshCellIds[i] of corresponding type is already written mesh. + */ +void MEDLoaderNS::appendCellProfileField(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, const int *thisMeshCellIdsPerType) +{ + //not implemented yet. + med_int numdt,numo; + med_float dt; + int nbComp=f->getNumberOfComponents(); + med_idt fid=appendFieldSimpleAtt(fileName,f,numdt,numo,dt); + std::list split; + prepareCellFieldDoubleForWriting(f,thisMeshCellIdsPerType,split); + const double *pt=f->getArray()->getConstPointer(); + int number=0; + for(std::list::const_iterator iter=split.begin();iter!=split.end();iter++) + { + char *nommaa=MEDLoaderBase::buildEmptyString(MED_TAILLE_NOM); + strcpy(nommaa,f->getMesh()->getName()); + char *profileName=MEDLoaderBase::buildEmptyString(MED_TAILLE_NOM); + std::ostringstream oss; oss << "Pfl" << f->getName() << "_" << number++; + strcpy(profileName,oss.str().c_str()); + const std::vector& ids=(*iter).getCellIdPerType(); + int *profile=new int [ids.size()]; + std::transform(ids.begin(),ids.end(),profile,std::bind2nd(std::plus(),1)); + MEDprofilEcr(fid,profile,ids.size(),profileName); + delete [] profile; + MEDchampEcr(fid,nommaa,(char *)f->getName(),(unsigned char*)pt,MED_FULL_INTERLACE,(*iter).getNbOfTuple(), + (char *)MED_NOGAUSS,MED_ALL,profileName,MED_COMPACT,MED_MAILLE, + typmai3[(int)(*iter).getType()],numdt,(char *)"",dt,numo); + delete [] profileName; + delete [] nommaa; + pt+=(*iter).getNbOfTuple()*nbComp; + } + MEDfermer(fid); +} + +/*! + * This method performs the classical job for fields before any values setting. + */ +med_idt MEDLoaderNS::appendFieldSimpleAtt(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, med_int& numdt, med_int& numo, med_float& dt) { med_idt fid=MEDouvrir((char *)fileName,MED_LECTURE_ECRITURE); int nbComp=f->getNumberOfComponents(); char *comp=MEDLoaderBase::buildEmptyString(nbComp*MED_TAILLE_PNOM); char *unit=MEDLoaderBase::buildEmptyString(nbComp*MED_TAILLE_PNOM); MEDchampCr(fid,(char *)f->getName(),MED_FLOAT64,comp,unit,nbComp); - med_int numdt,numo; - med_float dt; + ParaMEDMEM::TypeOfTimeDiscretization td=f->getTimeDiscretization(); if(td==ParaMEDMEM::NO_TIME) { @@ -1430,13 +1529,24 @@ void MEDLoaderNS::appendFieldDirectly(const char *fileName, ParaMEDMEM::MEDCoupl numdt=(med_int)tmp1; numo=(med_int)tmp2; dt=(med_float)tmp0; } + delete [] comp; + delete [] unit; + return fid; +} + +void MEDLoaderNS::appendFieldDirectly(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f) +{ + med_int numdt,numo; + med_float dt; + int nbComp=f->getNumberOfComponents(); + med_idt fid=appendFieldSimpleAtt(fileName,f,numdt,numo,dt); const double *pt=f->getArray()->getConstPointer(); switch(f->getTypeOfField()) { case ParaMEDMEM::ON_CELLS: { std::list split; - prepareCellFieldDoubleForWriting(f,split); + prepareCellFieldDoubleForWriting(f,0,split); for(std::list::const_iterator iter=split.begin();iter!=split.end();iter++) { char *nommaa=MEDLoaderBase::buildEmptyString(MED_TAILLE_NOM); @@ -1462,12 +1572,14 @@ void MEDLoaderNS::appendFieldDirectly(const char *fileName, ParaMEDMEM::MEDCoupl default: throw INTERP_KERNEL::Exception("Not managed this type of FIELD !"); } - delete [] comp; - delete [] unit; MEDfermer(fid); } -void MEDLoaderNS::prepareCellFieldDoubleForWriting(const ParaMEDMEM::MEDCouplingFieldDouble *f, std::list& split) +/*! + * This method splits field 'f' into types to be ready for writing. + * @param cellIdsPerType this parameter can be 0 if not in profile mode. If it is != 0 this array is of size f->getMesh()->getNumberOfCells(). + */ +void MEDLoaderNS::prepareCellFieldDoubleForWriting(const ParaMEDMEM::MEDCouplingFieldDouble *f, const int *cellIdsPerType, std::list& split) { int nbComp=f->getNumberOfComponents(); const MEDCouplingMesh *mesh=f->getMesh(); @@ -1480,11 +1592,18 @@ void MEDLoaderNS::prepareCellFieldDoubleForWriting(const ParaMEDMEM::MEDCoupling const int *conn=meshC->getNodalConnectivity()->getConstPointer(); int nbOfCells=meshC->getNumberOfCells(); INTERP_KERNEL::NormalizedCellType curType; + const int *wCellIdsPT=cellIdsPerType; for(const int *pt=connI;pt!=connI+nbOfCells;) { curType=(INTERP_KERNEL::NormalizedCellType)conn[*pt]; const int *pt2=std::find_if(pt+1,connI+nbOfCells,ConnReaderML(conn,(int)curType)); - split.push_back(MEDLoader::MEDFieldDoublePerCellType(curType,0,nbComp,pt2-pt)); + if(!cellIdsPerType) + split.push_back(MEDLoader::MEDFieldDoublePerCellType(curType,0,nbComp,pt2-pt,0)); + else + { + split.push_back(MEDLoader::MEDFieldDoublePerCellType(curType,0,nbComp,pt2-pt,wCellIdsPT)); + wCellIdsPT+=std::distance(pt,pt2); + } pt=pt2; } } @@ -1521,6 +1640,55 @@ void MEDLoaderNS::writeFieldTryingToFitExistingMesh(const char *fileName, ParaME throw INTERP_KERNEL::Exception(oss.str().c_str()); } MEDCouplingUMesh *m=MEDLoader::ReadUMeshFromFile(fileName,f->getMesh()->getName(),f2); + MEDCouplingUMesh *m2=MEDCouplingUMesh::mergeUMeshes(m,(MEDCouplingUMesh *)f->getMesh()); + bool areNodesMerged; + DataArrayInt *da=m2->mergeNodes(MEDLoader::_EPS_FOR_NODE_COMP,areNodesMerged); + if(!areNodesMerged || m2->getNumberOfNodes()!=m->getNumberOfNodes()) + { + da->decrRef(); + m2->decrRef(); + m->decrRef(); + std::ostringstream oss; oss << "Nodes in already written mesh \"" << f->getMesh()->getName() << "\" in file \"" << fileName << "\" does not fit coordinates of unstructured grid f->getMesh() !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } + switch(f->getTypeOfField()) + { + case ParaMEDMEM::ON_CELLS: + { + da->decrRef(); + DataArrayInt *da2=m2->zipConnectivityTraducer(MEDLoader::_COMP_FOR_CELL); + if(m2->getNumberOfCells()!=m->getNumberOfCells()) + { + da2->decrRef(); + m2->decrRef(); + m->decrRef(); + std::ostringstream oss1; oss1 << "Cells in already written mesh \"" << f->getMesh()->getName() << "\" in file \"" << fileName << "\" does not fit connectivity of unstructured grid f->getMesh() !"; + throw INTERP_KERNEL::Exception(oss1.str().c_str()); + } + da=m2->convertCellArrayPerGeoType(da2); + DataArrayInt *da3=da->substr(m2->getNumberOfCells()); + da2->decrRef(); + da2=m2->convertCellArrayPerGeoType(da3); + da3->decrRef(); + appendCellProfileField(fileName,f,da2->getConstPointer()); + da2->decrRef(); + break; + } + case ParaMEDMEM::ON_NODES: + { + appendNodeProfileField(fileName,f,da->getConstPointer()+m->getNumberOfNodes()); + break; + } + default: + { + da->decrRef(); + m2->decrRef(); + m->decrRef(); + throw INTERP_KERNEL::Exception("Not implemented other profile fitting from already written mesh for fields than on NODES and on CELLS."); + } + } + da->decrRef(); + m2->decrRef(); m->decrRef(); } diff --git a/src/MEDLoader/MEDLoader.hxx b/src/MEDLoader/MEDLoader.hxx index 965fa0584..83dc10ecb 100644 --- a/src/MEDLoader/MEDLoader.hxx +++ b/src/MEDLoader/MEDLoader.hxx @@ -37,7 +37,7 @@ namespace ParaMEDMEM class MEDLOADER_EXPORT MEDLoader { -public: + public: class MEDConnOfOneElemType { public: @@ -63,20 +63,24 @@ public: class MEDFieldDoublePerCellType { public: - MEDFieldDoublePerCellType(INTERP_KERNEL::NormalizedCellType type, double *values, int ncomp, int ntuple); + MEDFieldDoublePerCellType(INTERP_KERNEL::NormalizedCellType type, double *values, int ncomp, int ntuple, const int *cellIdPerType); INTERP_KERNEL::NormalizedCellType getType() const { return _type; } int getNbComp() const { return _ncomp; } int getNbOfTuple() const { return _ntuple; } int getNbOfValues() const { return _ncomp*_ntuple; } double *getArray() const { return _values; } + const std::vector& getCellIdPerType() const { return _cell_id_per_type; } void releaseArray(); private: int _ntuple; int _ncomp; double *_values; + std::vector _cell_id_per_type; INTERP_KERNEL::NormalizedCellType _type; }; // + static void setEpsilonForNodeComp(double val); + static void setCompPolicyForCell(int val); static std::vector GetMeshNames(const char *fileName); static std::vector GetMeshGroupsNames(const char *fileName, const char *meshName); static std::vector GetMeshFamilyNames(const char *fileName, const char *meshName); @@ -97,8 +101,11 @@ public: static void WriteUMeshes(const char *fileName, const char *meshName, const std::vector& meshes, bool writeFromScratch); static void WriteField(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f, bool writeFromScratch); static void WriteFieldUsingAlreadyWrittenMesh(const char *fileName, ParaMEDMEM::MEDCouplingFieldDouble *f); -private: + private: MEDLoader(); + public: + static double _EPS_FOR_NODE_COMP; + static int _COMP_FOR_CELL; }; #endif diff --git a/src/MEDLoader/Test/MEDLoaderTest.cxx b/src/MEDLoader/Test/MEDLoaderTest.cxx index 1ce22526d..341e52cee 100644 --- a/src/MEDLoader/Test/MEDLoaderTest.cxx +++ b/src/MEDLoader/Test/MEDLoaderTest.cxx @@ -331,6 +331,9 @@ void MEDLoaderTest::testFieldProfilRW1() { const char fileName[]="file12.med"; MEDCouplingUMesh *mesh1=build3DMesh_1(); + bool b; + DataArrayInt *da=mesh1->mergeNodes(1e-12,b); + da->decrRef(); MEDLoader::WriteUMesh(fileName,mesh1,true); const int part1[5]={1,2,4,13,15}; MEDCouplingUMesh *mesh2=(MEDCouplingUMesh *)mesh1->buildPartOfMySelf(part1,part1+5,true); @@ -352,6 +355,12 @@ void MEDLoaderTest::testFieldProfilRW1() f1->checkCoherency(); // MEDLoader::WriteField(fileName,f1,false);//<- false important for the test + // + MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldDoubleCell(fileName,f1->getMesh()->getName(),0,f1->getName(),2,7); + f2->checkCoherency(); + CPPUNIT_ASSERT(f1->isEqual(f2,1e-12,1e-12)); + // + f2->decrRef(); f1->decrRef(); mesh1->decrRef(); mesh2->decrRef(); -- 2.39.2