1 // Copyright (C) 2007-2015 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #include "ParaMEDMEMTest.hxx"
21 #include <cppunit/TestAssert.h>
23 #include "CommInterface.hxx"
24 #include "ProcessorGroup.hxx"
25 #include "MPIProcessorGroup.hxx"
26 #include "Topology.hxx"
27 #include "OverlapDEC.hxx"
28 #include "ParaMESH.hxx"
29 #include "ParaFIELD.hxx"
30 #include "ComponentTopology.hxx"
32 #include "MEDCouplingUMesh.hxx"
38 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
39 #include "MEDLoader.hxx"
40 #include "MEDLoaderBase.hxx"
41 #include "MEDCouplingFieldDouble.hxx"
42 #include "MEDCouplingMemArray.hxx"
43 #include "MEDCouplingRemapper.hxx"
45 using namespace ParaMEDMEM;
47 typedef MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> MUMesh;
48 typedef MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> MFDouble;
50 //void ParaMEDMEMTest::testOverlapDEC_LMEC_seq()
52 // // T_SC_Trio_src.med -- "SupportOf_"
53 // // T_SC_Trio_dst.med -- "SupportOf_T_SC_Trio"
54 // // h_TH_Trio_src.med -- "SupportOf_"
55 // // h_TH_Trio_dst.med -- "SupportOf_h_TH_Trio"
56 // string rep("/export/home/adrien/support/antoine_LMEC/");
57 // string src_mesh_nam(rep + string("T_SC_Trio_src.med"));
58 // string tgt_mesh_nam(rep + string("T_SC_Trio_dst.med"));
59 //// string src_mesh_nam(rep + string("h_TH_Trio_src.med"));
60 //// string tgt_mesh_nam(rep + string("h_TH_Trio_dst.med"));
61 // MUMesh src_mesh=MEDLoader::ReadUMeshFromFile(src_mesh_nam,"SupportOf_",0);
62 // MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_T_SC_Trio",0);
63 //// MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_h_TH_Trio",0);
65 // MFDouble srcField = MEDCouplingFieldDouble::New(ON_CELLS, ONE_TIME);
66 // srcField->setMesh(src_mesh);
67 // DataArrayDouble * dad = DataArrayDouble::New(); dad->alloc(src_mesh->getNumberOfCells(),1);
68 // dad->fillWithValue(1.0);
69 // srcField->setArray(dad);
70 // srcField->setNature(ConservativeVolumic);
72 // MEDCouplingRemapper remap;
73 // remap.setOrientation(2); // always consider surface intersections as absolute areas.
74 // remap.prepare(src_mesh, tgt_mesh, "P0P0");
75 // MFDouble tgtField = remap.transferField(srcField, 1.0e+300);
76 // tgtField->setName("result");
77 // string out_nam(rep + string("adrien.med"));
78 // MEDLoader::WriteField(out_nam,tgtField, true);
79 // cout << "wrote: " << out_nam << "\n";
80 // double integ1 = 0.0, integ2 = 0.0;
81 // srcField->integral(true, &integ1);
82 // tgtField->integral(true, &integ2);
83 //// tgtField->reprQuickOverview(cout);
84 // CPPUNIT_ASSERT_DOUBLES_EQUAL(integ1,integ2,1e-8);
89 //void ParaMEDMEMTest::testOverlapDEC_LMEC_para()
91 // using namespace ParaMEDMEM;
95 // MPI_Comm_size(MPI_COMM_WORLD,&size);
96 // MPI_Comm_rank(MPI_COMM_WORLD,&rank);
98 // if (size != 1) return ;
101 // std::set<int> procs;
103 // for (int i=0; i<nproc; i++)
106 // CommInterface interface;
107 // OverlapDEC dec(procs);
109 // ParaMESH* parameshS=0;
110 // ParaMESH* parameshT=0;
111 // ParaFIELD* parafieldS=0;
112 // ParaFIELD* parafieldT=0;
113 // MFDouble srcField;
115 // // **** FILE LOADING
116 // // T_SC_Trio_src.med -- "SupportOf_"
117 // // T_SC_Trio_dst.med -- "SupportOf_T_SC_Trio"
118 // // h_TH_Trio_src.med -- "SupportOf_"
119 // // h_TH_Trio_dst.med -- "SupportOf_h_TH_Trio"
120 // string rep("/export/home/adrien/support/antoine_LMEC/");
121 // string src_mesh_nam(rep + string("T_SC_Trio_src.med"));
122 // string tgt_mesh_nam(rep + string("T_SC_Trio_dst.med"));
125 // MPI_Barrier(MPI_COMM_WORLD);
128 // // string src_mesh_nam(rep + string("h_TH_Trio_src.med"));
129 // // string tgt_mesh_nam(rep + string("h_TH_Trio_dst.med"));
130 // MUMesh src_mesh=MEDLoader::ReadUMeshFromFile(src_mesh_nam,"SupportOf_",0);
131 // MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_T_SC_Trio",0);
132 // // MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_h_TH_Trio",0);
135 // srcField = MEDCouplingFieldDouble::New(ON_CELLS, ONE_TIME);
136 // srcField->setMesh(src_mesh);
137 // DataArrayDouble * dad = DataArrayDouble::New(); dad->alloc(src_mesh->getNumberOfCells(),1);
138 // dad->fillWithValue(1.0);
139 // srcField->setArray(dad);
140 // srcField->setNature(ConservativeVolumic);
142 // ComponentTopology comptopo;
143 // parameshS = new ParaMESH(src_mesh,*dec.getGroup(),"source mesh");
144 // parafieldS = new ParaFIELD(ON_CELLS,ONE_TIME,parameshS,comptopo);
145 // parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
146 // parafieldS->getField()->setArray(dad);
149 // parameshT=new ParaMESH(tgt_mesh,*dec.getGroup(),"target mesh");
150 // parafieldT=new ParaFIELD(ON_CELLS,ONE_TIME,parameshT,comptopo);
151 // parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
152 // parafieldT->getField()->getArray()->fillWithValue(1.0e300);
155 // dec.setOrientation(2);
156 // dec.attachSourceLocalField(parafieldS);
157 // dec.attachTargetLocalField(parafieldT);
158 // dec.synchronize();
159 // dec.sendRecvData(true);
163 // double integ1 = 0.0, integ2 = 0.0;
164 // MEDCouplingFieldDouble * tgtField;
166 // srcField->integral(true, &integ1);
167 // tgtField = parafieldT->getField();
168 //// tgtField->reprQuickOverview(cout);
169 // tgtField->integral(true, &integ2);
170 // tgtField->setName("result");
171 // string out_nam(rep + string("adrien_para.med"));
172 // MEDLoader::WriteField(out_nam,tgtField, true);
173 // cout << "wrote: " << out_nam << "\n";
174 // CPPUNIT_ASSERT_DOUBLES_EQUAL(integ1,integ2,1e-8);
176 // delete parafieldS;
177 // delete parafieldT;
181 // MPI_Barrier(MPI_COMM_WORLD);
184 void prepareData1(int rank, ProcessorGroup * grp,
185 MEDCouplingUMesh *& meshS, MEDCouplingUMesh *& meshT,
186 ParaMESH*& parameshS, ParaMESH*& parameshT,
187 ParaFIELD*& parafieldS, ParaFIELD*& parafieldT)
191 const double coordsS[10]={0.,0.,0.5,0.,1.,0.,0.,0.5,0.5,0.5};
192 const double coordsT[6]={0.,0.,1.,0.,1.,1.};
193 meshS=MEDCouplingUMesh::New();
194 meshS->setMeshDimension(2);
195 DataArrayDouble *myCoords=DataArrayDouble::New();
196 myCoords->alloc(5,2);
197 std::copy(coordsS,coordsS+10,myCoords->getPointer());
198 meshS->setCoords(myCoords);
200 int connS[7]={0,3,4,1, 1,4,2};
201 meshS->allocateCells(2);
202 meshS->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,connS);
203 meshS->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connS+4);
204 meshS->finishInsertingCells();
205 ComponentTopology comptopo;
206 parameshS=new ParaMESH(meshS, *grp,"source mesh");
207 parafieldS=new ParaFIELD(ON_CELLS,NO_TIME,parameshS,comptopo);
208 parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
209 double *valsS=parafieldS->getField()->getArray()->getPointer();
210 valsS[0]=7.; valsS[1]=8.;
212 meshT=MEDCouplingUMesh::New();
213 meshT->setMeshDimension(2);
214 myCoords=DataArrayDouble::New();
215 myCoords->alloc(3,2);
216 std::copy(coordsT,coordsT+6,myCoords->getPointer());
217 meshT->setCoords(myCoords);
219 int connT[3]={0,2,1};
220 meshT->allocateCells(1);
221 meshT->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT);
222 meshT->finishInsertingCells();
223 parameshT=new ParaMESH(meshT,*grp,"target mesh");
224 parafieldT=new ParaFIELD(ON_CELLS,NO_TIME,parameshT,comptopo);
225 parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
226 double *valsT=parafieldT->getField()->getArray()->getPointer();
232 const double coordsS[10]={1.,0.,0.5,0.5,1.,0.5,0.5,1.,1.,1.};
233 const double coordsT[6]={0.,0.,0.5,0.5,0.,1.};
234 meshS=MEDCouplingUMesh::New();
235 meshS->setMeshDimension(2);
236 DataArrayDouble *myCoords=DataArrayDouble::New();
237 myCoords->alloc(5,2);
238 std::copy(coordsS,coordsS+10,myCoords->getPointer());
239 meshS->setCoords(myCoords);
241 int connS[7]={0,1,2, 1,3,4,2};
242 meshS->allocateCells(2);
243 meshS->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connS);
244 meshS->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,connS+3);
245 meshS->finishInsertingCells();
246 ComponentTopology comptopo;
247 parameshS=new ParaMESH(meshS,*grp,"source mesh");
248 parafieldS=new ParaFIELD(ON_CELLS,NO_TIME,parameshS,comptopo);
249 parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
250 double *valsS=parafieldS->getField()->getArray()->getPointer();
254 meshT=MEDCouplingUMesh::New();
255 meshT->setMeshDimension(2);
256 myCoords=DataArrayDouble::New();
257 myCoords->alloc(3,2);
258 std::copy(coordsT,coordsT+6,myCoords->getPointer());
259 meshT->setCoords(myCoords);
261 int connT[3]={0,2,1};
262 meshT->allocateCells(1);
263 meshT->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT);
264 meshT->finishInsertingCells();
265 parameshT=new ParaMESH(meshT,*grp,"target mesh");
266 parafieldT=new ParaFIELD(ON_CELLS,NO_TIME,parameshT,comptopo);
267 parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
268 double *valsT=parafieldT->getField()->getArray()->getPointer();
274 const double coordsS[8]={0.,0.5, 0.5,0.5, 0.,1., 0.5,1.};
275 const double coordsT[6]={0.5,0.5,0.,1.,1.,1.};
276 meshS=MEDCouplingUMesh::New();
277 meshS->setMeshDimension(2);
278 DataArrayDouble *myCoords=DataArrayDouble::New();
279 myCoords->alloc(4,2);
280 std::copy(coordsS,coordsS+8,myCoords->getPointer());
281 meshS->setCoords(myCoords);
283 int connS[4]={0,2,3,1};
284 meshS->allocateCells(1);
285 meshS->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,connS);
286 meshS->finishInsertingCells();
287 ComponentTopology comptopo;
288 parameshS=new ParaMESH(meshS,*grp,"source mesh");
289 parafieldS=new ParaFIELD(ON_CELLS,NO_TIME,parameshS,comptopo);
290 parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
291 double *valsS=parafieldS->getField()->getArray()->getPointer();
294 meshT=MEDCouplingUMesh::New();
295 meshT->setMeshDimension(2);
296 myCoords=DataArrayDouble::New();
297 myCoords->alloc(3,2);
298 std::copy(coordsT,coordsT+6,myCoords->getPointer());
299 meshT->setCoords(myCoords);
301 int connT[3]={0,1,2};
302 meshT->allocateCells(1);
303 meshT->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT);
304 meshT->finishInsertingCells();
305 parameshT=new ParaMESH(meshT,*grp,"target mesh");
306 parafieldT=new ParaFIELD(ON_CELLS,NO_TIME,parameshT,comptopo);
307 parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
308 double *valsT=parafieldT->getField()->getArray()->getPointer();
314 /*! Test case from the official doc of the OverlapDEC.
315 * WARNING: bounding boxes are tweaked here to make the case more interesting (i.e. to avoid an all to all exchange
316 * between all procs).
318 void ParaMEDMEMTest::testOverlapDEC1()
321 MPI_Comm_size(MPI_COMM_WORLD,&size);
322 MPI_Comm_rank(MPI_COMM_WORLD,&rank);
323 // char hostname[256];
324 // printf("(%d) PID %d on localhost ready for attach\n", rank, getpid());
327 if (size != 3) return ;
332 for (int i=0; i<nproc; i++)
335 CommInterface interface;
336 OverlapDEC dec(procs);
337 ProcessorGroup * grp = dec.getGroup();
338 MEDCouplingUMesh* meshS=0, *meshT=0;
339 ParaMESH* parameshS=0, *parameshT=0;
340 ParaFIELD* parafieldS=0, *parafieldT=0;
342 MPI_Barrier(MPI_COMM_WORLD);
343 prepareData1(rank, grp, meshS, meshT, parameshS, parameshT, parafieldS, parafieldT);
345 /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
346 * HACK ON BOUNDING BOX TO MAKE THIS CASE SIMPLE AND USABLE IN DEBUG
347 * Bounding boxes are slightly smaller than should be thus localising the work to be done
348 * and avoiding every proc talking to everyone else.
349 * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
351 dec.setBoundingBoxAdjustmentAbs(-1.0e-12);
353 dec.attachSourceLocalField(parafieldS);
354 dec.attachTargetLocalField(parafieldT);
356 dec.sendRecvData(true);
360 CPPUNIT_ASSERT_DOUBLES_EQUAL(8.75,parafieldT->getField()->getArray()->getIJ(0,0),1e-12);
364 CPPUNIT_ASSERT_DOUBLES_EQUAL(8.5,parafieldT->getField()->getArray()->getIJ(0,0),1e-12);
368 CPPUNIT_ASSERT_DOUBLES_EQUAL(10.5,parafieldT->getField()->getArray()->getIJ(0,0),1e-12);
377 MPI_Barrier(MPI_COMM_WORLD);
381 * Same as testOverlapDEC1() but with regular bounding boxes. If you're looking for a nice debug case,
382 * testOverlapDEC1() is more appropriate.
384 void ParaMEDMEMTest::testOverlapDEC2()
387 MPI_Comm_size(MPI_COMM_WORLD,&size);
388 MPI_Comm_rank(MPI_COMM_WORLD,&rank);
390 if (size != 3) return ;
395 for (int i=0; i<nproc; i++)
398 CommInterface interface;
399 OverlapDEC dec(procs);
400 ProcessorGroup * grp = dec.getGroup();
401 MEDCouplingUMesh* meshS=0, *meshT=0;
402 ParaMESH* parameshS=0, *parameshT=0;
403 ParaFIELD* parafieldS=0, *parafieldT=0;
405 MPI_Barrier(MPI_COMM_WORLD);
406 prepareData1(rank, grp, meshS, meshT, parameshS, parameshT, parafieldS, parafieldT);
408 /* Normal bounding boxes here: */
409 dec.setBoundingBoxAdjustmentAbs(+1.0e-12);
411 dec.attachSourceLocalField(parafieldS);
412 dec.attachTargetLocalField(parafieldT);
414 dec.sendRecvData(true);
418 CPPUNIT_ASSERT_DOUBLES_EQUAL(8.75,parafieldT->getField()->getArray()->getIJ(0,0),1e-12);
422 CPPUNIT_ASSERT_DOUBLES_EQUAL(8.5,parafieldT->getField()->getArray()->getIJ(0,0),1e-12);
426 CPPUNIT_ASSERT_DOUBLES_EQUAL(10.5,parafieldT->getField()->getArray()->getIJ(0,0),1e-12);
435 MPI_Barrier(MPI_COMM_WORLD);