1import petsc4py 2from petsc4py import PETSc 3import unittest 4import os 5import filecmp 6import numpy as np 7 8# -------------------------------------------------------------------- 9 10ERR_ARG_OUTOFRANGE = 63 11 12class BaseTestPlex(object): 13 14 COMM = PETSc.COMM_WORLD 15 DIM = 1 16 CELLS = [[0, 1], [1, 2]] 17 COORDS = [[0.], [0.5], [1.]] 18 COMP = 1 19 DOFS = [1, 0] 20 21 def setUp(self): 22 self.plex = PETSc.DMPlex().createFromCellList(self.DIM, 23 self.CELLS, 24 self.COORDS, 25 comm=self.COMM) 26 27 def tearDown(self): 28 self.plex.destroy() 29 self.plex = None 30 31 def testTopology(self): 32 rank = self.COMM.rank 33 dim = self.plex.getDimension() 34 pStart, pEnd = self.plex.getChart() 35 cStart, cEnd = self.plex.getHeightStratum(0) 36 vStart, vEnd = self.plex.getDepthStratum(0) 37 numDepths = self.plex.getLabelSize("depth") 38 coords_raw = self.plex.getCoordinates().getArray() 39 coords = np.reshape(coords_raw, (vEnd - vStart, dim)) 40 self.assertEqual(dim, self.DIM) 41 self.assertEqual(numDepths, self.DIM+1) 42 if rank == 0 and self.CELLS is not None: 43 self.assertEqual(cEnd-cStart, len(self.CELLS)) 44 if rank == 0 and self.COORDS is not None: 45 self.assertEqual(vEnd-vStart, len(self.COORDS)) 46 self.assertTrue((coords == self.COORDS).all()) 47 48 def testClosure(self): 49 pStart, pEnd = self.plex.getChart() 50 for p in range(pStart, pEnd): 51 closure = self.plex.getTransitiveClosure(p)[0] 52 for c in closure: 53 cone = self.plex.getCone(c) 54 self.assertEqual(self.plex.getConeSize(c), len(cone)) 55 for i in cone: 56 self.assertIn(i, closure) 57 star = self.plex.getTransitiveClosure(p, useCone=False)[0] 58 for s in star: 59 support = self.plex.getSupport(s) 60 self.assertEqual(self.plex.getSupportSize(s), len(support)) 61 for i in support: 62 self.assertIn(i, star) 63 64 def testAdjacency(self): 65 PETSc.DMPlex.setAdjacencyUseAnchors(self.plex, False) 66 flag = PETSc.DMPlex.getAdjacencyUseAnchors(self.plex) 67 self.assertFalse(flag) 68 PETSc.DMPlex.setAdjacencyUseAnchors(self.plex, True) 69 flag = PETSc.DMPlex.getAdjacencyUseAnchors(self.plex) 70 self.assertTrue(flag) 71 PETSc.DMPlex.setBasicAdjacency(self.plex, False, False) 72 flagA, flagB = PETSc.DMPlex.getBasicAdjacency(self.plex) 73 self.assertFalse(flagA) 74 self.assertFalse(flagB) 75 PETSc.DMPlex.setBasicAdjacency(self.plex, True, True) 76 flagA, flagB = PETSc.DMPlex.getBasicAdjacency(self.plex) 77 self.assertTrue(flagA) 78 self.assertTrue(flagB) 79 pStart, pEnd = self.plex.getChart() 80 for p in range(pStart, pEnd): 81 adjacency = self.plex.getAdjacency(p) 82 self.assertTrue(p in adjacency) 83 self.assertTrue(len(adjacency) > 1) 84 85 def testSectionDofs(self): 86 self.plex.setNumFields(1) 87 section = self.plex.createSection([self.COMP], [self.DOFS]) 88 size = section.getStorageSize() 89 entity_dofs = [self.plex.getStratumSize("depth", d) * 90 self.DOFS[d] for d in range(self.DIM+1)] 91 self.assertEqual(sum(entity_dofs), size) 92 93 def testSectionClosure(self): 94 section = self.plex.createSection([self.COMP], [self.DOFS]) 95 self.plex.setSection(section) 96 vec = self.plex.createLocalVec() 97 pStart, pEnd = self.plex.getChart() 98 for p in range(pStart, pEnd): 99 for i in range(section.getDof(p)): 100 off = section.getOffset(p) 101 vec.setValue(off+i, p) 102 103 for p in range(pStart, pEnd): 104 point_closure = self.plex.getTransitiveClosure(p)[0] 105 dof_closure = self.plex.vecGetClosure(section, vec, p) 106 for p in dof_closure: 107 self.assertIn(p, point_closure) 108 109 def testBoundaryLabel(self): 110 pStart, pEnd = self.plex.getChart() 111 if (pEnd - pStart == 0): return 112 113 self.assertFalse(self.plex.hasLabel("boundary")) 114 self.plex.markBoundaryFaces("boundary") 115 self.assertTrue(self.plex.hasLabel("boundary")) 116 117 faces = self.plex.getStratumIS("boundary", 1) 118 for f in faces.getIndices(): 119 points, orient = self.plex.getTransitiveClosure(f, useCone=True) 120 for p in points: 121 self.plex.setLabelValue("boundary", p, 1) 122 123 for p in range(pStart, pEnd): 124 if self.plex.getLabelValue("boundary", p) != 1: 125 self.plex.setLabelValue("boundary", p, 2) 126 127 numBoundary = self.plex.getStratumSize("boundary", 1) 128 numInterior = self.plex.getStratumSize("boundary", 2) 129 self.assertNotEqual(numBoundary, pEnd - pStart) 130 self.assertNotEqual(numInterior, pEnd - pStart) 131 self.assertEqual(numBoundary + numInterior, pEnd - pStart) 132 133 def testMetric(self): 134 if self.DIM == 1: return 135 self.plex.distribute() 136 if self.CELLS is None and not self.plex.isSimplex(): return 137 self.plex.orient() 138 139 h_min = 1.0e-30 140 h_max = 1.0e+30 141 a_max = 1.0e+10 142 target = 8.0 143 p = 1.0 144 beta = 1.3 145 hausd = 0.01 146 self.plex.metricSetUniform(False) 147 self.plex.metricSetIsotropic(False) 148 self.plex.metricSetRestrictAnisotropyFirst(False) 149 self.plex.metricSetNoInsertion(False) 150 self.plex.metricSetNoSwapping(False) 151 self.plex.metricSetNoMovement(False) 152 self.plex.metricSetNoSurf(False) 153 self.plex.metricSetVerbosity(-1) 154 self.plex.metricSetNumIterations(3) 155 self.plex.metricSetMinimumMagnitude(h_min) 156 self.plex.metricSetMaximumMagnitude(h_max) 157 self.plex.metricSetMaximumAnisotropy(a_max) 158 self.plex.metricSetTargetComplexity(target) 159 self.plex.metricSetNormalizationOrder(p) 160 self.plex.metricSetGradationFactor(beta) 161 self.plex.metricSetHausdorffNumber(hausd) 162 163 self.assertFalse(self.plex.metricIsUniform()) 164 self.assertFalse(self.plex.metricIsIsotropic()) 165 self.assertFalse(self.plex.metricRestrictAnisotropyFirst()) 166 self.assertFalse(self.plex.metricNoInsertion()) 167 self.assertFalse(self.plex.metricNoSwapping()) 168 self.assertFalse(self.plex.metricNoMovement()) 169 self.assertFalse(self.plex.metricNoSurf()) 170 assert self.plex.metricGetVerbosity() == -1 171 assert self.plex.metricGetNumIterations() == 3 172 assert np.isclose(self.plex.metricGetMinimumMagnitude(), h_min) 173 assert np.isclose(self.plex.metricGetMaximumMagnitude(), h_max) 174 assert np.isclose(self.plex.metricGetMaximumAnisotropy(), a_max) 175 assert np.isclose(self.plex.metricGetTargetComplexity(), target) 176 assert np.isclose(self.plex.metricGetNormalizationOrder(), p) 177 assert np.isclose(self.plex.metricGetGradationFactor(), beta) 178 assert np.isclose(self.plex.metricGetHausdorffNumber(), hausd) 179 180 metric1 = self.plex.metricCreateUniform(0.5) 181 metric2 = self.plex.metricCreateUniform(1.0) 182 metric = self.plex.metricCreate() 183 det = self.plex.metricDeterminantCreate() 184 self.plex.metricAverage2(metric1, metric2, metric) 185 metric1.array[:] *= 1.5 186 assert np.allclose(metric.array, metric1.array) 187 self.plex.metricIntersection2(metric1, metric2, metric) 188 assert np.allclose(metric.array, metric2.array) 189 self.plex.metricEnforceSPD(metric, metric1, det) 190 assert np.allclose(metric.array, metric1.array) 191 self.plex.metricNormalize(metric, metric1, det, restrictSizes=False, restrictAnisotropy=False) 192 metric2.scale(pow(target, 2.0/self.DIM)) 193 assert np.allclose(metric1.array, metric2.array) 194 195 def testAdapt(self): 196 if self.DIM == 1: return 197 self.plex.orient() 198 plex = self.plex.refine() 199 plex.distribute() 200 if self.CELLS is None and not plex.isSimplex(): return 201 if sum(self.DOFS) > 1: return 202 metric = plex.metricCreateUniform(9.0) 203 try: 204 newplex = plex.adaptMetric(metric,"") 205 except PETSc.Error as exc: 206 if exc.ierr != ERR_ARG_OUTOFRANGE: raise 207 208 209# -------------------------------------------------------------------- 210 211class BaseTestPlex_2D(BaseTestPlex): 212 DIM = 2 213 CELLS = [[0, 1, 3], [1, 3, 4], [1, 2, 4], [2, 4, 5], 214 [3, 4, 6], [4, 6, 7], [4, 5, 7], [5, 7, 8]] 215 COORDS = [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0], 216 [0.0, 0.5], [0.5, 0.5], [1.0, 0.5], 217 [0.0, 1.0], [0.5, 1.0], [1.0, 1.0]] 218 DOFS = [1, 0, 0] 219 220class BaseTestPlex_3D(BaseTestPlex): 221 DIM = 3 222 CELLS = [[0, 2, 3, 7], [0, 2, 6, 7], [0, 4, 6, 7], 223 [0, 1, 3, 7], [0, 1, 5, 7], [0, 4, 5, 7]] 224 COORDS = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [1., 1., 0.], 225 [0., 0., 1.], [1., 0., 1.], [0., 1., 1.], [1., 1., 1.]] 226 DOFS = [1, 0, 0, 0] 227 228# -------------------------------------------------------------------- 229 230class TestPlex_1D(BaseTestPlex, unittest.TestCase): 231 pass 232 233class TestPlex_2D(BaseTestPlex_2D, unittest.TestCase): 234 pass 235 236class TestPlex_3D(BaseTestPlex_3D, unittest.TestCase): 237 pass 238 239class TestPlex_2D_P3(BaseTestPlex_2D, unittest.TestCase): 240 DOFS = [1, 2, 1] 241 242class TestPlex_3D_P3(BaseTestPlex_3D, unittest.TestCase): 243 DOFS = [1, 2, 1, 0] 244 245class TestPlex_3D_P4(BaseTestPlex_3D, unittest.TestCase): 246 DOFS = [1, 3, 3, 1] 247 248class TestPlex_2D_BoxTensor(BaseTestPlex_2D, unittest.TestCase): 249 CELLS = None 250 COORDS = None 251 def setUp(self): 252 self.plex = PETSc.DMPlex().createBoxMesh([3,3], simplex=False) 253 254class TestPlex_3D_BoxTensor(BaseTestPlex_3D, unittest.TestCase): 255 CELLS = None 256 COORDS = None 257 def setUp(self): 258 self.plex = PETSc.DMPlex().createBoxMesh([3,3,3], simplex=False) 259 260try: 261 raise PETSc.Error 262 PETSc.DMPlex().createBoxMesh([2,2], simplex=True, comm=PETSc.COMM_SELF).destroy() 263except PETSc.Error: 264 pass 265else: 266 class TestPlex_2D_Box(BaseTestPlex_2D, unittest.TestCase): 267 CELLS = None 268 COORDS = None 269 def setUp(self): 270 self.plex = PETSc.DMPlex().createBoxMesh([1,1], simplex=True) 271 272 class TestPlex_2D_Boundary(BaseTestPlex_2D, unittest.TestCase): 273 CELLS = None 274 COORDS = None 275 def setUp(self): 276 boundary = PETSc.DMPlex().create(self.COMM) 277 boundary.createSquareBoundary([0., 0.], [1., 1.], [2, 2]) 278 boundary.setDimension(self.DIM-1) 279 self.plex = PETSc.DMPlex().generate(boundary) 280 281 class TestPlex_3D_Box(BaseTestPlex_3D, unittest.TestCase): 282 CELLS = None 283 COORDS = None 284 def setUp(self): 285 self.plex = PETSc.DMPlex().createBoxMesh([1,1,1], simplex=True) 286 287 class TestPlex_3D_Boundary(BaseTestPlex_3D, unittest.TestCase): 288 CELLS = None 289 COORDS = None 290 def setUp(self): 291 boundary = PETSc.DMPlex().create(self.COMM) 292 boundary.createCubeBoundary([0., 0., 0.], [1., 1., 1.], [1, 1, 1]) 293 boundary.setDimension(self.DIM-1) 294 self.plex = PETSc.DMPlex().generate(boundary) 295 296# -------------------------------------------------------------------- 297 298PETSC_DIR = petsc4py.get_config()['PETSC_DIR'] 299 300def check_dtype(method): 301 def wrapper(self, *args, **kwargs): 302 if PETSc.ScalarType is PETSc.ComplexType: 303 return 304 else: 305 return method(self, *args, **kwargs) 306 return wrapper 307 308def check_package(method): 309 def wrapper(self, *args, **kwargs): 310 if not PETSc.Sys.hasExternalPackage("hdf5"): 311 return 312 elif self.PARTITIONERTYPE != "simple" and \ 313 not PETSc.Sys.hasExternalPackage(self.PARTITIONERTYPE): 314 return 315 else: 316 return method(self, *args, **kwargs) 317 return wrapper 318 319def check_nsize(method): 320 def wrapper(self, *args, **kwargs): 321 if PETSc.COMM_WORLD.size != self.NSIZE: 322 return 323 else: 324 return method(self, *args, **kwargs) 325 return wrapper 326 327class BaseTestPlexHDF5(object): 328 NSIZE = 4 329 NTIMES = 3 330 331 def setUp(self): 332 self.txtvwr = PETSc.Viewer() 333 334 def tearDown(self): 335 if not PETSc.COMM_WORLD.rank: 336 if os.path.exists(self.outfile()): 337 os.remove(self.outfile()) 338 if os.path.exists(self.tmp_output_file()): 339 os.remove(self.tmp_output_file()) 340 self.txtvwr = None 341 342 def _name(self): 343 return "%s_outformat-%s_%s" % (self.SUFFIX, 344 self.OUTFORMAT, 345 self.PARTITIONERTYPE) 346 347 def infile(self): 348 return os.path.join(PETSC_DIR, "share/petsc/datafiles/", 349 "meshes/blockcylinder-50.h5") 350 351 def outfile(self): 352 return os.path.join("./temp_test_dmplex_%s.h5" % self._name()) 353 354 def informat(self): 355 return PETSc.Viewer.Format.HDF5_XDMF 356 357 def outformat(self): 358 d = {"hdf5_petsc": PETSc.Viewer.Format.HDF5_PETSC, 359 "hdf5_xdmf": PETSc.Viewer.Format.HDF5_XDMF} 360 return d[self.OUTFORMAT] 361 362 def partitionerType(self): 363 d = {"simple": PETSc.Partitioner.Type.SIMPLE, 364 "ptscotch": PETSc.Partitioner.Type.PTSCOTCH, 365 "parmetis": PETSc.Partitioner.Type.PARMETIS} 366 return d[self.PARTITIONERTYPE] 367 368 def ref_output_file(self): 369 return os.path.join(PETSC_DIR, "src/dm/impls/plex/tutorials/", 370 "output/ex5_%s.out" % self._name()) 371 372 def tmp_output_file(self): 373 return os.path.join("./temp_test_dmplex_%s.out" % self._name()) 374 375 def outputText(self, msg, comm): 376 if not comm.rank: 377 with open(self.tmp_output_file(), 'a') as f: 378 f.write(msg) 379 380 def outputPlex(self, plex): 381 self.txtvwr.createASCII(self.tmp_output_file(), 382 mode='a', comm=plex.comm) 383 plex.view(viewer=self.txtvwr) 384 self.txtvwr.destroy() 385 386 @check_dtype 387 @check_package 388 @check_nsize 389 def testViewLoadCycle(self): 390 grank = PETSc.COMM_WORLD.rank 391 for i in range(self.NTIMES): 392 if i == 0: 393 infname = self.infile() 394 informt = self.informat() 395 else: 396 infname = self.outfile() 397 informt = self.outformat() 398 if self.HETEROGENEOUS: 399 mycolor = (grank > self.NTIMES - i) 400 else: 401 mycolor = 0 402 try: 403 import mpi4py 404 except ImportError: 405 self.skipTest('mpi4py') # throws special exception to signal test skip 406 mpicomm = PETSc.COMM_WORLD.tompi4py() 407 comm = PETSc.Comm(comm=mpicomm.Split(color=mycolor, key=grank)) 408 if mycolor == 0: 409 self.outputText("Begin cycle %d\n" % i, comm) 410 plex = PETSc.DMPlex() 411 vwr = PETSc.ViewerHDF5() 412 # Create plex 413 plex.create(comm=comm) 414 plex.setName("DMPlex Object") 415 # Load data from XDMF into dm in parallel 416 vwr.create(infname, mode='r', comm=comm) 417 vwr.pushFormat(format=informt) 418 plex.load(viewer=vwr) 419 plex.setOptionsPrefix("loaded_") 420 plex.distributeSetDefault(False) 421 plex.setFromOptions() 422 vwr.popFormat() 423 vwr.destroy() 424 self.outputPlex(plex) 425 # Test DM is indeed distributed 426 flg = plex.isDistributed() 427 self.outputText("Loaded mesh distributed? %s\n" % 428 str(flg).upper(), comm) 429 # Interpolate 430 plex.interpolate() 431 plex.setOptionsPrefix("interpolated_") 432 plex.setFromOptions() 433 self.outputPlex(plex) 434 # Redistribute 435 part = plex.getPartitioner() 436 part.setType(self.partitionerType()) 437 _ = plex.distribute(overlap=0) 438 plex.setName("DMPlex Object") 439 plex.setOptionsPrefix("redistributed_") 440 plex.setFromOptions() 441 self.outputPlex(plex) 442 # Save redistributed dm to XDMF in parallel 443 vwr.create(self.outfile(), mode='w', comm=comm) 444 vwr.pushFormat(format=self.outformat()) 445 plex.setName("DMPlex Object") 446 plex.view(viewer=vwr) 447 vwr.popFormat() 448 vwr.destroy() 449 # Destroy plex 450 plex.destroy() 451 self.outputText("End cycle %d\n--------\n" % i, comm) 452 PETSc.COMM_WORLD.Barrier() 453 # Check that the output is identical to that of plex/tutorial/ex5.c. 454 self.assertTrue(filecmp.cmp(self.tmp_output_file(), 455 self.ref_output_file(), shallow=False), 456 'Contents of the files not the same.') 457 PETSc.COMM_WORLD.Barrier() 458 459class BaseTestPlexHDF5Homogeneous(BaseTestPlexHDF5): 460 """Test save on N / load on N.""" 461 SUFFIX = 0 462 HETEROGENEOUS = False 463 464class BaseTestPlexHDF5Heterogeneous(BaseTestPlexHDF5): 465 """Test save on N / load on M.""" 466 SUFFIX = 1 467 HETEROGENEOUS = True 468 469class TestPlexHDF5PETSCSimpleHomogeneous(BaseTestPlexHDF5Homogeneous, 470 unittest.TestCase): 471 OUTFORMAT = "hdf5_petsc" 472 PARTITIONERTYPE = "simple" 473 474""" 475Skipping. PTScotch produces different distributions when run 476in a sequence in a single session. 477 478class TestPlexHDF5PETSCPTScotchHomogeneous(BaseTestPlexHDF5Homogeneous, 479 unittest.TestCase): 480 OUTFORMAT = "hdf5_petsc" 481 PARTITIONERTYPE = "ptscotch" 482""" 483 484class TestPlexHDF5PETSCParmetisHomogeneous(BaseTestPlexHDF5Homogeneous, 485 unittest.TestCase): 486 OUTFORMAT = "hdf5_petsc" 487 PARTITIONERTYPE = "parmetis" 488 489class TestPlexHDF5XDMFSimpleHomogeneous(BaseTestPlexHDF5Homogeneous, 490 unittest.TestCase): 491 OUTFORMAT = "hdf5_xdmf" 492 PARTITIONERTYPE = "simple" 493 494""" 495Skipping. PTScotch produces different distributions when run 496in a sequence in a single session. 497 498class TestPlexHDF5XDMFPTScotchHomogeneous(BaseTestPlexHDF5Homogeneous, 499 unittest.TestCase): 500 OUTFORMAT = "hdf5_xdmf" 501 PARTITIONERTYPE = "ptscotch" 502""" 503 504class TestPlexHDF5XDMFParmetisHomogeneous(BaseTestPlexHDF5Homogeneous, 505 unittest.TestCase): 506 OUTFORMAT = "hdf5_xdmf" 507 PARTITIONERTYPE = "parmetis" 508 509class TestPlexHDF5PETSCSimpleHeterogeneous(BaseTestPlexHDF5Heterogeneous, 510 unittest.TestCase): 511 OUTFORMAT = "hdf5_petsc" 512 PARTITIONERTYPE = "simple" 513 514""" 515Skipping. PTScotch produces different distributions when run 516in a sequence in a single session. 517 518class TestPlexHDF5PETSCPTScotchHeterogeneous(BaseTestPlexHDF5Heterogeneous, 519 unittest.TestCase): 520 OUTFORMAT = "hdf5_petsc" 521 PARTITIONERTYPE = "ptscotch" 522""" 523 524class TestPlexHDF5PETSCParmetisHeterogeneous(BaseTestPlexHDF5Heterogeneous, 525 unittest.TestCase): 526 OUTFORMAT = "hdf5_petsc" 527 PARTITIONERTYPE = "parmetis" 528 529class TestPlexHDF5XDMFSimpleHeterogeneous(BaseTestPlexHDF5Heterogeneous, 530 unittest.TestCase): 531 OUTFORMAT = "hdf5_xdmf" 532 PARTITIONERTYPE = "simple" 533 534class TestPlexHDF5XDMFPTScotchHeterogeneous(BaseTestPlexHDF5Heterogeneous, 535 unittest.TestCase): 536 OUTFORMAT = "hdf5_xdmf" 537 PARTITIONERTYPE = "ptscotch" 538 539class TestPlexHDF5XDMFParmetisHeterogeneous(BaseTestPlexHDF5Heterogeneous, 540 unittest.TestCase): 541 OUTFORMAT = "hdf5_xdmf" 542 PARTITIONERTYPE = "parmetis" 543 544# -------------------------------------------------------------------- 545 546if __name__ == '__main__': 547 unittest.main() 548