From 647dcbb377639c94c4f8f29aba4ca57c6d3ea31d Mon Sep 17 00:00:00 2001 From: ageay Date: Thu, 10 Oct 2013 06:48:42 +0000 Subject: [PATCH] DataArrayDouble::getIdsNotInRange, DataArrayInt::getIdsNotInRange --- src/MEDCoupling/MEDCouplingMemArray.cxx | 70 ++++++++++++++++--- src/MEDCoupling/MEDCouplingMemArray.hxx | 2 + src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 38 ++++++++++ src/MEDCoupling_Swig/MEDCouplingMemArray.i | 4 ++ 4 files changed, 106 insertions(+), 8 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index fec8fe9db..ed638f94e 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -4469,6 +4469,8 @@ DataArrayDoubleIterator *DataArrayDouble::iterator() * needed. * \throw If \a this->getNumberOfComponents() != 1. * + * \sa DataArrayDouble::getIdsNotInRange + * * \ref cpp_mcdataarraydouble_getidsinrange "Here is a C++ example".
* \ref py_mcdataarraydouble_getidsinrange "Here is a Python example". */ @@ -4477,15 +4479,41 @@ DataArrayInt *DataArrayDouble::getIdsInRange(double vmin, double vmax) const checkAllocated(); if(getNumberOfComponents()!=1) throw INTERP_KERNEL::Exception("DataArrayDouble::getIdsInRange : this must have exactly one component !"); - const double *cptr=getConstPointer(); - MEDCouplingAutoRefCountObjectPtr ret=DataArrayInt::New(); ret->alloc(0,1); - int nbOfTuples=getNumberOfTuples(); + const double *cptr(begin()); + MEDCouplingAutoRefCountObjectPtr ret(DataArrayInt::New()); ret->alloc(0,1); + int nbOfTuples(getNumberOfTuples()); for(int i=0;i=vmin && *cptr<=vmax) ret->pushBackSilent(i); return ret.retn(); } +/*! + * Returns a new DataArrayInt contating indices of tuples of \a this one-dimensional + * array whose values are not within a given range. Textual data is not copied. + * \param [in] vmin - a lowest not acceptable value (excluded). + * \param [in] vmax - a greatest not acceptable value (excluded). + * \return DataArrayInt * - the new instance of DataArrayInt. + * The caller is to delete this result array using decrRef() as it is no more + * needed. + * \throw If \a this->getNumberOfComponents() != 1. + * + * \sa DataArrayDouble::getIdsInRange + */ +DataArrayInt *DataArrayDouble::getIdsNotInRange(double vmin, double vmax) const +{ + checkAllocated(); + if(getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DataArrayDouble::getIdsNotInRange : this must have exactly one component !"); + const double *cptr(begin()); + MEDCouplingAutoRefCountObjectPtr ret(DataArrayInt::New()); ret->alloc(0,1); + int nbOfTuples(getNumberOfTuples()); + for(int i=0;ivmax) + ret->pushBackSilent(i); + return ret.retn(); +} + /*! * Returns a new DataArrayDouble by concatenating two given arrays, so that (1) the number * of tuples in the result array is a sum of the number of tuples of given arrays and (2) @@ -9016,29 +9044,55 @@ void DataArrayInt::applyModulus(int val) * \param [in] vmin begin of range. This value is included in range (included). * \param [in] vmax end of range. This value is \b not included in range (excluded). * \return a newly allocated data array that the caller should deal with. + * + * \sa DataArrayInt::getIdsNotInRange */ DataArrayInt *DataArrayInt::getIdsInRange(int vmin, int vmax) const { checkAllocated(); if(getNumberOfComponents()!=1) throw INTERP_KERNEL::Exception("DataArrayInt::getIdsInRange : this must have exactly one component !"); - const int *cptr=getConstPointer(); - MEDCouplingAutoRefCountObjectPtr ret=DataArrayInt::New(); ret->alloc(0,1); - int nbOfTuples=getNumberOfTuples(); + const int *cptr(begin()); + MEDCouplingAutoRefCountObjectPtr ret(DataArrayInt::New()); ret->alloc(0,1); + int nbOfTuples(getNumberOfTuples()); for(int i=0;i=vmin && *cptrpushBackSilent(i); return ret.retn(); } +/*! + * This method works only on data array with one component. + * This method returns a newly allocated array storing stored ascendantly tuple ids in \b this so that + * this[*id] \b not in [\b vmin,\b vmax) + * + * \param [in] vmin begin of range. This value is \b not included in range (excluded). + * \param [in] vmax end of range. This value is included in range (included). + * \return a newly allocated data array that the caller should deal with. + * + * \sa DataArrayInt::getIdsInRange + */ +DataArrayInt *DataArrayInt::getIdsNotInRange(int vmin, int vmax) const +{ + checkAllocated(); + if(getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DataArrayInt::getIdsNotInRange : this must have exactly one component !"); + const int *cptr(getConstPointer()); + MEDCouplingAutoRefCountObjectPtr ret(DataArrayInt::New()); ret->alloc(0,1); + int nbOfTuples(getNumberOfTuples()); + for(int i=0;i=vmax) + ret->pushBackSilent(i); + return ret.retn(); +} + /*! * This method works only on data array with one component. * This method checks that all ids in \b this are in [ \b vmin, \b vmax ). If there is at least one element in \a this not in [ \b vmin, \b vmax ) an exception will be thrown. * * \param [in] vmin begin of range. This value is included in range (included). * \param [in] vmax end of range. This value is \b not included in range (excluded). - * \return if all ids in \a this are so that (*this)[i]==i for all i in [ 0, \c this->getNumberOfTuples() ). - */ + * \return if all ids in \a this are so that (*this)[i]==i for all i in [ 0, \c this->getNumberOfTuples() ). */ bool DataArrayInt::checkAllIdsInRange(int vmin, int vmax) const { checkAllocated(); diff --git a/src/MEDCoupling/MEDCouplingMemArray.hxx b/src/MEDCoupling/MEDCouplingMemArray.hxx index 1fa5194fb..8280413d2 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.hxx +++ b/src/MEDCoupling/MEDCouplingMemArray.hxx @@ -345,6 +345,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void applyFuncFast32(const char *func); MEDCOUPLING_EXPORT void applyFuncFast64(const char *func); MEDCOUPLING_EXPORT DataArrayInt *getIdsInRange(double vmin, double vmax) const; + MEDCOUPLING_EXPORT DataArrayInt *getIdsNotInRange(double vmin, double vmax) const; MEDCOUPLING_EXPORT static DataArrayDouble *Aggregate(const DataArrayDouble *a1, const DataArrayDouble *a2); MEDCOUPLING_EXPORT static DataArrayDouble *Aggregate(const std::vector& arr); MEDCOUPLING_EXPORT static DataArrayDouble *Meld(const DataArrayDouble *a1, const DataArrayDouble *a2); @@ -558,6 +559,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void applyPow(int val); MEDCOUPLING_EXPORT void applyRPow(int val); MEDCOUPLING_EXPORT DataArrayInt *getIdsInRange(int vmin, int vmax) const; + MEDCOUPLING_EXPORT DataArrayInt *getIdsNotInRange(int vmin, int vmax) const; MEDCOUPLING_EXPORT bool checkAllIdsInRange(int vmin, int vmax) const; MEDCOUPLING_EXPORT static DataArrayInt *Aggregate(const DataArrayInt *a1, const DataArrayInt *a2, int offsetA2); MEDCOUPLING_EXPORT static DataArrayInt *Aggregate(const std::vector& arr); diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index ee96fd4f6..b820e9192 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -13907,6 +13907,44 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertTrue(DataArrayDouble([0.58093333350930543,0.58093333350930543]).isEqual(m.getMeasureField(False).getArray(),1e-12)) pass + def testSwig2Hexa8HavingFacesWarped1(self): + """ This test is bases on a "error" of interpolation detected. After investigation cell #3 of src is warped that leads to the fact that when trg is + intersected with src the sum of intersection volume is greater than the volume of the trg cell. + A test that can be done is to split the cell #3 of src into tetrohedrons and by summing all the volumes it does not fit the volume computed of cell#3 unsplitted (expect for + GENERAL_24). + """ + srcCoo=DataArrayDouble([0.15694071546650565,0.09383333333333337,6.920842121738133,0.15774332475430292,0.185486666666667,6.920682472824616,0.1585459340420992,0.27713999999999994,6.9205228239111,0.07427195882345167,0.05782666666666668,6.937285959830335,0.06343673343819695,0.11347333333333297,6.939441220162809,0.05260150805294228,0.16911999999999996,6.941596480495282,0.014076262238703396,0.04800666666666667,6.949259628344076,0.014076262238703396,0.07092000000000007,6.949259628344076,0.15407499632681992,0.09383333333333338,6.897607484780063,0.15489234394181514,0.18548666666666702,6.897567331066572,0.15570969155680933,0.27714,6.897527177353081,0.06988819198237989,0.05782666666666669,6.901743317269663,0.05885399917995321,0.11347333333333298,6.9022853924017955,0.047819806377526586,0.16912,6.902827467533927,0.0085871208577874,0.048006666666666684,6.9047548457815076,0.0085871208577874,0.07092000000000008,6.9047548457815076,0.153883333333333,0.09383333333333338,6.820902,0.154701666666667,0.18548666666666702,6.820902,0.15551999999999996,0.27714,6.820902,0.06959499999999999,0.05782666666666669,6.820902,0.058547499999999975,0.11347333333333298,6.820902,0.04749999999999999,0.16912,6.820902],22,3) + src=MEDCouplingUMesh("TBmesh3D",3) ; src.setCoords(srcCoo) + src.allocateCells() + src.insertNextCell(NORM_HEXA8,[0,1,4,3,8,9,12,11]) + src.insertNextCell(NORM_HEXA8,[1,2,5,4,9,10,13,12]) + src.insertNextCell(NORM_HEXA8,[4,5,7,6,12,13,15,14]) + src.insertNextCell(NORM_HEXA8,[8,9,12,11,16,17,20,19]) + src.insertNextCell(NORM_HEXA8,[9,10,13,12,17,18,21,20]) + src.checkCoherency2() + # trg is useless here but I keep it in case of MEDCouplingRemapper were expected to do something about warped NORM_HEXA8 + trgCoo=DataArrayDouble([0.0960891897852753,0.105088620541845,6.8598,0.0599574480546212,0.118434267436059,6.8598,0.113514510609589,0.14874473653263,6.8598,0.0831322609794463,0.167319109733883,6.8598,0.0960891897852753,0.105088620541845,6.92146666666667,0.0599574480546212,0.118434267436059,6.92146666666667,0.113514510609589,0.14874473653263,6.92146666666667,0.0831322609794463,0.167319109733883,6.92146666666667],8,3) + trg=MEDCouplingUMesh("MESH",3) ; trg.setCoords(trgCoo) + trg.allocateCells() + trg.insertNextCell(NORM_HEXA8,[0,1,3,2,4,5,7,6]) + # + srcFace=src.buildDescendingConnectivity()[0] + conn=MEDCoupling1SGTUMesh(srcFace).getNodalConnectivity() ; conn.rearrange(4) + eqFaces=srcFace.computePlaneEquationOf3DFaces() + nodeIdInCell=3 + e=(srcFace.getCoords()[conn[:,nodeIdInCell]]*eqFaces[:,:-1]).sumPerTuple()+eqFaces[:,3]# e represent the error between the expected 'a*X+b*Y+c*Z+d' in eqFaces and 0. Closer e to 0. is closer the 4th point is to the plane built with the 3 first points + lambd=-e/(eqFaces[:,:3]**2).sumPerTuple() + pts=lambd*eqFaces[:,:-1]+srcFace.getCoords()[conn[:,nodeIdInCell]]#pts represent the projection of the last points of each NORM_QUAD4 to the plane defined by the 3 first points of the NORM_QUAD4 cell + shouldBeZero=(pts*eqFaces[:,:-1]).sumPerTuple()+eqFaces[:,3]# this line is useless only to be sure that pts are on the plane. + check=(pts-srcFace.getCoords()[conn[:,nodeIdInCell]]).magnitude() # check contains the distance of the last point to its plane + idsToTest=check.getIdsNotInRange(0.,1e-10) + self.assertTrue(idsToTest.isEqual(DataArrayInt([17,18,19,20,22,23,24]))) + idsToTest2=idsToTest.getIdsNotInRange(18,22) + self.assertTrue(idsToTest2.isEqual(DataArrayInt([0,4,5,6]))) + idsToTest2.rearrange(2) + self.assertTrue(idsToTest2.sumPerTuple().isEqual(DataArrayInt([4,11]))) + pass + def setUp(self): pass pass diff --git a/src/MEDCoupling_Swig/MEDCouplingMemArray.i b/src/MEDCoupling_Swig/MEDCouplingMemArray.i index f3e7a7655..9039ac147 100644 --- a/src/MEDCoupling_Swig/MEDCouplingMemArray.i +++ b/src/MEDCoupling_Swig/MEDCouplingMemArray.i @@ -70,6 +70,7 @@ %newobject ParaMEDMEM::DataArrayInt::sumPerTuple; %newobject ParaMEDMEM::DataArrayInt::negate; %newobject ParaMEDMEM::DataArrayInt::getIdsInRange; +%newobject ParaMEDMEM::DataArrayInt::getIdsNotInRange; %newobject ParaMEDMEM::DataArrayInt::Aggregate; %newobject ParaMEDMEM::DataArrayInt::AggregateIndexes; %newobject ParaMEDMEM::DataArrayInt::Meld; @@ -148,6 +149,7 @@ %newobject ParaMEDMEM::DataArrayDouble::changeNbOfComponents; %newobject ParaMEDMEM::DataArrayDouble::accumulatePerChunck; %newobject ParaMEDMEM::DataArrayDouble::getIdsInRange; +%newobject ParaMEDMEM::DataArrayDouble::getIdsNotInRange; %newobject ParaMEDMEM::DataArrayDouble::negate; %newobject ParaMEDMEM::DataArrayDouble::applyFunc; %newobject ParaMEDMEM::DataArrayDouble::applyFunc2; @@ -574,6 +576,7 @@ namespace ParaMEDMEM void applyFuncFast32(const char *func) throw(INTERP_KERNEL::Exception); void applyFuncFast64(const char *func) throw(INTERP_KERNEL::Exception); DataArrayInt *getIdsInRange(double vmin, double vmax) const throw(INTERP_KERNEL::Exception); + DataArrayInt *getIdsNotInRange(double vmin, double vmax) const throw(INTERP_KERNEL::Exception); static DataArrayDouble *Aggregate(const DataArrayDouble *a1, const DataArrayDouble *a2) throw(INTERP_KERNEL::Exception); static DataArrayDouble *Meld(const DataArrayDouble *a1, const DataArrayDouble *a2) throw(INTERP_KERNEL::Exception); static DataArrayDouble *Dot(const DataArrayDouble *a1, const DataArrayDouble *a2) throw(INTERP_KERNEL::Exception); @@ -2560,6 +2563,7 @@ namespace ParaMEDMEM void applyPow(int val) throw(INTERP_KERNEL::Exception); void applyRPow(int val) throw(INTERP_KERNEL::Exception); DataArrayInt *getIdsInRange(int vmin, int vmax) const throw(INTERP_KERNEL::Exception); + DataArrayInt *getIdsNotInRange(int vmin, int vmax) const throw(INTERP_KERNEL::Exception); bool checkAllIdsInRange(int vmin, int vmax) const throw(INTERP_KERNEL::Exception); static DataArrayInt *Aggregate(const DataArrayInt *a1, const DataArrayInt *a2, int offsetA2) throw(INTERP_KERNEL::Exception); static DataArrayInt *Meld(const DataArrayInt *a1, const DataArrayInt *a2) throw(INTERP_KERNEL::Exception); -- 2.30.2