1 static char help[] = "Test FEM layout and GlobalToNaturalSF\n\n"; 2 3 /* 4 In order to see the vectors which are being tested, use 5 6 -ua_vec_view -s_vec_view 7 */ 8 9 #include <petsc.h> 10 #include <exodusII.h> 11 12 #include <petsc/private/dmpleximpl.h> 13 14 int main(int argc, char **argv) { 15 DM dm, pdm, dmU, dmA, dmS, dmUA, dmUA2, *dmList; 16 DM seqdmU, seqdmA, seqdmS, seqdmUA, seqdmUA2, *seqdmList; 17 Vec X, U, A, S, UA, UA2; 18 Vec seqX, seqU, seqA, seqS, seqUA, seqUA2; 19 IS isU, isA, isS, isUA; 20 IS seqisU, seqisA, seqisS, seqisUA; 21 PetscSection section; 22 const PetscInt fieldU = 0; 23 const PetscInt fieldA = 2; 24 const PetscInt fieldS = 1; 25 const PetscInt fieldUA[2] = {0, 2}; 26 char ifilename[PETSC_MAX_PATH_LEN]; 27 IS csIS; 28 const PetscInt *csID; 29 PetscInt *pStartDepth, *pEndDepth; 30 PetscInt order = 1; 31 PetscInt sdim, d, pStart, pEnd, p, numCS, set; 32 PetscMPIInt rank, size; 33 PetscSF migrationSF; 34 35 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 36 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 37 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 38 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FEM Layout Options", "ex64"); 39 PetscCall(PetscOptionsString("-i", "Filename to read", "ex64", ifilename, ifilename, sizeof(ifilename), NULL)); 40 PetscOptionsEnd(); 41 42 /* Read the mesh from a file in any supported format */ 43 PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename, "ex64_plex", PETSC_TRUE, &dm)); 44 PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE)); 45 PetscCall(DMSetFromOptions(dm)); 46 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 47 PetscCall(DMGetDimension(dm, &sdim)); 48 49 /* Create the main section containing all fields */ 50 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion)); 51 PetscCall(PetscSectionSetNumFields(section, 3)); 52 PetscCall(PetscSectionSetFieldName(section, fieldU, "U")); 53 PetscCall(PetscSectionSetFieldName(section, fieldA, "Alpha")); 54 PetscCall(PetscSectionSetFieldName(section, fieldS, "Sigma")); 55 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 56 PetscCall(PetscSectionSetChart(section, pStart, pEnd)); 57 PetscCall(PetscMalloc2(sdim + 1, &pStartDepth, sdim + 1, &pEndDepth)); 58 for (d = 0; d <= sdim; ++d) PetscCall(DMPlexGetDepthStratum(dm, d, &pStartDepth[d], &pEndDepth[d])); 59 /* Vector field U, Scalar field Alpha, Tensor field Sigma */ 60 PetscCall(PetscSectionSetFieldComponents(section, fieldU, sdim)); 61 PetscCall(PetscSectionSetFieldComponents(section, fieldA, 1)); 62 PetscCall(PetscSectionSetFieldComponents(section, fieldS, sdim * (sdim + 1) / 2)); 63 64 /* Going through cell sets then cells, and setting up storage for the sections */ 65 PetscCall(DMGetLabelSize(dm, "Cell Sets", &numCS)); 66 PetscCall(DMGetLabelIdIS(dm, "Cell Sets", &csIS)); 67 if (csIS) PetscCall(ISGetIndices(csIS, &csID)); 68 for (set = 0; set < numCS; set++) { 69 IS cellIS; 70 const PetscInt *cellID; 71 PetscInt numCells, cell, closureSize, *closureA = NULL; 72 73 PetscCall(DMGetStratumSize(dm, "Cell Sets", csID[set], &numCells)); 74 PetscCall(DMGetStratumIS(dm, "Cell Sets", csID[set], &cellIS)); 75 if (numCells > 0) { 76 /* dof layout ordered by increasing height in the DAG: cell, face, edge, vertex */ 77 PetscInt dofUP1Tri[] = {2, 0, 0}; 78 PetscInt dofAP1Tri[] = {1, 0, 0}; 79 PetscInt dofUP2Tri[] = {2, 2, 0}; 80 PetscInt dofAP2Tri[] = {1, 1, 0}; 81 PetscInt dofUP1Quad[] = {2, 0, 0}; 82 PetscInt dofAP1Quad[] = {1, 0, 0}; 83 PetscInt dofUP2Quad[] = {2, 2, 2}; 84 PetscInt dofAP2Quad[] = {1, 1, 1}; 85 PetscInt dofS2D[] = {0, 0, 3}; 86 PetscInt dofUP1Tet[] = {3, 0, 0, 0}; 87 PetscInt dofAP1Tet[] = {1, 0, 0, 0}; 88 PetscInt dofUP2Tet[] = {3, 3, 0, 0}; 89 PetscInt dofAP2Tet[] = {1, 1, 0, 0}; 90 PetscInt dofUP1Hex[] = {3, 0, 0, 0}; 91 PetscInt dofAP1Hex[] = {1, 0, 0, 0}; 92 PetscInt dofUP2Hex[] = {3, 3, 3, 3}; 93 PetscInt dofAP2Hex[] = {1, 1, 1, 1}; 94 PetscInt dofS3D[] = {0, 0, 0, 6}; 95 PetscInt *dofU, *dofA, *dofS; 96 97 switch (sdim) { 98 case 2: dofS = dofS2D; break; 99 case 3: dofS = dofS3D; break; 100 default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No layout for dimension %" PetscInt_FMT, sdim); 101 } 102 103 /* Identify cell type based on closure size only. This works for Tri/Tet/Quad/Hex meshes 104 It will not be enough to identify more exotic elements like pyramid or prisms... */ 105 PetscCall(ISGetIndices(cellIS, &cellID)); 106 PetscCall(DMPlexGetTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA)); 107 switch (closureSize) { 108 case 7: /* Tri */ 109 if (order == 1) { 110 dofU = dofUP1Tri; 111 dofA = dofAP1Tri; 112 } else { 113 dofU = dofUP2Tri; 114 dofA = dofAP2Tri; 115 } 116 break; 117 case 9: /* Quad */ 118 if (order == 1) { 119 dofU = dofUP1Quad; 120 dofA = dofAP1Quad; 121 } else { 122 dofU = dofUP2Quad; 123 dofA = dofAP2Quad; 124 } 125 break; 126 case 15: /* Tet */ 127 if (order == 1) { 128 dofU = dofUP1Tet; 129 dofA = dofAP1Tet; 130 } else { 131 dofU = dofUP2Tet; 132 dofA = dofAP2Tet; 133 } 134 break; 135 case 27: /* Hex */ 136 if (order == 1) { 137 dofU = dofUP1Hex; 138 dofA = dofAP1Hex; 139 } else { 140 dofU = dofUP2Hex; 141 dofA = dofAP2Hex; 142 } 143 break; 144 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Unknown element with closure size %" PetscInt_FMT, closureSize); 145 } 146 PetscCall(DMPlexRestoreTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA)); 147 148 for (cell = 0; cell < numCells; cell++) { 149 PetscInt *closure = NULL; 150 151 PetscCall(DMPlexGetTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure)); 152 for (p = 0; p < closureSize; ++p) { 153 /* Find depth of p */ 154 for (d = 0; d <= sdim; ++d) { 155 if ((closure[2 * p] >= pStartDepth[d]) && (closure[2 * p] < pEndDepth[d])) { 156 PetscCall(PetscSectionSetDof(section, closure[2 * p], dofU[d] + dofA[d] + dofS[d])); 157 PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldU, dofU[d])); 158 PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldA, dofA[d])); 159 PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldS, dofS[d])); 160 } 161 } 162 } 163 PetscCall(DMPlexRestoreTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure)); 164 } 165 PetscCall(ISRestoreIndices(cellIS, &cellID)); 166 PetscCall(ISDestroy(&cellIS)); 167 } 168 } 169 if (csIS) PetscCall(ISRestoreIndices(csIS, &csID)); 170 PetscCall(ISDestroy(&csIS)); 171 PetscCall(PetscSectionSetUp(section)); 172 PetscCall(DMSetLocalSection(dm, section)); 173 PetscCall(PetscObjectViewFromOptions((PetscObject)section, NULL, "-dm_section_view")); 174 PetscCall(PetscSectionDestroy(§ion)); 175 176 /* Get DM and IS for each field of dm */ 177 PetscCall(DMCreateSubDM(dm, 1, &fieldU, &seqisU, &seqdmU)); 178 PetscCall(DMCreateSubDM(dm, 1, &fieldA, &seqisA, &seqdmA)); 179 PetscCall(DMCreateSubDM(dm, 1, &fieldS, &seqisS, &seqdmS)); 180 PetscCall(DMCreateSubDM(dm, 2, fieldUA, &seqisUA, &seqdmUA)); 181 182 PetscCall(PetscMalloc1(2, &seqdmList)); 183 seqdmList[0] = seqdmU; 184 seqdmList[1] = seqdmA; 185 186 PetscCall(DMCreateSuperDM(seqdmList, 2, NULL, &seqdmUA2)); 187 PetscCall(PetscFree(seqdmList)); 188 189 PetscCall(DMGetGlobalVector(dm, &seqX)); 190 PetscCall(DMGetGlobalVector(seqdmU, &seqU)); 191 PetscCall(DMGetGlobalVector(seqdmA, &seqA)); 192 PetscCall(DMGetGlobalVector(seqdmS, &seqS)); 193 PetscCall(DMGetGlobalVector(seqdmUA, &seqUA)); 194 PetscCall(DMGetGlobalVector(seqdmUA2, &seqUA2)); 195 196 PetscCall(PetscObjectSetName((PetscObject)seqX, "seqX")); 197 PetscCall(PetscObjectSetName((PetscObject)seqU, "seqU")); 198 PetscCall(PetscObjectSetName((PetscObject)seqA, "seqAlpha")); 199 PetscCall(PetscObjectSetName((PetscObject)seqS, "seqSigma")); 200 PetscCall(PetscObjectSetName((PetscObject)seqUA, "seqUAlpha")); 201 PetscCall(PetscObjectSetName((PetscObject)seqUA2, "seqUAlpha2")); 202 PetscCall(VecSet(seqX, -111.)); 203 204 /* Setting u to [x,y,z] and alpha to x^2+y^2+z^2 by writing in UAlpha then restricting to U and Alpha */ 205 { 206 PetscSection sectionUA; 207 Vec UALoc; 208 PetscSection coordSection; 209 Vec coord; 210 PetscScalar *cval, *xyz; 211 PetscInt clSize, i, j; 212 213 PetscCall(DMGetLocalSection(seqdmUA, §ionUA)); 214 PetscCall(DMGetLocalVector(seqdmUA, &UALoc)); 215 PetscCall(VecGetArray(UALoc, &cval)); 216 PetscCall(DMGetCoordinateSection(seqdmUA, &coordSection)); 217 PetscCall(DMGetCoordinatesLocal(seqdmUA, &coord)); 218 PetscCall(DMPlexGetChart(seqdmUA, &pStart, &pEnd)); 219 220 for (p = pStart; p < pEnd; ++p) { 221 PetscInt dofUA, offUA; 222 223 PetscCall(PetscSectionGetDof(sectionUA, p, &dofUA)); 224 if (dofUA > 0) { 225 xyz = NULL; 226 PetscCall(PetscSectionGetOffset(sectionUA, p, &offUA)); 227 PetscCall(DMPlexVecGetClosure(seqdmUA, coordSection, coord, p, &clSize, &xyz)); 228 cval[offUA + sdim] = 0.; 229 for (i = 0; i < sdim; ++i) { 230 cval[offUA + i] = 0; 231 for (j = 0; j < clSize / sdim; ++j) cval[offUA + i] += xyz[j * sdim + i]; 232 cval[offUA + i] = cval[offUA + i] * sdim / clSize; 233 cval[offUA + sdim] += PetscSqr(cval[offUA + i]); 234 } 235 PetscCall(DMPlexVecRestoreClosure(seqdmUA, coordSection, coord, p, &clSize, &xyz)); 236 } 237 } 238 PetscCall(VecRestoreArray(UALoc, &cval)); 239 PetscCall(DMLocalToGlobalBegin(seqdmUA, UALoc, INSERT_VALUES, seqUA)); 240 PetscCall(DMLocalToGlobalEnd(seqdmUA, UALoc, INSERT_VALUES, seqUA)); 241 PetscCall(DMRestoreLocalVector(seqdmUA, &UALoc)); 242 243 /* Update X */ 244 PetscCall(VecISCopy(seqX, seqisUA, SCATTER_FORWARD, seqUA)); 245 PetscCall(VecViewFromOptions(seqUA, NULL, "-ua_vec_view")); 246 247 /* Restrict to U and Alpha */ 248 PetscCall(VecISCopy(seqX, seqisU, SCATTER_REVERSE, seqU)); 249 PetscCall(VecISCopy(seqX, seqisA, SCATTER_REVERSE, seqA)); 250 251 /* restrict to UA2 */ 252 PetscCall(VecISCopy(seqX, seqisUA, SCATTER_REVERSE, seqUA2)); 253 PetscCall(VecViewFromOptions(seqUA2, NULL, "-ua2_vec_view")); 254 } 255 256 { 257 PetscInt ovlp = 0; 258 PetscPartitioner part; 259 260 PetscCall(DMSetUseNatural(dm, PETSC_TRUE)); 261 PetscCall(DMPlexGetPartitioner(dm, &part)); 262 PetscCall(PetscPartitionerSetFromOptions(part)); 263 PetscCall(DMPlexDistribute(dm, ovlp, &migrationSF, &pdm)); 264 if (!pdm) pdm = dm; 265 if (migrationSF) { 266 PetscCall(DMPlexSetMigrationSF(pdm, migrationSF)); 267 PetscCall(PetscSFDestroy(&migrationSF)); 268 } 269 PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view")); 270 } 271 272 /* Get DM and IS for each field of dm */ 273 PetscCall(DMCreateSubDM(pdm, 1, &fieldU, &isU, &dmU)); 274 PetscCall(DMCreateSubDM(pdm, 1, &fieldA, &isA, &dmA)); 275 PetscCall(DMCreateSubDM(pdm, 1, &fieldS, &isS, &dmS)); 276 PetscCall(DMCreateSubDM(pdm, 2, fieldUA, &isUA, &dmUA)); 277 278 PetscCall(PetscMalloc1(2, &dmList)); 279 dmList[0] = dmU; 280 dmList[1] = dmA; 281 282 PetscCall(DMCreateSuperDM(dmList, 2, NULL, &dmUA2)); 283 PetscCall(PetscFree(dmList)); 284 285 PetscCall(DMGetGlobalVector(pdm, &X)); 286 PetscCall(DMGetGlobalVector(dmU, &U)); 287 PetscCall(DMGetGlobalVector(dmA, &A)); 288 PetscCall(DMGetGlobalVector(dmS, &S)); 289 PetscCall(DMGetGlobalVector(dmUA, &UA)); 290 PetscCall(DMGetGlobalVector(dmUA2, &UA2)); 291 292 PetscCall(PetscObjectSetName((PetscObject)X, "X")); 293 PetscCall(PetscObjectSetName((PetscObject)U, "U")); 294 PetscCall(PetscObjectSetName((PetscObject)A, "Alpha")); 295 PetscCall(PetscObjectSetName((PetscObject)S, "Sigma")); 296 PetscCall(PetscObjectSetName((PetscObject)UA, "UAlpha")); 297 PetscCall(PetscObjectSetName((PetscObject)UA2, "UAlpha2")); 298 PetscCall(VecSet(X, -111.)); 299 300 /* Setting u to [x,y,z] and alpha to x^2+y^2+z^2 by writing in UAlpha then restricting to U and Alpha */ 301 { 302 PetscSection sectionUA; 303 Vec UALoc; 304 PetscSection coordSection; 305 Vec coord; 306 PetscScalar *cval, *xyz; 307 PetscInt clSize, i, j; 308 309 PetscCall(DMGetLocalSection(dmUA, §ionUA)); 310 PetscCall(DMGetLocalVector(dmUA, &UALoc)); 311 PetscCall(VecGetArray(UALoc, &cval)); 312 PetscCall(DMGetCoordinateSection(dmUA, &coordSection)); 313 PetscCall(DMGetCoordinatesLocal(dmUA, &coord)); 314 PetscCall(DMPlexGetChart(dmUA, &pStart, &pEnd)); 315 316 for (p = pStart; p < pEnd; ++p) { 317 PetscInt dofUA, offUA; 318 319 PetscCall(PetscSectionGetDof(sectionUA, p, &dofUA)); 320 if (dofUA > 0) { 321 xyz = NULL; 322 PetscCall(PetscSectionGetOffset(sectionUA, p, &offUA)); 323 PetscCall(DMPlexVecGetClosure(dmUA, coordSection, coord, p, &clSize, &xyz)); 324 cval[offUA + sdim] = 0.; 325 for (i = 0; i < sdim; ++i) { 326 cval[offUA + i] = 0; 327 for (j = 0; j < clSize / sdim; ++j) cval[offUA + i] += xyz[j * sdim + i]; 328 cval[offUA + i] = cval[offUA + i] * sdim / clSize; 329 cval[offUA + sdim] += PetscSqr(cval[offUA + i]); 330 } 331 PetscCall(DMPlexVecRestoreClosure(dmUA, coordSection, coord, p, &clSize, &xyz)); 332 } 333 } 334 PetscCall(VecRestoreArray(UALoc, &cval)); 335 PetscCall(DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA)); 336 PetscCall(DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA)); 337 PetscCall(DMRestoreLocalVector(dmUA, &UALoc)); 338 339 /* Update X */ 340 PetscCall(VecISCopy(X, isUA, SCATTER_FORWARD, UA)); 341 PetscCall(VecViewFromOptions(UA, NULL, "-ua_vec_view")); 342 343 /* Restrict to U and Alpha */ 344 PetscCall(VecISCopy(X, isU, SCATTER_REVERSE, U)); 345 PetscCall(VecISCopy(X, isA, SCATTER_REVERSE, A)); 346 347 /* restrict to UA2 */ 348 PetscCall(VecISCopy(X, isUA, SCATTER_REVERSE, UA2)); 349 PetscCall(VecViewFromOptions(UA2, NULL, "-ua2_vec_view")); 350 } 351 352 /* Verification */ 353 354 Vec Xnat, Unat, Anat, UAnat, Snat, UA2nat; 355 PetscReal norm; 356 PetscCall(DMGetGlobalVector(dm, &Xnat)); 357 PetscCall(DMGetGlobalVector(seqdmU, &Unat)); 358 PetscCall(DMGetGlobalVector(seqdmA, &Anat)); 359 PetscCall(DMGetGlobalVector(seqdmUA, &UAnat)); 360 PetscCall(DMGetGlobalVector(seqdmS, &Snat)); 361 PetscCall(DMGetGlobalVector(seqdmUA2, &UA2nat)); 362 363 PetscCall(DMPlexGlobalToNaturalBegin(pdm, X, Xnat)); 364 PetscCall(DMPlexGlobalToNaturalEnd(pdm, X, Xnat)); 365 PetscCall(VecAXPY(Xnat, -1.0, seqX)); 366 PetscCall(VecNorm(Xnat, NORM_INFINITY, &norm)); 367 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "X ||Vin - Vout|| = %g", (double)norm); 368 369 PetscCall(DMPlexGlobalToNaturalBegin(dmU, U, Unat)); 370 PetscCall(DMPlexGlobalToNaturalEnd(dmU, U, Unat)); 371 PetscCall(VecAXPY(Unat, -1.0, seqU)); 372 PetscCall(VecNorm(Unat, NORM_INFINITY, &norm)); 373 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "U ||Vin - Vout|| = %g", (double)norm); 374 375 PetscCall(DMPlexGlobalToNaturalBegin(dmA, A, Anat)); 376 PetscCall(DMPlexGlobalToNaturalEnd(dmA, A, Anat)); 377 PetscCall(VecAXPY(Anat, -1.0, seqA)); 378 PetscCall(VecNorm(Anat, NORM_INFINITY, &norm)); 379 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "A ||Vin - Vout|| = %g", (double)norm); 380 381 PetscCall(DMPlexGlobalToNaturalBegin(dmUA, UA, UAnat)); 382 PetscCall(DMPlexGlobalToNaturalEnd(dmUA, UA, UAnat)); 383 PetscCall(VecAXPY(UAnat, -1.0, seqUA)); 384 PetscCall(VecNorm(UAnat, NORM_INFINITY, &norm)); 385 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "UA ||Vin - Vout|| = %g", (double)norm); 386 387 PetscCall(DMPlexGlobalToNaturalBegin(dmS, S, Snat)); 388 PetscCall(DMPlexGlobalToNaturalEnd(dmS, S, Snat)); 389 PetscCall(VecAXPY(Snat, -1.0, seqS)); 390 PetscCall(VecNorm(Snat, NORM_INFINITY, &norm)); 391 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "S ||Vin - Vout|| = %g", (double)norm); 392 393 PetscCall(DMPlexGlobalToNaturalBegin(dmUA2, UA2, UA2nat)); 394 PetscCall(DMPlexGlobalToNaturalEnd(dmUA2, UA2, UA2nat)); 395 PetscCall(VecAXPY(UA2nat, -1.0, seqUA2)); 396 PetscCall(VecNorm(UA2nat, NORM_INFINITY, &norm)); 397 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "UAlpha2 ||Vin - Vout|| = %g", (double)norm); 398 399 PetscCall(DMRestoreGlobalVector(dm, &Xnat)); 400 PetscCall(DMRestoreGlobalVector(seqdmU, &Unat)); 401 PetscCall(DMRestoreGlobalVector(seqdmA, &Anat)); 402 PetscCall(DMRestoreGlobalVector(seqdmUA, &UAnat)); 403 PetscCall(DMRestoreGlobalVector(seqdmS, &Snat)); 404 PetscCall(DMRestoreGlobalVector(seqdmUA2, &UA2nat)); 405 406 PetscCall(DMRestoreGlobalVector(seqdmUA2, &seqUA2)); 407 PetscCall(DMRestoreGlobalVector(seqdmUA, &seqUA)); 408 PetscCall(DMRestoreGlobalVector(seqdmS, &seqS)); 409 PetscCall(DMRestoreGlobalVector(seqdmA, &seqA)); 410 PetscCall(DMRestoreGlobalVector(seqdmU, &seqU)); 411 PetscCall(DMRestoreGlobalVector(dm, &seqX)); 412 PetscCall(DMDestroy(&seqdmU)); 413 PetscCall(ISDestroy(&seqisU)); 414 PetscCall(DMDestroy(&seqdmA)); 415 PetscCall(ISDestroy(&seqisA)); 416 PetscCall(DMDestroy(&seqdmS)); 417 PetscCall(ISDestroy(&seqisS)); 418 PetscCall(DMDestroy(&seqdmUA)); 419 PetscCall(ISDestroy(&seqisUA)); 420 PetscCall(DMDestroy(&seqdmUA2)); 421 422 PetscCall(DMRestoreGlobalVector(dmUA2, &UA2)); 423 PetscCall(DMRestoreGlobalVector(dmUA, &UA)); 424 PetscCall(DMRestoreGlobalVector(dmS, &S)); 425 PetscCall(DMRestoreGlobalVector(dmA, &A)); 426 PetscCall(DMRestoreGlobalVector(dmU, &U)); 427 PetscCall(DMRestoreGlobalVector(pdm, &X)); 428 PetscCall(DMDestroy(&dmU)); 429 PetscCall(ISDestroy(&isU)); 430 PetscCall(DMDestroy(&dmA)); 431 PetscCall(ISDestroy(&isA)); 432 PetscCall(DMDestroy(&dmS)); 433 PetscCall(ISDestroy(&isS)); 434 PetscCall(DMDestroy(&dmUA)); 435 PetscCall(ISDestroy(&isUA)); 436 PetscCall(DMDestroy(&dmUA2)); 437 PetscCall(DMDestroy(&pdm)); 438 if (size > 1) PetscCall(DMDestroy(&dm)); 439 PetscCall(PetscFree2(pStartDepth, pEndDepth)); 440 PetscCall(PetscFinalize()); 441 return 0; 442 } 443 444 /*TEST 445 446 build: 447 requires: !complex parmetis exodusii pnetcdf 448 # 2D seq 449 test: 450 suffix: 0 451 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -dm_view -petscpartitioner_type parmetis 452 test: 453 suffix: 1 454 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -dm_view -petscpartitioner_type parmetis 455 test: 456 suffix: 2 457 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -dm_view -petscpartitioner_type parmetis 458 459 # 2D par 460 test: 461 suffix: 3 462 nsize: 2 463 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -dm_view -petscpartitioner_type parmetis 464 test: 465 suffix: 4 466 nsize: 2 467 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -dm_view -petscpartitioner_type parmetis 468 test: 469 suffix: 5 470 nsize: 2 471 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -dm_view -petscpartitioner_type parmetis 472 473 #3d seq 474 test: 475 suffix: 6 476 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -dm_view -petscpartitioner_type parmetis 477 test: 478 suffix: 7 479 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -dm_view -petscpartitioner_type parmetis 480 481 #3d par 482 test: 483 suffix: 8 484 nsize: 2 485 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -dm_view -petscpartitioner_type parmetis 486 test: 487 suffix: 9 488 nsize: 2 489 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -dm_view -petscpartitioner_type parmetis 490 491 TEST*/ 492