From: abn Date: Mon, 30 Nov 2015 12:48:19 +0000 (+0100) Subject: OverlapDEC: new test for default values and multi-compo fields. X-Git-Tag: simple_cmake~4^2~3 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=8baa4e58323c1dd796265204111fc5c7c6f8b77a;p=tools%2Fmedcoupling.git OverlapDEC: new test for default values and multi-compo fields. --- diff --git a/src/ParaMEDMEM/OverlapDEC.cxx b/src/ParaMEDMEM/OverlapDEC.cxx index 6c2cfcd1b..12374f912 100644 --- a/src/ParaMEDMEM/OverlapDEC.cxx +++ b/src/ParaMEDMEM/OverlapDEC.cxx @@ -299,7 +299,7 @@ namespace ParaMEDMEM _interpolation_matrix->addContribution(src,srcIds,srcMeth,(*it).first,trg,trgIds,trgMeth,(*it).second); } _interpolation_matrix->prepare(_locator->getProcsToSendFieldData()); - _interpolation_matrix->computeDeno(); + _interpolation_matrix->computeSurfacesAndDeno(); } void OverlapDEC::attachSourceLocalField(ParaFIELD *field, bool ownPt) diff --git a/src/ParaMEDMEM/OverlapInterpolationMatrix.cxx b/src/ParaMEDMEM/OverlapInterpolationMatrix.cxx index 52f89e47e..7b84d60a4 100644 --- a/src/ParaMEDMEM/OverlapInterpolationMatrix.cxx +++ b/src/ParaMEDMEM/OverlapInterpolationMatrix.cxx @@ -241,7 +241,6 @@ namespace ParaMEDMEM _mapping.addContributionST(sparse_matrix_part,srcIds,srcProcId,trgIds,trgProcId); } - /*! * 'procsToSendField' gives the list of procs field data has to be sent to. * See OverlapElementLocator::computeBoundingBoxesAndTodoList() @@ -254,12 +253,29 @@ namespace ParaMEDMEM _mapping.prepare(procsToSendField,0); } - void OverlapInterpolationMatrix::computeDeno() + void OverlapInterpolationMatrix::computeSurfacesAndDeno() { if(_target_field->getField()->getNature()==ConservativeVolumic) _mapping.computeDenoConservativeVolumic(_target_field->getField()->getNumberOfTuplesExpected()); else - throw INTERP_KERNEL::Exception("Policy Not implemented yet : only ConservativeVolumic defined !"); + throw INTERP_KERNEL::Exception("OverlapDEC: Policy not implemented yet: only ConservativeVolumic!"); +// { +// if(_target_field->getField()->getNature()==RevIntegral) +// { +// MEDCouplingAutoRefCountObjectPtr f; +// int orient = getOrientation(); // From InterpolationOptions inheritance +// if(orient == 2) // absolute areas +// f = _target_support->getMeasureField(true); +// else +// if(orient == 0) // relative areas +// f = _target_support->getMeasureField(false); +// else +// throw INTERP_KERNEL::Exception("OverlapDEC: orientation policy not impl. yet!"); +// _mapping.computeDenoRevIntegral(*f->getArray()); +// } +// else +// throw INTERP_KERNEL::Exception("OverlapDEC: Policy not implemented yet: only ConservativeVolumic and RevIntegral defined!"); +// } } void OverlapInterpolationMatrix::multiply(double default_val) @@ -272,11 +288,6 @@ namespace ParaMEDMEM _mapping.transposeMultiply(_target_field->getField(),_source_field->getField()); } -// bool OverlapInterpolationMatrix::isSurfaceComputationNeeded(const std::string& method) const -// { -// return method=="P0"; -// } - void OverlapInterpolationMatrix::TransposeMatrix(const std::vector& matIn, int nbColsMatIn, std::vector& matOut) { diff --git a/src/ParaMEDMEM/OverlapInterpolationMatrix.hxx b/src/ParaMEDMEM/OverlapInterpolationMatrix.hxx index 30319196a..de7fb905b 100644 --- a/src/ParaMEDMEM/OverlapInterpolationMatrix.hxx +++ b/src/ParaMEDMEM/OverlapInterpolationMatrix.hxx @@ -57,7 +57,7 @@ namespace ParaMEDMEM void prepare(const std::vector< int > & procsToSendField); - void computeDeno(); + void computeSurfacesAndDeno(); void multiply(double default_val); diff --git a/src/ParaMEDMEM/OverlapMapping.cxx b/src/ParaMEDMEM/OverlapMapping.cxx index 28e1f4b3e..1ef168ff4 100644 --- a/src/ParaMEDMEM/OverlapMapping.cxx +++ b/src/ParaMEDMEM/OverlapMapping.cxx @@ -156,8 +156,7 @@ void OverlapMapping::prepare(const std::vector< int >& procsToSendField, int nbO updateZipSourceIdsForMultiply(); //finish to fill _the_matrix_st with already in place matrix in _matrixes_st (local computation) finishToFillFinalMatrixST(); - // Prepare proc lists for future field data exchange (mutliply()): - _proc_ids_to_recv_vector_st = _the_matrix_st_source_proc_id; + // Prepare proc list for future field data exchange (mutliply()): _proc_ids_to_send_vector_st = procsToSendField; // Make some space on local proc: _matrixes_st.clear(); @@ -167,54 +166,79 @@ void OverlapMapping::prepare(const std::vector< int >& procsToSendField, int nbO #endif } -///** -// * Fill the members _proc_ids_to_send_vector_st and _proc_ids_to_recv_vector_st -// * indicating for each proc to which proc it should send its source field data -// * and from which proc it should receive source field data. -// * -// * Should be called once finishToFillFinalMatrixST() has been executed, i.e. when the -// * local matrices are completely ready. +///*! +// * Compute denominators for IntegralGlobConstraint interp. +// * TO BE REVISED: needs another communication since some bits are held non locally // */ -//void OverlapMapping::fillProcToSendRcvForMultiply(const std::vector< int >& procsToSendField) +//void OverlapMapping::computeDenoGlobConstraint() //{ -//// _proc_ids_to_send_vector_st.clear(); -// _proc_ids_to_recv_vector_st.clear(); -// // Receiving side is easy - just inspect non-void terms in the final matrix -// int i=0; -// std::vector< std::vector< SparseDoubleVec > >::const_iterator it; -// std::vector< SparseDoubleVec >::const_iterator it2; -// for(it=_the_matrix_st.begin(); it!=_the_matrix_st.end(); it++, i++) -// _proc_ids_to_recv_vector_st.push_back(_the_matrix_st_source_proc_id[i]); -// -// // Sending side was computed in OverlapElementLocator::computeBoundingBoxesAndTodoList() -// _proc_ids_to_send_vector_st = procsToSendField; +// _the_deno_st.clear(); +// std::size_t sz1=_the_matrix_st.size(); +// _the_deno_st.resize(sz1); +// for(std::size_t i=0;i deno(nbOfTuplesTrg); + // Fills in the vector indexed by target cell ID: for(std::size_t i=0;i& mat=_the_matrix_st[i]; int curSrcId=_the_matrix_st_source_proc_id[i]; std::vector::iterator isItem1=std::find(_sent_trg_proc_st2.begin(),_sent_trg_proc_st2.end(),curSrcId); int rowId=0; - if(isItem1==_sent_trg_proc_st2.end() || curSrcId==myProcId)//item1 of step2 main algo. Simple, because rowId of mat are directly target ids. + if(isItem1==_sent_trg_proc_st2.end() || curSrcId==myProcId) // Local computation: simple, because rowId of mat are directly target cell ids. { for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++) for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) deno[rowId]+=(*it2).second; } - else - {//item0 of step2 main algo. More complicated. + else // matrix was received, remote computation + { std::vector::iterator fnd=isItem1;//std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),curSrcId); int locId=std::distance(_sent_trg_proc_st2.begin(),fnd); const DataArrayInt *trgIds=_sent_trg_ids_st2[locId]; @@ -249,7 +274,7 @@ void OverlapMapping::computeDenoConservativeVolumic(int nbOfTuplesTrg) deno[trgIds2[rowId]]+=(*it2).second; } } - // + // Broadcast the vector into a structure similar to the initial sparse matrix of numerators: for(std::size_t i=0;i receive 'interp source IDs' set of field values */ + const std::vector< int > & _proc_ids_to_recv_vector_st = _the_matrix_st_source_proc_id; if (procID == myProcID) nbrecv[procID] = 0; else diff --git a/src/ParaMEDMEM/OverlapMapping.hxx b/src/ParaMEDMEM/OverlapMapping.hxx index b30afe418..58dab05c4 100644 --- a/src/ParaMEDMEM/OverlapMapping.hxx +++ b/src/ParaMEDMEM/OverlapMapping.hxx @@ -26,6 +26,7 @@ #include #include +//#define DEC_DEBUG namespace ParaMEDMEM { @@ -52,7 +53,9 @@ namespace ParaMEDMEM void addContributionST(const std::vector< SparseDoubleVec >& matrixST, const DataArrayInt *srcIds, int srcProcId, const DataArrayInt *trgIds, int trgProcId); void prepare(const std::vector< int >& procsToSendField, int nbOfTrgElems); void computeDenoConservativeVolumic(int nbOfTuplesTrg); - void computeDenoGlobConstraint(); +// void computeDenoIntegralGlobConstraint(); +// void computeDenoIntegral(); + void computeDenoRevIntegral(const DataArrayDouble & targetAreas); // void multiply(const MEDCouplingFieldDouble *fieldInput, MEDCouplingFieldDouble *fieldOutput, double default_val) const; void transposeMultiply(const MEDCouplingFieldDouble *fieldInput, MEDCouplingFieldDouble *fieldOutput); @@ -125,7 +128,7 @@ namespace ParaMEDMEM //! Proc IDs to which data will be sent (originating this current proc) for matrix-vector computation std::vector< int > _proc_ids_to_send_vector_st; //! Proc IDs from which data will be received (on this current proc) for matrix-vector computation - std::vector< int > _proc_ids_to_recv_vector_st; + //std::vector< int > _proc_ids_to_recv_vector_st; // directly equal to _the_matrix_st_source_proc_id // Denominators (computed from the numerator matrix). As for _the_matrix_st it is paired with _the_matrix_st_source_proc_id std::vector< std::vector< SparseDoubleVec > > _the_deno_st; diff --git a/src/ParaMEDMEMTest/ParaMEDMEMTest.hxx b/src/ParaMEDMEMTest/ParaMEDMEMTest.hxx index 1273da391..b15b53fd1 100644 --- a/src/ParaMEDMEMTest/ParaMEDMEMTest.hxx +++ b/src/ParaMEDMEMTest/ParaMEDMEMTest.hxx @@ -50,8 +50,7 @@ class ParaMEDMEMTest : public CppUnit::TestFixture CPPUNIT_TEST(testOverlapDEC1); // 3 procs CPPUNIT_TEST(testOverlapDEC2); // 3 procs CPPUNIT_TEST(testOverlapDEC3); // 2 procs -// CPPUNIT_TEST(testOverlapDEC4); - + CPPUNIT_TEST(testOverlapDEC4); // 2 procs CPPUNIT_TEST(testSynchronousEqualInterpKernelWithoutInterpNativeDEC_2D);// 5 procs CPPUNIT_TEST(testSynchronousEqualInterpKernelWithoutInterpDEC_2D); // 5 procs @@ -111,10 +110,10 @@ public: void testInterpKernelDEC3DSurfEmptyBBox(); void testOverlapDEC1(); void testOverlapDEC2(); +// void testOverlapDEC2_bis(); void testOverlapDEC3(); +// void testOverlapDEC3_bis(); void testOverlapDEC4(); - void testOverlapDEC_LMEC_seq(); - void testOverlapDEC_LMEC_para(); #ifdef MED_ENABLE_FVM void testNonCoincidentDEC_2D(); void testNonCoincidentDEC_3D(); diff --git a/src/ParaMEDMEMTest/ParaMEDMEMTest_OverlapDEC.cxx b/src/ParaMEDMEMTest/ParaMEDMEMTest_OverlapDEC.cxx index d0f41d12c..63bfc9f9b 100644 --- a/src/ParaMEDMEMTest/ParaMEDMEMTest_OverlapDEC.cxx +++ b/src/ParaMEDMEMTest/ParaMEDMEMTest_OverlapDEC.cxx @@ -386,7 +386,7 @@ void ParaMEDMEMTest::testOverlapDEC1() /*! * Same as testOverlapDEC1() but with regular bounding boxes. If you're looking for a nice debug case, - * testOverlapDEC1() is more appropriate. + * testOverlapDEC1() is identical in terms of geometry and field values, and more appropriate. */ void ParaMEDMEMTest::testOverlapDEC2() { @@ -472,6 +472,7 @@ void prepareData2_buildOneSquare(MEDCouplingUMesh* & meshS_0, MEDCouplingUMesh* meshT_0->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT+3); meshT_0->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT+6); meshT_0->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT+9); + } /** @@ -487,7 +488,9 @@ void prepareData2_buildOneSquare(MEDCouplingUMesh* & meshS_0, MEDCouplingUMesh* void prepareData2(int rank, ProcessorGroup * grp, NatureOfField nature, MEDCouplingUMesh *& meshS, MEDCouplingUMesh *& meshT, ParaMESH*& parameshS, ParaMESH*& parameshT, - ParaFIELD*& parafieldS, ParaFIELD*& parafieldT) + ParaFIELD*& parafieldS, ParaFIELD*& parafieldT, + bool stripPartOfSource=false, + int fieldCompoNum=1) { MEDCouplingUMesh *meshS_0 = 0, *meshT_0 = 0; prepareData2_buildOneSquare(meshS_0, meshT_0); @@ -502,16 +505,26 @@ void prepareData2(int rank, ProcessorGroup * grp, NatureOfField nature, meshS_2->translate(tr2); std::vector vec; - vec.push_back(meshS_0);vec.push_back(meshS_1);vec.push_back(meshS_2); + vec.push_back(meshS_0);vec.push_back(meshS_1); + if (!stripPartOfSource) + vec.push_back(meshS_2); meshS = MEDCouplingUMesh::MergeUMeshes(vec); meshS_1->decrRef(); meshS_2->decrRef(); - ComponentTopology comptopo; + ComponentTopology comptopo(fieldCompoNum); parameshS=new ParaMESH(meshS, *grp,"source mesh"); parafieldS=new ParaFIELD(ON_CELLS,ONE_TIME,parameshS,comptopo); parafieldS->getField()->setNature(nature); double *valsS=parafieldS->getField()->getArray()->getPointer(); - valsS[0]=1.; valsS[1]=2.; valsS[2]=3.; + for(int i=0; i < fieldCompoNum; i++) + { + valsS[i] = 1. * (10^i); + valsS[fieldCompoNum+i] = 2. * (10^i); + if (!stripPartOfSource) + { + valsS[2*fieldCompoNum+i] = 3. * (10^i); + } + } // const double tr3[] = {0.0, -1.5}; @@ -541,12 +554,16 @@ void prepareData2(int rank, ProcessorGroup * grp, NatureOfField nature, meshS = MEDCouplingUMesh::MergeUMeshes(vec); meshS_3->decrRef(); meshS_4->decrRef(); - ComponentTopology comptopo; + ComponentTopology comptopo(fieldCompoNum); parameshS=new ParaMESH(meshS, *grp,"source mesh"); parafieldS=new ParaFIELD(ON_CELLS,ONE_TIME,parameshS,comptopo); parafieldS->getField()->setNature(nature); double *valsS=parafieldS->getField()->getArray()->getPointer(); - valsS[0]=4.; valsS[1]=5.; + for(int i=0; i < fieldCompoNum; i++) + { + valsS[i] = 4. * (10^i); + valsS[fieldCompoNum+i] = 5. * (10^i); + } // const double tr5[] = {1.5, 0.0}; @@ -595,12 +612,6 @@ void ParaMEDMEMTest::testOverlapDEC3() ParaFIELD* parafieldS=0, *parafieldT=0; prepareData2(rank, grp, ConservativeVolumic, meshS, meshT, parameshS, parameshT, parafieldS, parafieldT); -// if (rank == 1) -// { -// int i=1, j=0; -// while (i!=0) -// j=2; -// } dec.attachSourceLocalField(parafieldS); dec.attachTargetLocalField(parafieldT); @@ -637,8 +648,89 @@ void ParaMEDMEMTest::testOverlapDEC3() } /*! - * Test default values and multi-component fields + * Tests: + * - default value + * - multi-component fields */ void ParaMEDMEMTest::testOverlapDEC4() { + int size, rank; + MPI_Comm_size(MPI_COMM_WORLD,&size); + MPI_Comm_rank(MPI_COMM_WORLD,&rank); + + int nproc = 2; + if (size != nproc) return ; + std::set procs; + for (int i=0; igetField(); + if(rank==0) + { + CPPUNIT_ASSERT_EQUAL(8, resField->getNumberOfTuples()); + for(int i=0;i<4;i++) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,resField->getArray()->getIJ(i*2,0),1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(10.0,resField->getArray()->getIJ(i*2+1,0),1e-12); + } + for(int i=4;i<8;i++) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0,resField->getArray()->getIJ(i*2,0),1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(40.0,resField->getArray()->getIJ(i*2+1,0),1e-12); + } + } + if(rank==1) + { + CPPUNIT_ASSERT_EQUAL(12, resField->getNumberOfTuples()); + for(int i=0;i<4;i++) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(2.0,resField->getArray()->getIJ(i*2,0),1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(20.0,resField->getArray()->getIJ(i*2+1,0),1e-12); + } + // Default value should be here: + for(int i=4;i<8;i++) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(defVal,resField->getArray()->getIJ(i*2,0),1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(defVal,resField->getArray()->getIJ(i*2+1,0),1e-12); + } + for(int i=8;i<12;i++) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(5.0,resField->getArray()->getIJ(i*2,0),1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(50.0,resField->getArray()->getIJ(i*2+1,0),1e-12); + } + } + delete parafieldS; + delete parafieldT; + delete parameshS; + delete parameshT; + meshS->decrRef(); + meshT->decrRef(); + + MPI_Barrier(MPI_COMM_WORLD); + } +