1 static char help[] = "Tests for DMLabel\n\n"; 2 3 #include <petscdmplex.h> 4 #include <petsc/private/dmimpl.h> 5 6 static PetscErrorCode TestInsertion() 7 { 8 DMLabel label, label2; 9 const PetscInt values[5] = {0, 3, 4, -1, 176}, N = 10000; 10 PetscInt i, v; 11 12 PetscFunctionBegin; 13 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Test Label", &label)); 14 PetscCall(DMLabelSetDefaultValue(label, -100)); 15 for (i = 0; i < N; ++i) PetscCall(DMLabelSetValue(label, i, values[i % 5])); 16 /* Test get in hash mode */ 17 for (i = 0; i < N; ++i) { 18 PetscInt val; 19 20 PetscCall(DMLabelGetValue(label, i, &val)); 21 PetscCheck(val == values[i % 5], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Value %" PetscInt_FMT " for point %" PetscInt_FMT " should be %" PetscInt_FMT, val, i, values[i % 5]); 22 } 23 /* Test stratum */ 24 for (v = 0; v < 5; ++v) { 25 IS stratum; 26 const PetscInt *points; 27 PetscInt n; 28 29 PetscCall(DMLabelGetStratumIS(label, values[v], &stratum)); 30 PetscCheck(stratum, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Stratum %" PetscInt_FMT " is empty!", v); 31 PetscCall(ISGetIndices(stratum, &points)); 32 PetscCall(ISGetLocalSize(stratum, &n)); 33 for (i = 0; i < n; ++i) PetscCheck(points[i] == i * 5 + v, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " should be %" PetscInt_FMT, points[i], i * 5 + v); 34 PetscCall(ISRestoreIndices(stratum, &points)); 35 PetscCall(ISDestroy(&stratum)); 36 } 37 /* Test get in array mode */ 38 for (i = 0; i < N; ++i) { 39 PetscInt val; 40 41 PetscCall(DMLabelGetValue(label, i, &val)); 42 PetscCheck(val == values[i % 5], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Value %" PetscInt_FMT " should be %" PetscInt_FMT, val, values[i % 5]); 43 } 44 /* Test Duplicate */ 45 PetscCall(DMLabelDuplicate(label, &label2)); 46 for (i = 0; i < N; ++i) { 47 PetscInt val; 48 49 PetscCall(DMLabelGetValue(label2, i, &val)); 50 PetscCheck(val == values[i % 5], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Value %" PetscInt_FMT " should be %" PetscInt_FMT, val, values[i % 5]); 51 } 52 PetscCall(DMLabelDestroy(&label2)); 53 PetscCall(DMLabelDestroy(&label)); 54 PetscFunctionReturn(PETSC_SUCCESS); 55 } 56 57 static PetscErrorCode TestEmptyStrata(MPI_Comm comm) 58 { 59 DM dm, dmDist; 60 PetscPartitioner part; 61 PetscInt c0[6] = {2, 3, 6, 7, 9, 11}; 62 PetscInt c1[6] = {4, 5, 7, 8, 10, 12}; 63 PetscInt c2[4] = {13, 15, 19, 21}; 64 PetscInt c3[4] = {14, 16, 20, 22}; 65 PetscInt c4[4] = {15, 17, 21, 23}; 66 PetscInt c5[4] = {16, 18, 22, 24}; 67 PetscInt c6[4] = {13, 14, 19, 20}; 68 PetscInt c7[4] = {15, 16, 21, 22}; 69 PetscInt c8[4] = {17, 18, 23, 24}; 70 PetscInt c9[4] = {13, 14, 15, 16}; 71 PetscInt c10[4] = {15, 16, 17, 18}; 72 PetscInt c11[4] = {19, 20, 21, 22}; 73 PetscInt c12[4] = {21, 22, 23, 24}; 74 PetscInt dim = 3; 75 PetscMPIInt rank; 76 77 PetscFunctionBegin; 78 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 79 /* A 3D box with two adjacent cells, sharing one face and four vertices */ 80 PetscCall(DMCreate(comm, &dm)); 81 PetscCall(DMSetType(dm, DMPLEX)); 82 PetscCall(DMSetDimension(dm, dim)); 83 if (rank == 0) { 84 PetscCall(DMPlexSetChart(dm, 0, 25)); 85 PetscCall(DMPlexSetConeSize(dm, 0, 6)); 86 PetscCall(DMPlexSetConeSize(dm, 1, 6)); 87 PetscCall(DMPlexSetConeSize(dm, 2, 4)); 88 PetscCall(DMPlexSetConeSize(dm, 3, 4)); 89 PetscCall(DMPlexSetConeSize(dm, 4, 4)); 90 PetscCall(DMPlexSetConeSize(dm, 5, 4)); 91 PetscCall(DMPlexSetConeSize(dm, 6, 4)); 92 PetscCall(DMPlexSetConeSize(dm, 7, 4)); 93 PetscCall(DMPlexSetConeSize(dm, 8, 4)); 94 PetscCall(DMPlexSetConeSize(dm, 9, 4)); 95 PetscCall(DMPlexSetConeSize(dm, 10, 4)); 96 PetscCall(DMPlexSetConeSize(dm, 11, 4)); 97 PetscCall(DMPlexSetConeSize(dm, 12, 4)); 98 } 99 PetscCall(DMSetUp(dm)); 100 if (rank == 0) { 101 PetscCall(DMPlexSetCone(dm, 0, c0)); 102 PetscCall(DMPlexSetCone(dm, 1, c1)); 103 PetscCall(DMPlexSetCone(dm, 2, c2)); 104 PetscCall(DMPlexSetCone(dm, 3, c3)); 105 PetscCall(DMPlexSetCone(dm, 4, c4)); 106 PetscCall(DMPlexSetCone(dm, 5, c5)); 107 PetscCall(DMPlexSetCone(dm, 6, c6)); 108 PetscCall(DMPlexSetCone(dm, 7, c7)); 109 PetscCall(DMPlexSetCone(dm, 8, c8)); 110 PetscCall(DMPlexSetCone(dm, 9, c9)); 111 PetscCall(DMPlexSetCone(dm, 10, c10)); 112 PetscCall(DMPlexSetCone(dm, 11, c11)); 113 PetscCall(DMPlexSetCone(dm, 12, c12)); 114 } 115 PetscCall(DMPlexSymmetrize(dm)); 116 /* Create a user managed depth label, so that we can leave out edges */ 117 { 118 DMLabel label; 119 PetscInt numValues, maxValues = 0, v; 120 121 PetscCall(DMCreateLabel(dm, "depth")); 122 PetscCall(DMPlexGetDepthLabel(dm, &label)); 123 if (rank == 0) { 124 PetscInt i; 125 126 for (i = 0; i < 25; ++i) { 127 if (i < 2) PetscCall(DMLabelSetValue(label, i, 3)); 128 else if (i < 13) PetscCall(DMLabelSetValue(label, i, 2)); 129 else { 130 if (i == 13) PetscCall(DMLabelAddStratum(label, 1)); 131 PetscCall(DMLabelSetValue(label, i, 0)); 132 } 133 } 134 } 135 PetscCall(DMLabelGetNumValues(label, &numValues)); 136 PetscCallMPI(MPIU_Allreduce(&numValues, &maxValues, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm))); 137 for (v = numValues; v < maxValues; ++v) PetscCall(DMLabelAddStratum(label, v)); 138 } 139 { 140 DMLabel label; 141 PetscCall(DMPlexGetDepthLabel(dm, &label)); 142 PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_(comm))); 143 } 144 PetscCall(DMPlexGetPartitioner(dm, &part)); 145 PetscCall(PetscPartitionerSetFromOptions(part)); 146 PetscCall(DMPlexDistribute(dm, 1, NULL, &dmDist)); 147 if (dmDist) { 148 PetscCall(DMDestroy(&dm)); 149 dm = dmDist; 150 } 151 { 152 DMLabel label; 153 PetscCall(DMPlexGetDepthLabel(dm, &label)); 154 PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_(comm))); 155 } 156 /* Create a cell vector */ 157 { 158 Vec v; 159 PetscSection s; 160 PetscInt numComp[] = {1}; 161 PetscInt dof[] = {0, 0, 0, 1}; 162 PetscInt N; 163 164 PetscCall(DMSetNumFields(dm, 1)); 165 PetscCall(DMPlexCreateSection(dm, NULL, numComp, dof, 0, NULL, NULL, NULL, NULL, &s)); 166 PetscCall(DMSetLocalSection(dm, s)); 167 PetscCall(PetscSectionDestroy(&s)); 168 PetscCall(DMCreateGlobalVector(dm, &v)); 169 PetscCall(VecGetSize(v, &N)); 170 if (N != 2) { 171 PetscCall(DMView(dm, PETSC_VIEWER_STDOUT_(comm))); 172 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "FAIL: Vector size %" PetscInt_FMT " != 2", N); 173 } 174 PetscCall(VecDestroy(&v)); 175 } 176 PetscCall(DMDestroy(&dm)); 177 PetscFunctionReturn(PETSC_SUCCESS); 178 } 179 180 static PetscErrorCode TestDistribution(MPI_Comm comm) 181 { 182 DM dm, dmDist; 183 PetscPartitioner part; 184 DMLabel label; 185 char filename[PETSC_MAX_PATH_LEN]; 186 const char *name = "test label"; 187 PetscInt overlap = 0, cStart, cEnd, c; 188 PetscMPIInt rank; 189 PetscBool flg; 190 191 PetscFunctionBegin; 192 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 193 PetscCall(PetscOptionsGetString(NULL, NULL, "-filename", filename, sizeof(filename), &flg)); 194 if (!flg) PetscFunctionReturn(PETSC_SUCCESS); 195 PetscCall(PetscOptionsGetInt(NULL, NULL, "-overlap", &overlap, NULL)); 196 PetscCall(DMPlexCreateFromFile(comm, filename, "ex11_plex", PETSC_TRUE, &dm)); 197 PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE)); 198 PetscCall(DMCreateLabel(dm, name)); 199 PetscCall(DMGetLabel(dm, name, &label)); 200 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 201 for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelSetValue(label, c, c)); 202 PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD)); 203 PetscCall(DMPlexGetPartitioner(dm, &part)); 204 PetscCall(PetscPartitionerSetFromOptions(part)); 205 PetscCall(DMPlexDistribute(dm, overlap, NULL, &dmDist)); 206 if (dmDist) { 207 PetscCall(DMDestroy(&dm)); 208 dm = dmDist; 209 } 210 PetscCall(PetscObjectSetName((PetscObject)dm, "Mesh")); 211 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 212 PetscCall(DMGetLabel(dm, name, &label)); 213 PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD)); 214 PetscCall(DMDestroy(&dm)); 215 PetscFunctionReturn(PETSC_SUCCESS); 216 } 217 218 static PetscErrorCode TestUniversalLabel(MPI_Comm comm) 219 { 220 DM dm1, dm2; 221 DMLabel bd1, bd2, ulabel; 222 DMUniversalLabel universal; 223 PetscInt pStart, pEnd, p; 224 PetscBool run = PETSC_FALSE, notFile; 225 226 PetscFunctionBeginUser; 227 PetscCall(PetscOptionsGetBool(NULL, NULL, "-universal", &run, NULL)); 228 if (!run) PetscFunctionReturn(PETSC_SUCCESS); 229 230 char filename[PETSC_MAX_PATH_LEN]; 231 PetscBool flg; 232 233 PetscCall(PetscOptionsGetString(NULL, NULL, "-filename", filename, sizeof(filename), &flg)); 234 if (flg) { 235 PetscCall(DMPlexCreateFromFile(comm, filename, "ex11_plex", PETSC_TRUE, &dm1)); 236 } else { 237 PetscCall(DMCreate(comm, &dm1)); 238 PetscCall(DMSetType(dm1, DMPLEX)); 239 PetscCall(DMSetFromOptions(dm1)); 240 } 241 PetscCall(DMHasLabel(dm1, "marker", ¬File)); 242 if (notFile) { 243 PetscCall(DMCreateLabel(dm1, "Boundary Faces")); 244 PetscCall(DMGetLabel(dm1, "Boundary Faces", &bd1)); 245 PetscCall(DMPlexMarkBoundaryFaces(dm1, 13, bd1)); 246 PetscCall(DMCreateLabel(dm1, "Boundary")); 247 PetscCall(DMGetLabel(dm1, "Boundary", &bd2)); 248 PetscCall(DMPlexMarkBoundaryFaces(dm1, 121, bd2)); 249 PetscCall(DMPlexLabelComplete(dm1, bd2)); 250 } 251 PetscCall(PetscObjectSetName((PetscObject)dm1, "First Mesh")); 252 PetscCall(DMViewFromOptions(dm1, NULL, "-dm_view")); 253 254 PetscCall(DMUniversalLabelCreate(dm1, &universal)); 255 PetscCall(DMUniversalLabelGetLabel(universal, &ulabel)); 256 PetscCall(PetscObjectViewFromOptions((PetscObject)ulabel, NULL, "-universal_view")); 257 258 if (!notFile) { 259 PetscInt Nl, l; 260 261 PetscCall(DMClone(dm1, &dm2)); 262 PetscCall(DMGetNumLabels(dm2, &Nl)); 263 for (l = Nl - 1; l >= 0; --l) { 264 PetscBool isdepth, iscelltype; 265 const char *name; 266 267 PetscCall(DMGetLabelName(dm2, l, &name)); 268 PetscCall(PetscStrncmp(name, "depth", 6, &isdepth)); 269 PetscCall(PetscStrncmp(name, "celltype", 9, &iscelltype)); 270 if (!isdepth && !iscelltype) PetscCall(DMRemoveLabel(dm2, name, NULL)); 271 } 272 } else { 273 PetscCall(DMCreate(comm, &dm2)); 274 PetscCall(DMSetType(dm2, DMPLEX)); 275 PetscCall(DMSetFromOptions(dm2)); 276 } 277 PetscCall(PetscObjectSetName((PetscObject)dm2, "Second Mesh")); 278 PetscCall(DMUniversalLabelCreateLabels(universal, PETSC_TRUE, dm2)); 279 PetscCall(DMPlexGetChart(dm2, &pStart, &pEnd)); 280 for (p = pStart; p < pEnd; ++p) { 281 PetscInt val; 282 283 PetscCall(DMLabelGetValue(ulabel, p, &val)); 284 if (val < 0) continue; 285 PetscCall(DMUniversalLabelSetLabelValue(universal, dm2, PETSC_TRUE, p, val)); 286 } 287 PetscCall(DMViewFromOptions(dm2, NULL, "-dm_view")); 288 289 PetscCall(DMUniversalLabelDestroy(&universal)); 290 PetscCall(DMDestroy(&dm1)); 291 PetscCall(DMDestroy(&dm2)); 292 PetscFunctionReturn(PETSC_SUCCESS); 293 } 294 295 int main(int argc, char **argv) 296 { 297 PetscFunctionBeginUser; 298 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 299 /*PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));*/ 300 PetscCall(TestInsertion()); 301 PetscCall(TestEmptyStrata(PETSC_COMM_WORLD)); 302 PetscCall(TestDistribution(PETSC_COMM_WORLD)); 303 PetscCall(TestUniversalLabel(PETSC_COMM_WORLD)); 304 PetscCall(PetscFinalize()); 305 return 0; 306 } 307 308 /*TEST 309 310 test: 311 suffix: 0 312 requires: triangle 313 test: 314 suffix: 1 315 requires: triangle 316 nsize: 2 317 args: -petscpartitioner_type simple 318 319 testset: 320 suffix: gmsh 321 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square.msh -petscpartitioner_type simple 322 test: 323 suffix: 1 324 nsize: 1 325 test: 326 suffix: 2 327 nsize: 2 328 329 testset: 330 suffix: exodusii 331 requires: exodusii 332 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/2Dgrd.exo -petscpartitioner_type simple 333 test: 334 suffix: 1 335 nsize: 1 336 test: 337 suffix: 2 338 nsize: 2 339 340 test: 341 suffix: univ 342 requires: triangle 343 args: -universal -dm_view -universal_view 344 345 test: 346 # Note that the labels differ because we have multiply-marked some points during EGADS creation 347 suffix: univ_egads_sphere 348 requires: egads datafilespath 349 args: -universal -dm_plex_filename ${DATAFILESPATH}/meshes/cad/sphere_example.egadslite -dm_view -universal_view 350 351 test: 352 # Note that the labels differ because we have multiply-marked some points during EGADS creation 353 suffix: univ_egads_ball 354 requires: egads ctetgen datafilespath 355 args: -universal -dm_plex_boundary_filename ${DATAFILESPATH}/meshes/cad/sphere_example.egadslite -dm_view -universal_view 356 357 TEST*/ 358