1from petsc4py import PETSc 2import unittest 3import numpy as np 4 5# -------------------------------------------------------------------- 6 7ERR_SUP = 56 8 9class BaseTestPlex(object): 10 11 COMM = PETSc.COMM_WORLD 12 DIM = 1 13 CELLS = [[0, 1], [1, 2]] 14 COORDS = [[0.], [0.5], [1.]] 15 COMP = 1 16 DOFS = [1, 0] 17 18 def setUp(self): 19 self.plex = PETSc.DMPlex().createFromCellList(self.DIM, 20 self.CELLS, 21 self.COORDS, 22 comm=self.COMM) 23 24 def tearDown(self): 25 self.plex.destroy() 26 self.plex = None 27 28 def testTopology(self): 29 dim = self.plex.getDimension() 30 pStart, pEnd = self.plex.getChart() 31 cStart, cEnd = self.plex.getHeightStratum(0) 32 vStart, vEnd = self.plex.getDepthStratum(0) 33 numDepths = self.plex.getLabelSize("depth") 34 coords_raw = self.plex.getCoordinates().getArray() 35 coords = np.reshape(coords_raw, (vEnd - vStart, dim)) 36 self.assertEqual(dim, self.DIM) 37 self.assertEqual(numDepths, self.DIM+1) 38 if self.CELLS is not None: 39 self.assertEqual(cEnd-cStart, len(self.CELLS)) 40 if self.COORDS is not None: 41 self.assertEqual(vEnd-vStart, len(self.COORDS)) 42 self.assertTrue((coords == self.COORDS).all()) 43 44 def testClosure(self): 45 pStart, pEnd = self.plex.getChart() 46 for p in range(pStart, pEnd): 47 closure = self.plex.getTransitiveClosure(p)[0] 48 for c in closure: 49 cone = self.plex.getCone(c) 50 self.assertEqual(self.plex.getConeSize(c), len(cone)) 51 for i in cone: 52 self.assertIn(i, closure) 53 star = self.plex.getTransitiveClosure(p, useCone=False)[0] 54 for s in star: 55 support = self.plex.getSupport(s) 56 self.assertEqual(self.plex.getSupportSize(s), len(support)) 57 for i in support: 58 self.assertIn(i, star) 59 60 def testAdjacency(self): 61 PETSc.DMPlex.setAdjacencyUseAnchors(self.plex, False) 62 flag = PETSc.DMPlex.getAdjacencyUseAnchors(self.plex) 63 self.assertFalse(flag) 64 PETSc.DMPlex.setAdjacencyUseAnchors(self.plex, True) 65 flag = PETSc.DMPlex.getAdjacencyUseAnchors(self.plex) 66 self.assertTrue(flag) 67 PETSc.DMPlex.setBasicAdjacency(self.plex, False, False) 68 flagA, flagB = PETSc.DMPlex.getBasicAdjacency(self.plex) 69 self.assertFalse(flagA) 70 self.assertFalse(flagB) 71 PETSc.DMPlex.setBasicAdjacency(self.plex, True, True) 72 flagA, flagB = PETSc.DMPlex.getBasicAdjacency(self.plex) 73 self.assertTrue(flagA) 74 self.assertTrue(flagB) 75 pStart, pEnd = self.plex.getChart() 76 for p in range(pStart, pEnd): 77 adjacency = self.plex.getAdjacency(p) 78 self.assertTrue(p in adjacency) 79 self.assertTrue(len(adjacency) > 1) 80 81 def testSectionDofs(self): 82 self.plex.setNumFields(1) 83 section = self.plex.createSection([self.COMP], [self.DOFS]) 84 size = section.getStorageSize() 85 entity_dofs = [self.plex.getStratumSize("depth", d) * 86 self.DOFS[d] for d in range(self.DIM+1)] 87 self.assertEqual(sum(entity_dofs), size) 88 89 def testSectionClosure(self): 90 section = self.plex.createSection([self.COMP], [self.DOFS]) 91 self.plex.setSection(section) 92 vec = self.plex.createLocalVec() 93 pStart, pEnd = self.plex.getChart() 94 for p in range(pStart, pEnd): 95 for i in range(section.getDof(p)): 96 off = section.getOffset(p) 97 vec.setValue(off+i, p) 98 99 for p in range(pStart, pEnd): 100 point_closure = self.plex.getTransitiveClosure(p)[0] 101 dof_closure = self.plex.vecGetClosure(section, vec, p) 102 for p in dof_closure: 103 self.assertIn(p, point_closure) 104 105 def testBoundaryLabel(self): 106 pStart, pEnd = self.plex.getChart() 107 if (pEnd - pStart == 0): return 108 109 self.assertFalse(self.plex.hasLabel("boundary")) 110 self.plex.markBoundaryFaces("boundary") 111 self.assertTrue(self.plex.hasLabel("boundary")) 112 113 faces = self.plex.getStratumIS("boundary", 1) 114 for f in faces.getIndices(): 115 points, orient = self.plex.getTransitiveClosure(f, useCone=True) 116 for p in points: 117 self.plex.setLabelValue("boundary", p, 1) 118 119 for p in range(pStart, pEnd): 120 if self.plex.getLabelValue("boundary", p) != 1: 121 self.plex.setLabelValue("boundary", p, 2) 122 123 numBoundary = self.plex.getStratumSize("boundary", 1) 124 numInterior = self.plex.getStratumSize("boundary", 2) 125 self.assertNotEqual(numBoundary, pEnd - pStart) 126 self.assertNotEqual(numInterior, pEnd - pStart) 127 self.assertEqual(numBoundary + numInterior, pEnd - pStart) 128 129 130 def testAdapt(self): 131 dim = self.plex.getDimension() 132 if dim == 1: return 133 vStart, vEnd = self.plex.getDepthStratum(0) 134 numVertices = vEnd-vStart 135 metric_array = np.zeros([numVertices,dim,dim]) 136 for met in metric_array: 137 met[:,:] = np.diag([9]*dim) 138 metric = PETSc.Vec().createWithArray(metric_array) 139 try: 140 newplex = self.plex.adaptMetric(metric,"") 141 except PETSc.Error as exc: 142 if exc.ierr != ERR_SUP: raise 143 144 145# -------------------------------------------------------------------- 146 147class BaseTestPlex_2D(BaseTestPlex): 148 DIM = 2 149 CELLS = [[0, 1, 3], [1, 3, 4], [1, 2, 4], [2, 4, 5], 150 [3, 4, 6], [4, 6, 7], [4, 5, 7], [5, 7, 8]] 151 COORDS = [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0], 152 [0.0, 0.5], [0.5, 0.5], [1.0, 0.5], 153 [0.0, 1.0], [0.5, 1.0], [1.0, 1.0]] 154 DOFS = [1, 0, 0] 155 156class BaseTestPlex_3D(BaseTestPlex): 157 DIM = 3 158 CELLS = [[0, 2, 3, 7], [0, 2, 6, 7], [0, 4, 6, 7], 159 [0, 1, 3, 7], [0, 1, 5, 7], [0, 4, 5, 7]] 160 COORDS = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [1., 1., 0.], 161 [0., 0., 1.], [1., 0., 1.], [0., 1., 1.], [1., 1., 1.]] 162 DOFS = [1, 0, 0, 0] 163 164# -------------------------------------------------------------------- 165 166class TestPlex_1D(BaseTestPlex, unittest.TestCase): 167 pass 168 169class TestPlex_2D(BaseTestPlex_2D, unittest.TestCase): 170 pass 171 172class TestPlex_3D(BaseTestPlex_3D, unittest.TestCase): 173 pass 174 175class TestPlex_2D_P3(BaseTestPlex_2D, unittest.TestCase): 176 DOFS = [1, 2, 1] 177 178class TestPlex_3D_P3(BaseTestPlex_3D, unittest.TestCase): 179 DOFS = [1, 2, 1, 0] 180 181class TestPlex_3D_P4(BaseTestPlex_3D, unittest.TestCase): 182 DOFS = [1, 3, 3, 1] 183 184class TestPlex_2D_BoxTensor(BaseTestPlex_2D, unittest.TestCase): 185 CELLS = None 186 COORDS = None 187 def setUp(self): 188 self.plex = PETSc.DMPlex().createBoxMesh([3,3], simplex=False) 189 190class TestPlex_3D_BoxTensor(BaseTestPlex_3D, unittest.TestCase): 191 CELLS = None 192 COORDS = None 193 def setUp(self): 194 self.plex = PETSc.DMPlex().createBoxMesh([3,3,3], simplex=False) 195 196import sys 197try: 198 raise PETSc.Error 199 PETSc.DMPlex().createBoxMesh([2,2], simplex=True, comm=PETSc.COMM_SELF).destroy() 200except PETSc.Error: 201 pass 202else: 203 class TestPlex_2D_Box(BaseTestPlex_2D, unittest.TestCase): 204 CELLS = None 205 COORDS = None 206 def setUp(self): 207 self.plex = PETSc.DMPlex().createBoxMesh([1,1], simplex=True) 208 209 class TestPlex_2D_Boundary(BaseTestPlex_2D, unittest.TestCase): 210 CELLS = None 211 COORDS = None 212 def setUp(self): 213 boundary = PETSc.DMPlex().create(self.COMM) 214 boundary.createSquareBoundary([0., 0.], [1., 1.], [2, 2]) 215 boundary.setDimension(self.DIM-1) 216 self.plex = PETSc.DMPlex().generate(boundary) 217 218 class TestPlex_3D_Box(BaseTestPlex_3D, unittest.TestCase): 219 CELLS = None 220 COORDS = None 221 def setUp(self): 222 self.plex = PETSc.DMPlex().createBoxMesh([1,1,1], simplex=True) 223 224 class TestPlex_3D_Boundary(BaseTestPlex_3D, unittest.TestCase): 225 CELLS = None 226 COORDS = None 227 def setUp(self): 228 boundary = PETSc.DMPlex().create(self.COMM) 229 boundary.createCubeBoundary([0., 0., 0.], [1., 1., 1.], [1, 1, 1]) 230 boundary.setDimension(self.DIM-1) 231 self.plex = PETSc.DMPlex().generate(boundary) 232 233# -------------------------------------------------------------------- 234 235if __name__ == '__main__': 236 unittest.main() 237