_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)
_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()
_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<MEDCouplingFieldDouble> 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)
_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<SparseDoubleVec >& matIn,
int nbColsMatIn, std::vector<SparseDoubleVec >& matOut)
{
void prepare(const std::vector< int > & procsToSendField);
- void computeDeno();
+ void computeSurfacesAndDeno();
void multiply(double default_val);
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();
#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<sz1;i++)
+// {
+// std::size_t sz2=_the_matrix_st[i].size();
+// _the_deno_st[i].resize(sz2);
+// for(std::size_t j=0;j<sz2;j++)
+// {
+// double sum=0;
+// SparseDoubleVec& mToFill=_the_deno_st[i][j];
+// const SparseDoubleVec& m=_the_matrix_st[i][j];
+// for(SparseDoubleVec::const_iterator it=m.begin();it!=m.end();it++)
+// sum+=(*it).second;
+// for(SparseDoubleVec::const_iterator it=m.begin();it!=m.end();it++)
+// mToFill[(*it).first]=sum;
+// }
+// }
+// printDenoMatrix();
//}
-/*!
- * Compute denominators for IntegralGlobConstraint interp.
- */
-void OverlapMapping::computeDenoGlobConstraint()
+///*! Compute integral denominators
+// * TO BE REVISED: needs another communication since some source areas are held non locally
+// */
+//void OverlapMapping::computeDenoIntegral()
+//{
+// _the_deno_st.clear();
+// std::size_t sz1=_the_matrix_st.size();
+// _the_deno_st.resize(sz1);
+// for(std::size_t i=0;i<sz1;i++)
+// {
+// std::size_t sz2=_the_matrix_st[i].size();
+// _the_deno_st[i].resize(sz2);
+// for(std::size_t j=0;j<sz2;j++)
+// {
+// SparseDoubleVec& mToFill=_the_deno_st[i][j];
+// for(SparseDoubleVec::const_iterator it=mToFill.begin();it!=mToFill.end();it++)
+// mToFill[(*it).first] = sourceAreas;
+// }
+// }
+// printDenoMatrix();
+//}
+
+/*! Compute rev integral denominators
+ */
+void OverlapMapping::computeDenoRevIntegral(const DataArrayDouble & targetAreas)
{
_the_deno_st.clear();
std::size_t sz1=_the_matrix_st.size();
_the_deno_st.resize(sz1);
+ const double * targetAreasP = targetAreas.getConstPointer();
for(std::size_t i=0;i<sz1;i++)
{
std::size_t sz2=_the_matrix_st[i].size();
_the_deno_st[i].resize(sz2);
for(std::size_t j=0;j<sz2;j++)
{
- double sum=0;
SparseDoubleVec& mToFill=_the_deno_st[i][j];
- const SparseDoubleVec& m=_the_matrix_st[i][j];
- for(SparseDoubleVec::const_iterator it=m.begin();it!=m.end();it++)
- sum+=(*it).second;
- for(SparseDoubleVec::const_iterator it=m.begin();it!=m.end();it++)
- mToFill[(*it).first]=sum;
+ SparseDoubleVec& mToIterate=_the_matrix_st[i][j];
+ for(SparseDoubleVec::const_iterator it=mToIterate.begin();it!=mToIterate.end();it++)
+ mToFill[(*it).first] = targetAreasP[j];
}
}
+// printDenoMatrix();
}
+
/*!
* Compute denominators for ConvervativeVolumic interp.
*/
std::size_t sz1=_the_matrix_st.size();
_the_deno_st.resize(sz1);
std::vector<double> deno(nbOfTuplesTrg);
+ // Fills in the vector indexed by target cell ID:
for(std::size_t i=0;i<sz1;i++)
{
const std::vector< SparseDoubleVec >& mat=_the_matrix_st[i];
int curSrcId=_the_matrix_st_source_proc_id[i];
std::vector<int>::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<int>::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];
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<sz1;i++)
{
int rowId=0;
}
}
+
/*!
* This method performs a transpose multiply of 'fieldInput' and put the result into 'fieldOutput'.
* 'fieldInput' is expected to be the sourcefield and 'fieldOutput' the targetfield.
* % else (=we did NOT compute the job, hence procID has, and knows the matrix)
* => 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
#include <vector>
#include <map>
+//#define DEC_DEBUG
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);
//! 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;
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
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();
/*!
* 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()
{
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);
+
}
/**
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);
meshS_2->translate(tr2);
std::vector<const MEDCouplingUMesh*> 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};
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};
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);
}
/*!
- * 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<int> procs;
+ for (int i=0; i<nproc; i++)
+ procs.insert(i);
+
+ CommInterface interface;
+ OverlapDEC dec(procs);
+ ProcessorGroup * grp = dec.getGroup();
+ MEDCouplingUMesh* meshS=0, *meshT=0;
+ ParaMESH* parameshS=0, *parameshT=0;
+ ParaFIELD* parafieldS=0, *parafieldT=0;
+
+ // As before, except than one of the source cell is removed, and that the field now has 2 components
+ prepareData2(rank, grp, ConservativeVolumic, meshS, meshT, parameshS, parameshT, parafieldS, parafieldT,
+ true, 2);
+// if (rank == 1)
+// {
+// int i=1, j=0;
+// while (i!=0)
+// j=2;
+// }
+
+ dec.attachSourceLocalField(parafieldS);
+ dec.attachTargetLocalField(parafieldT);
+ double defVal = -300.0;
+ dec.setDefaultValue(defVal);
+ dec.synchronize();
+ dec.sendRecvData(true);
+ //
+ MEDCouplingFieldDouble * resField = parafieldT->getField();
+ 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);
+
}
+
