2 # -*- coding: iso-8859-1 -*-
3 # Copyright (C) 2007-2023 CEA, EDF
5 # This library is free software; you can redistribute it and/or
6 # modify it under the terms of the GNU Lesser General Public
7 # License as published by the Free Software Foundation; either
8 # version 2.1 of the License, or (at your option) any later version.
10 # This library is distributed in the hope that it will be useful,
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 # Lesser General Public License for more details.
15 # You should have received a copy of the GNU Lesser General Public
16 # License along with this library; if not, write to the Free Software
17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 from medcoupling import *
23 from ParaMEDMEMTestTools import WriteInTmpDir
27 from mpi4py import MPI
30 class ParaMEDMEM_IK_DEC_Tests(unittest.TestCase):
31 """ See test_StructuredCoincidentDEC_py_1() for a quick start.
33 def generateFullSource(self):
34 """ The complete source mesh: 4 squares each divided in 2 diagonaly (so 8 cells in total) """
35 msh = self.generateFullTarget()
37 msh.setName("src_mesh")
38 fld = MEDCouplingFieldDouble(ON_CELLS, ONE_TIME)
39 fld.setMesh(msh); fld.setName("source_F");
40 da = DataArrayDouble(msh.getNumberOfCells())
46 def generateFullTarget(self):
47 """ The complete target mesh: 4 squares """
48 m1 = MEDCouplingCMesh("tgt_msh")
49 da = DataArrayDouble([0,1,2])
51 msh = m1.buildUnstructured()
55 # Below, the two functions emulating the set up of a piece of the source and target mesh
56 # on each proc. Obviously in real world problems, this comes from your code and is certainly
57 # not computed by cuting again from scratch the full-size mesh!!
59 def getPartialSource(self, rank):
60 """ Will return an empty mesh piece for rank=2 and 3 """
61 msh, f = self.generateFullSource()
63 sub_m, sub_f = msh[0:4], f[0:4]
65 sub_m, sub_f = msh[4:8], f[4:8]
69 def getPartialTarget(self, rank):
70 """ One square for each rank """
71 msh = self.generateFullTarget()
77 # Receiving side must prepare an empty field that will be filled by DEC:
78 fld = MEDCouplingFieldDouble(ON_CELLS, ONE_TIME)
79 da = DataArrayDouble(sub_m.getNumberOfCells())
85 def testInterpKernelDEC_ctor(self):
86 """ Test the various Python ctors """
87 size = MPI.COMM_WORLD.size
89 print("Should be run on 4 procs!")
91 # Define two processor groups
93 l1, l2 = range(nproc_source), range(size - nproc_source, size)
95 i1 = InterpKernelDEC.New(l1, l2)
96 # Should also work directly:
97 i2 = InterpKernelDEC(l1, l2)
99 interface = CommInterface()
100 source_group = MPIProcessorGroup(interface, list(l1))
101 target_group = MPIProcessorGroup(interface, list(l2))
102 i3 = InterpKernelDEC.New(source_group, target_group)
103 # Should also work directly:
104 i4 = InterpKernelDEC(source_group, target_group)
105 # With 2 iterables and a custom comm:
106 i5 = InterpKernelDEC.New(l1, l2, MPI.COMM_WORLD)
107 # Work directly with the **hack**
108 i6 = InterpKernelDEC(l1, l2, MPI._addressof(MPI.COMM_WORLD))
109 # Should fail with 2 proc groups **and** a communicator
110 self.assertRaises(InterpKernelException, InterpKernelDEC.New, source_group, target_group, MPI.COMM_WORLD)
111 self.assertRaises(Exception, InterpKernelDEC, source_group, target_group, MPI.COMM_WORLD)
112 i6.release(); i5.release(); i4.release(); i3.release(); i2.release(); i1.release()
113 source_group.release()
114 target_group.release()
117 def testInterpKernelDEC_2D_py_1(self):
118 """ This test illustrates a basic use of the InterpKernelDEC.
119 Look at the C++ documentation of the class for more informations.
121 size = MPI.COMM_WORLD.size
122 rank = MPI.COMM_WORLD.rank
124 print("Should be run on 4 procs!")
127 # Define two processor groups
129 procs_source = list(range(nproc_source))
130 procs_target = list(range(size - nproc_source, size))
132 interface = CommInterface()
133 source_group = MPIProcessorGroup(interface, procs_source)
134 target_group = MPIProcessorGroup(interface, procs_target)
135 idec = InterpKernelDEC(source_group, target_group)
137 # Write out full size meshes/fields for inspection
139 _, fld = self.generateFullSource()
140 mshT = self.generateFullTarget()
141 WriteField("./source_field_FULL.med", fld, True)
142 WriteUMesh("./target_mesh_FULL.med", mshT, True)
144 MPI.COMM_WORLD.Barrier() # really necessary??
147 # OK, let's go DEC !!
149 if source_group.containsMyRank():
150 _, fieldS = self.getPartialSource(rank)
151 fieldS.setNature(IntensiveMaximum) # The only policy supported for now ...
152 WriteField("./source_field_part_%d.med" % rank, fieldS, True)
153 idec.attachLocalField(fieldS)
157 if target_group.containsMyRank():
158 mshT, fieldT = self.getPartialTarget(rank)
159 fieldT.setNature(IntensiveMaximum)
160 WriteUMesh("./target_mesh_part_%d.med" % rank, mshT, True)
161 idec.attachLocalField(fieldT)
164 # Now the actual checks:
166 self.assertEqual(fieldT.getArray().getValues(), [1.0, 9.0])
168 self.assertEqual(fieldT.getArray().getValues(), [5.0, 13.0])
170 # Release DEC (this involves MPI exchanges -- notably the release of the communicator -- so better be done before MPI.Finalize()
172 source_group.release()
173 target_group.release()
174 MPI.COMM_WORLD.Barrier()
177 def test_InterpKernelDEC_2D_py_2(self):
178 """ More involved test using Para* objects.
180 size = MPI.COMM_WORLD.size
181 rank = MPI.COMM_WORLD.rank
183 print("Should be run on 5 procs!")
188 procs_source = list(range(nproc_source))
189 procs_target = list(range(size - nproc_source + 1, size))
191 interface = CommInterface()
192 target_group = MPIProcessorGroup(interface, procs_target)
193 source_group = MPIProcessorGroup(interface, procs_source)
194 dec = InterpKernelDEC(source_group, target_group)
196 data_dir = os.path.join(os.environ['MEDCOUPLING_ROOT_DIR'], "share", "resources", "med")
197 if not os.path.isdir(data_dir):
198 data_dir = os.environ.get('MED_RESOURCES_DIR',"::").split(":")[1]
200 filename_xml1 = os.path.join(data_dir, "square1_split")
201 filename_xml2 = os.path.join(data_dir, "square2_split")
203 MPI.COMM_WORLD.Barrier()
204 if source_group.containsMyRank():
205 filename = filename_xml1 + str(rank+1) + ".med"
206 meshname = "Mesh_2_" + str(rank+1)
207 mesh=ReadUMeshFromFile(filename,meshname,0)
208 paramesh=ParaMESH(mesh,source_group,"source mesh")
209 comptopo = ComponentTopology()
210 parafield = ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo)
211 parafield.getField().setNature(IntensiveMaximum)
212 nb_local=mesh.getNumberOfCells()
213 value = [1.0]*nb_local
214 parafield.getField().setValues(value)
215 icocofield = ICoCoMEDDoubleField(parafield.getField())
216 dec.attachLocalField(icocofield)
219 filename = filename_xml2 + str(rank - nproc_source + 1) + ".med"
220 meshname = "Mesh_3_" + str(rank - nproc_source + 1)
221 mesh=ReadUMeshFromFile(filename,meshname,0)
222 paramesh=ParaMESH(mesh,target_group,"target mesh")
223 comptopo = ComponentTopology()
224 parafield = ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo)
225 parafield.getField().setNature(IntensiveMaximum)
226 nb_local=mesh.getNumberOfCells()
227 value = [0.0]*nb_local
228 parafield.getField().setValues(value)
229 icocofield = ICoCoMEDDoubleField(parafield.getField())
230 dec.attachLocalField(icocofield)
233 if source_group.containsMyRank():
234 field_before_int = parafield.getVolumeIntegral(0,True)
236 dec.setForcedRenormalization(False)
239 field_after_int=parafield.getVolumeIntegral(0,True);
240 self.assertTrue(math.fabs(field_after_int-field_before_int)<1e-8)
244 dec.setForcedRenormalization(False)
250 # Some clean up that still needs MPI communication, so to be done before MPI_Finalize()
254 target_group.release()
255 source_group.release()
256 MPI.COMM_WORLD.Barrier()
258 def test_InterpKernelDEC_default(self):
260 [EDF27375] : Put a default value when non intersecting case
262 size = MPI.COMM_WORLD.size
263 rank = MPI.COMM_WORLD.rank
265 print("Should be run on 4 procs!")
268 procs_source = list(range(nproc_source))
269 procs_target = list(range(size - nproc_source, size))
271 interface = CommInterface()
272 target_group = MPIProcessorGroup(interface, procs_target)
273 source_group = MPIProcessorGroup(interface, procs_source)
274 dec = InterpKernelDEC(source_group, target_group)
276 MPI.COMM_WORLD.Barrier()
277 if source_group.containsMyRank():
278 mesh = eval("Source_Proc_{}".format(rank))()
279 nb_local=mesh.getNumberOfCells()
280 field = MEDCouplingFieldDouble(ON_CELLS)
281 field.setNature(IntensiveMaximum)
283 arr = DataArrayDouble(nb_local) ; arr.iota() ; arr += rank
285 dec.attachLocalField(field)
286 dec.synchronizeWithDefaultValue(-2000.0)
291 self.assertTrue(field.getArray().isEqual(DataArrayDouble([0.6,0.6,-2000]),1e-12))
292 self.assertTrue( dec.retrieveNonFetchedIds().isEqual(DataArrayInt([2])) )
294 self.assertTrue(field.getArray().isEqual(DataArrayDouble([1.0,-2000]),1e-12))
295 self.assertTrue( dec.retrieveNonFetchedIds().isEqual(DataArrayInt([1])) )
297 mesh = eval("Target_Proc_{}".format(rank))()
298 nb_local=mesh.getNumberOfCells()
299 field = MEDCouplingFieldDouble(ON_CELLS)
300 field.setNature(IntensiveMaximum)
302 arr = DataArrayDouble(nb_local) ; arr[:] = -20
304 dec.attachLocalField(field)
305 dec.synchronizeWithDefaultValue(-1000.0)
308 # matrix S0 / T2 = [[(0,S0,1),(1,S0,1.5)]
309 # IntensiveMaximum => [[(0,S0,1/2.5),(1,S0,1.5/2.5)]
311 self.assertTrue(field.getArray().isEqual(DataArrayDouble([0.6]),1e-12))
312 self.assertTrue( dec.retrieveNonFetchedIds().isEqual(DataArrayInt([])) )
314 # matrix S1 / T3 = [[],[(0,S1,1.0)],[(0,S1,2.0)],[]]
315 # IntensiveMaximum => [[],[(0,S1,1.0/1.0)],[(0,S1,2.0/2.0)],[]]
316 self.assertTrue(field.getArray().isEqual(DataArrayDouble([-1000.0, 1.0, 1.0, -1000.0]),1e-8))
317 self.assertTrue( dec.retrieveNonFetchedIds().isEqual(DataArrayInt([0,3])) )
321 # Some clean up that still needs MPI communication, so to be done before MPI_Finalize()
323 target_group.release()
324 source_group.release()
325 MPI.COMM_WORLD.Barrier()
328 coo = DataArrayDouble([(0,2),(2,2),(4,2),(0,4),(2,4),(4,4),(0,6),(2,6)])
329 m = MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
330 m.insertNextCell(NORM_QUAD4,[0,1,4,3])
331 m.insertNextCell(NORM_QUAD4,[1,2,5,4])
332 m.insertNextCell(NORM_QUAD4,[3,4,7,6])
336 coo = DataArrayDouble([(6,2),(8,2),(10,2),(6,4),(8,4),(10,4)])
337 m = MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
338 m.insertNextCell(NORM_QUAD4,[0,1,4,3])
339 m.insertNextCell(NORM_QUAD4,[1,2,5,4])
343 coo = DataArrayDouble([(1,0),(3.5,0),(1,3),(3.5,3)])
344 m = MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
345 m.insertNextCell(NORM_QUAD4,[0,1,3,2])
349 coo = DataArrayDouble([(6,0),(7,0),(8,0),(9,0),(10,0),
350 (6,1),(7,1),(9,1),(10,1),
353 m = MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
354 m.insertNextCell(NORM_QUAD4,[0,1,6,5])
355 m.insertNextCell(NORM_QUAD4,[1,2,10,9])
356 m.insertNextCell(NORM_QUAD4,[5,6,12,11])
357 m.insertNextCell(NORM_QUAD4,[3,4,8,7])
360 if __name__ == "__main__":