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) { 232 cval[offUA+i] += xyz[j*sdim+i]; 233 } 234 cval[offUA+i] = cval[offUA+i] * sdim / clSize; 235 cval[offUA+sdim] += PetscSqr(cval[offUA+i]); 236 } 237 PetscCall(DMPlexVecRestoreClosure(seqdmUA, coordSection, coord, p, &clSize, &xyz)); 238 } 239 } 240 PetscCall(VecRestoreArray(UALoc, &cval)); 241 PetscCall(DMLocalToGlobalBegin(seqdmUA, UALoc, INSERT_VALUES, seqUA)); 242 PetscCall(DMLocalToGlobalEnd(seqdmUA, UALoc, INSERT_VALUES, seqUA)); 243 PetscCall(DMRestoreLocalVector(seqdmUA, &UALoc)); 244 245 /* Update X */ 246 PetscCall(VecISCopy(seqX, seqisUA, SCATTER_FORWARD, seqUA)); 247 PetscCall(VecViewFromOptions(seqUA, NULL, "-ua_vec_view")); 248 249 /* Restrict to U and Alpha */ 250 PetscCall(VecISCopy(seqX, seqisU, SCATTER_REVERSE, seqU)); 251 PetscCall(VecISCopy(seqX, seqisA, SCATTER_REVERSE, seqA)); 252 253 /* restrict to UA2 */ 254 PetscCall(VecISCopy(seqX, seqisUA, SCATTER_REVERSE, seqUA2)); 255 PetscCall(VecViewFromOptions(seqUA2, NULL, "-ua2_vec_view")); 256 } 257 258 { 259 PetscInt ovlp = 0; 260 PetscPartitioner part; 261 262 PetscCall(DMSetUseNatural(dm,PETSC_TRUE)); 263 PetscCall(DMPlexGetPartitioner(dm,&part)); 264 PetscCall(PetscPartitionerSetFromOptions(part)); 265 PetscCall(DMPlexDistribute(dm,ovlp,&migrationSF,&pdm)); 266 if (!pdm) pdm = dm; 267 if (migrationSF) { 268 PetscCall(DMPlexSetMigrationSF(pdm, migrationSF)); 269 PetscCall(PetscSFDestroy(&migrationSF)); 270 } 271 PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view")); 272 } 273 274 /* Get DM and IS for each field of dm */ 275 PetscCall(DMCreateSubDM(pdm, 1, &fieldU, &isU, &dmU)); 276 PetscCall(DMCreateSubDM(pdm, 1, &fieldA, &isA, &dmA)); 277 PetscCall(DMCreateSubDM(pdm, 1, &fieldS, &isS, &dmS)); 278 PetscCall(DMCreateSubDM(pdm, 2, fieldUA, &isUA, &dmUA)); 279 280 PetscCall(PetscMalloc1(2,&dmList)); 281 dmList[0] = dmU; 282 dmList[1] = dmA; 283 284 PetscCall(DMCreateSuperDM(dmList,2,NULL,&dmUA2)); 285 PetscCall(PetscFree(dmList)); 286 287 PetscCall(DMGetGlobalVector(pdm, &X)); 288 PetscCall(DMGetGlobalVector(dmU, &U)); 289 PetscCall(DMGetGlobalVector(dmA, &A)); 290 PetscCall(DMGetGlobalVector(dmS, &S)); 291 PetscCall(DMGetGlobalVector(dmUA, &UA)); 292 PetscCall(DMGetGlobalVector(dmUA2, &UA2)); 293 294 PetscCall(PetscObjectSetName((PetscObject) X, "X")); 295 PetscCall(PetscObjectSetName((PetscObject) U, "U")); 296 PetscCall(PetscObjectSetName((PetscObject) A, "Alpha")); 297 PetscCall(PetscObjectSetName((PetscObject) S, "Sigma")); 298 PetscCall(PetscObjectSetName((PetscObject) UA, "UAlpha")); 299 PetscCall(PetscObjectSetName((PetscObject) UA2, "UAlpha2")); 300 PetscCall(VecSet(X, -111.)); 301 302 /* 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 */ 303 { 304 PetscSection sectionUA; 305 Vec UALoc; 306 PetscSection coordSection; 307 Vec coord; 308 PetscScalar *cval, *xyz; 309 PetscInt clSize, i, j; 310 311 PetscCall(DMGetLocalSection(dmUA, §ionUA)); 312 PetscCall(DMGetLocalVector(dmUA, &UALoc)); 313 PetscCall(VecGetArray(UALoc, &cval)); 314 PetscCall(DMGetCoordinateSection(dmUA, &coordSection)); 315 PetscCall(DMGetCoordinatesLocal(dmUA, &coord)); 316 PetscCall(DMPlexGetChart(dmUA, &pStart, &pEnd)); 317 318 for (p = pStart; p < pEnd; ++p) { 319 PetscInt dofUA, offUA; 320 321 PetscCall(PetscSectionGetDof(sectionUA, p, &dofUA)); 322 if (dofUA > 0) { 323 xyz=NULL; 324 PetscCall(PetscSectionGetOffset(sectionUA, p, &offUA)); 325 PetscCall(DMPlexVecGetClosure(dmUA, coordSection, coord, p, &clSize, &xyz)); 326 cval[offUA+sdim] = 0.; 327 for (i = 0; i < sdim; ++i) { 328 cval[offUA+i] = 0; 329 for (j = 0; j < clSize/sdim; ++j) { 330 cval[offUA+i] += xyz[j*sdim+i]; 331 } 332 cval[offUA+i] = cval[offUA+i] * sdim / clSize; 333 cval[offUA+sdim] += PetscSqr(cval[offUA+i]); 334 } 335 PetscCall(DMPlexVecRestoreClosure(dmUA, coordSection, coord, p, &clSize, &xyz)); 336 } 337 } 338 PetscCall(VecRestoreArray(UALoc, &cval)); 339 PetscCall(DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA)); 340 PetscCall(DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA)); 341 PetscCall(DMRestoreLocalVector(dmUA, &UALoc)); 342 343 /* Update X */ 344 PetscCall(VecISCopy(X, isUA, SCATTER_FORWARD, UA)); 345 PetscCall(VecViewFromOptions(UA, NULL, "-ua_vec_view")); 346 347 /* Restrict to U and Alpha */ 348 PetscCall(VecISCopy(X, isU, SCATTER_REVERSE, U)); 349 PetscCall(VecISCopy(X, isA, SCATTER_REVERSE, A)); 350 351 /* restrict to UA2 */ 352 PetscCall(VecISCopy(X, isUA, SCATTER_REVERSE, UA2)); 353 PetscCall(VecViewFromOptions(UA2, NULL, "-ua2_vec_view")); 354 } 355 356 /* Verification */ 357 358 Vec Xnat, Unat, Anat, UAnat, Snat, UA2nat; 359 PetscReal norm; 360 PetscCall(DMGetGlobalVector(dm, &Xnat)); 361 PetscCall(DMGetGlobalVector(seqdmU, &Unat)); 362 PetscCall(DMGetGlobalVector(seqdmA, &Anat)); 363 PetscCall(DMGetGlobalVector(seqdmUA, &UAnat)); 364 PetscCall(DMGetGlobalVector(seqdmS, &Snat)); 365 PetscCall(DMGetGlobalVector(seqdmUA2, &UA2nat)); 366 367 PetscCall(DMPlexGlobalToNaturalBegin(pdm, X, Xnat)); 368 PetscCall(DMPlexGlobalToNaturalEnd(pdm, X, Xnat)); 369 PetscCall(VecAXPY(Xnat, -1.0, seqX)); 370 PetscCall(VecNorm(Xnat, NORM_INFINITY, &norm)); 371 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "X ||Vin - Vout|| = %g", (double) norm); 372 373 PetscCall(DMPlexGlobalToNaturalBegin(dmU, U, Unat)); 374 PetscCall(DMPlexGlobalToNaturalEnd(dmU, U, Unat)); 375 PetscCall(VecAXPY(Unat, -1.0, seqU)); 376 PetscCall(VecNorm(Unat, NORM_INFINITY, &norm)); 377 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "U ||Vin - Vout|| = %g", (double) norm); 378 379 PetscCall(DMPlexGlobalToNaturalBegin(dmA, A, Anat)); 380 PetscCall(DMPlexGlobalToNaturalEnd(dmA, A, Anat)); 381 PetscCall(VecAXPY(Anat, -1.0, seqA)); 382 PetscCall(VecNorm(Anat, NORM_INFINITY, &norm)); 383 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "A ||Vin - Vout|| = %g", (double) norm); 384 385 PetscCall(DMPlexGlobalToNaturalBegin(dmUA, UA, UAnat)); 386 PetscCall(DMPlexGlobalToNaturalEnd(dmUA, UA, UAnat)); 387 PetscCall(VecAXPY(UAnat, -1.0, seqUA)); 388 PetscCall(VecNorm(UAnat, NORM_INFINITY, &norm)); 389 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "UA ||Vin - Vout|| = %g", (double) norm); 390 391 PetscCall(DMPlexGlobalToNaturalBegin(dmS, S, Snat)); 392 PetscCall(DMPlexGlobalToNaturalEnd(dmS, S, Snat)); 393 PetscCall(VecAXPY(Snat, -1.0, seqS)); 394 PetscCall(VecNorm(Snat, NORM_INFINITY, &norm)); 395 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "S ||Vin - Vout|| = %g", (double) norm); 396 397 PetscCall(DMPlexGlobalToNaturalBegin(dmUA2, UA2, UA2nat)); 398 PetscCall(DMPlexGlobalToNaturalEnd(dmUA2, UA2, UA2nat)); 399 PetscCall(VecAXPY(UA2nat, -1.0, seqUA2)); 400 PetscCall(VecNorm(UA2nat, NORM_INFINITY, &norm)); 401 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "UAlpha2 ||Vin - Vout|| = %g", (double) norm); 402 403 PetscCall(DMRestoreGlobalVector(dm, &Xnat)); 404 PetscCall(DMRestoreGlobalVector(seqdmU, &Unat)); 405 PetscCall(DMRestoreGlobalVector(seqdmA, &Anat)); 406 PetscCall(DMRestoreGlobalVector(seqdmUA, &UAnat)); 407 PetscCall(DMRestoreGlobalVector(seqdmS, &Snat)); 408 PetscCall(DMRestoreGlobalVector(seqdmUA2, &UA2nat)); 409 410 PetscCall(DMRestoreGlobalVector(seqdmUA2, &seqUA2)); 411 PetscCall(DMRestoreGlobalVector(seqdmUA, &seqUA)); 412 PetscCall(DMRestoreGlobalVector(seqdmS, &seqS)); 413 PetscCall(DMRestoreGlobalVector(seqdmA, &seqA)); 414 PetscCall(DMRestoreGlobalVector(seqdmU, &seqU)); 415 PetscCall(DMRestoreGlobalVector(dm, &seqX)); 416 PetscCall(DMDestroy(&seqdmU));PetscCall(ISDestroy(&seqisU)); 417 PetscCall(DMDestroy(&seqdmA));PetscCall(ISDestroy(&seqisA)); 418 PetscCall(DMDestroy(&seqdmS));PetscCall(ISDestroy(&seqisS)); 419 PetscCall(DMDestroy(&seqdmUA));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)); PetscCall(ISDestroy(&isU)); 429 PetscCall(DMDestroy(&dmA)); PetscCall(ISDestroy(&isA)); 430 PetscCall(DMDestroy(&dmS)); PetscCall(ISDestroy(&isS)); 431 PetscCall(DMDestroy(&dmUA));PetscCall(ISDestroy(&isUA)); 432 PetscCall(DMDestroy(&dmUA2)); 433 PetscCall(DMDestroy(&pdm)); 434 if (size > 1) PetscCall(DMDestroy(&dm)); 435 PetscCall(PetscFree2(pStartDepth, pEndDepth)); 436 PetscCall(PetscFinalize()); 437 return 0; 438 } 439 440 /*TEST 441 442 build: 443 requires: !complex parmetis exodusii pnetcdf 444 # 2D seq 445 test: 446 suffix: 0 447 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -dm_view -petscpartitioner_type parmetis 448 test: 449 suffix: 1 450 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -dm_view -petscpartitioner_type parmetis 451 test: 452 suffix: 2 453 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -dm_view -petscpartitioner_type parmetis 454 455 # 2D par 456 test: 457 suffix: 3 458 nsize: 2 459 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -dm_view -petscpartitioner_type parmetis 460 test: 461 suffix: 4 462 nsize: 2 463 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -dm_view -petscpartitioner_type parmetis 464 test: 465 suffix: 5 466 nsize: 2 467 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -dm_view -petscpartitioner_type parmetis 468 469 #3d seq 470 test: 471 suffix: 6 472 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -dm_view -petscpartitioner_type parmetis 473 test: 474 suffix: 7 475 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -dm_view -petscpartitioner_type parmetis 476 477 #3d par 478 test: 479 suffix: 8 480 nsize: 2 481 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -dm_view -petscpartitioner_type parmetis 482 test: 483 suffix: 9 484 nsize: 2 485 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -dm_view -petscpartitioner_type parmetis 486 487 TEST*/ 488