1 static char help[] = "Test FEM layout with DM and ExodusII storage\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 { 16 DM dm, pdm, dmU, dmA, dmS, dmUA, dmUA2, *dmList; 17 Vec X, U, A, S, UA, UA2; 18 IS isU, isA, isS, isUA; 19 PetscSection section; 20 const PetscInt fieldU = 0; 21 const PetscInt fieldA = 2; 22 const PetscInt fieldS = 1; 23 const PetscInt fieldUA[2] = {0, 2}; 24 char ifilename[PETSC_MAX_PATH_LEN], ofilename[PETSC_MAX_PATH_LEN]; 25 int exoid = -1; 26 IS csIS; 27 const PetscInt *csID; 28 PetscInt *pStartDepth, *pEndDepth; 29 PetscInt order = 1; 30 PetscInt sdim, d, pStart, pEnd, p, numCS, set; 31 PetscMPIInt rank, size; 32 PetscViewer viewer; 33 34 PetscFunctionBeginUser; 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", "ex26"); 39 PetscCall(PetscOptionsString("-i", "Filename to read", "ex26", ifilename, ifilename, sizeof(ifilename), NULL)); 40 PetscCall(PetscOptionsString("-o", "Filename to write", "ex26", ofilename, ofilename, sizeof(ofilename), NULL)); 41 PetscCall(PetscOptionsBoundedInt("-order", "FEM polynomial order", "ex26", order, &order, NULL, 1)); 42 PetscOptionsEnd(); 43 PetscCheck((order >= 1) && (order <= 2), PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported polynomial order %" PetscInt_FMT " not in [1, 2]", order); 44 45 /* Read the mesh from a file in any supported format */ 46 PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename, NULL, PETSC_TRUE, &dm)); 47 PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE)); 48 PetscCall(DMSetFromOptions(dm)); 49 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 50 PetscCall(DMGetDimension(dm, &sdim)); 51 52 /* Create the exodus result file */ 53 { 54 PetscInt numstep = 3, step; 55 char *nodalVarName[4]; 56 char *zonalVarName[6]; 57 int *truthtable; 58 PetscInt numNodalVar, numZonalVar, i; 59 60 /* enable exodus debugging information */ 61 ex_opts(EX_VERBOSE | EX_DEBUG); 62 /* Create the exodus file */ 63 PetscCall(PetscViewerExodusIIOpen(PETSC_COMM_WORLD, ofilename, FILE_MODE_WRITE, &viewer)); 64 /* The long way would be */ 65 /* 66 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD,&viewer)); 67 PetscCall(PetscViewerSetType(viewer,PETSCVIEWEREXODUSII)); 68 PetscCall(PetscViewerFileSetMode(viewer,FILE_MODE_APPEND)); 69 PetscCall(PetscViewerFileSetName(viewer,ofilename)); 70 */ 71 /* set the mesh order */ 72 PetscCall(PetscViewerExodusIISetOrder(viewer, order)); 73 PetscCall(PetscViewerView(viewer, PETSC_VIEWER_STDOUT_WORLD)); 74 /* 75 Notice how the exodus file is actually NOT open at this point (exoid is -1) 76 Since we are overwriting the file (mode is FILE_MODE_WRITE), we are going to have to 77 write the geometry (the DM), which can only be done on a brand new file. 78 */ 79 80 /* Save the geometry to the file, erasing all previous content */ 81 PetscCall(DMView(dm, viewer)); 82 PetscCall(PetscViewerView(viewer, PETSC_VIEWER_STDOUT_WORLD)); 83 /* 84 Note how the exodus file is now open 85 */ 86 87 /* "Format" the exodus result file, i.e. allocate space for nodal and zonal variables */ 88 switch (sdim) { 89 case 2: 90 numNodalVar = 3; 91 nodalVarName[0] = (char *)"U_x"; 92 nodalVarName[1] = (char *)"U_y"; 93 nodalVarName[2] = (char *)"Alpha"; 94 numZonalVar = 3; 95 zonalVarName[0] = (char *)"Sigma_11"; 96 zonalVarName[1] = (char *)"Sigma_22"; 97 zonalVarName[2] = (char *)"Sigma_12"; 98 break; 99 case 3: 100 numNodalVar = 4; 101 nodalVarName[0] = (char *)"U_x"; 102 nodalVarName[1] = (char *)"U_y"; 103 nodalVarName[2] = (char *)"U_z"; 104 nodalVarName[3] = (char *)"Alpha"; 105 numZonalVar = 6; 106 zonalVarName[0] = (char *)"Sigma_11"; 107 zonalVarName[1] = (char *)"Sigma_22"; 108 zonalVarName[2] = (char *)"Sigma_33"; 109 zonalVarName[3] = (char *)"Sigma_23"; 110 zonalVarName[4] = (char *)"Sigma_13"; 111 zonalVarName[5] = (char *)"Sigma_12"; 112 break; 113 default: 114 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No layout for dimension %" PetscInt_FMT, sdim); 115 } 116 PetscCall(PetscViewerExodusIIGetId(viewer, &exoid)); 117 PetscCallExternal(ex_put_variable_param, exoid, EX_ELEM_BLOCK, numZonalVar); 118 PetscCallExternal(ex_put_variable_names, exoid, EX_ELEM_BLOCK, numZonalVar, zonalVarName); 119 PetscCallExternal(ex_put_variable_param, exoid, EX_NODAL, numNodalVar); 120 PetscCallExternal(ex_put_variable_names, exoid, EX_NODAL, numNodalVar, nodalVarName); 121 numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK); 122 123 /* 124 An exodusII truth table specifies which fields are saved at which time step 125 It speeds up I/O but reserving space for fieldsin the file ahead of time. 126 */ 127 PetscCall(PetscMalloc1(numZonalVar * numCS, &truthtable)); 128 for (i = 0; i < numZonalVar * numCS; ++i) truthtable[i] = 1; 129 PetscCallExternal(ex_put_truth_table, exoid, EX_ELEM_BLOCK, numCS, numZonalVar, truthtable); 130 PetscCall(PetscFree(truthtable)); 131 132 /* Writing time step information in the file. Note that this is currently broken in the exodus library for netcdf4 (HDF5-based) files */ 133 for (step = 0; step < numstep; ++step) { 134 PetscReal time = step; 135 PetscCallExternal(ex_put_time, exoid, step + 1, &time); 136 } 137 } 138 139 /* Create the main section containing all fields */ 140 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion)); 141 PetscCall(PetscSectionSetNumFields(section, 3)); 142 PetscCall(PetscSectionSetFieldName(section, fieldU, "U")); 143 PetscCall(PetscSectionSetFieldName(section, fieldA, "Alpha")); 144 PetscCall(PetscSectionSetFieldName(section, fieldS, "Sigma")); 145 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 146 PetscCall(PetscSectionSetChart(section, pStart, pEnd)); 147 PetscCall(PetscMalloc2(sdim + 1, &pStartDepth, sdim + 1, &pEndDepth)); 148 for (d = 0; d <= sdim; ++d) PetscCall(DMPlexGetDepthStratum(dm, d, &pStartDepth[d], &pEndDepth[d])); 149 /* Vector field U, Scalar field Alpha, Tensor field Sigma */ 150 PetscCall(PetscSectionSetFieldComponents(section, fieldU, sdim)); 151 PetscCall(PetscSectionSetFieldComponents(section, fieldA, 1)); 152 PetscCall(PetscSectionSetFieldComponents(section, fieldS, sdim * (sdim + 1) / 2)); 153 154 /* Going through cell sets then cells, and setting up storage for the sections */ 155 PetscCall(DMGetLabelSize(dm, "Cell Sets", &numCS)); 156 PetscCall(DMGetLabelIdIS(dm, "Cell Sets", &csIS)); 157 if (csIS) PetscCall(ISGetIndices(csIS, &csID)); 158 for (set = 0; set < numCS; set++) { 159 IS cellIS; 160 const PetscInt *cellID; 161 PetscInt numCells, cell, closureSize, *closureA = NULL; 162 163 PetscCall(DMGetStratumSize(dm, "Cell Sets", csID[set], &numCells)); 164 PetscCall(DMGetStratumIS(dm, "Cell Sets", csID[set], &cellIS)); 165 if (numCells > 0) { 166 /* dof layout ordered by increasing height in the DAG: cell, face, edge, vertex */ 167 PetscInt dofUP1Tri[] = {2, 0, 0}; 168 PetscInt dofAP1Tri[] = {1, 0, 0}; 169 PetscInt dofUP2Tri[] = {2, 2, 0}; 170 PetscInt dofAP2Tri[] = {1, 1, 0}; 171 PetscInt dofUP1Quad[] = {2, 0, 0}; 172 PetscInt dofAP1Quad[] = {1, 0, 0}; 173 PetscInt dofUP2Quad[] = {2, 2, 2}; 174 PetscInt dofAP2Quad[] = {1, 1, 1}; 175 PetscInt dofS2D[] = {0, 0, 3}; 176 PetscInt dofUP1Tet[] = {3, 0, 0, 0}; 177 PetscInt dofAP1Tet[] = {1, 0, 0, 0}; 178 PetscInt dofUP2Tet[] = {3, 3, 0, 0}; 179 PetscInt dofAP2Tet[] = {1, 1, 0, 0}; 180 PetscInt dofUP1Hex[] = {3, 0, 0, 0}; 181 PetscInt dofAP1Hex[] = {1, 0, 0, 0}; 182 PetscInt dofUP2Hex[] = {3, 3, 3, 3}; 183 PetscInt dofAP2Hex[] = {1, 1, 1, 1}; 184 PetscInt dofS3D[] = {0, 0, 0, 6}; 185 PetscInt *dofU, *dofA, *dofS; 186 187 switch (sdim) { 188 case 2: 189 dofS = dofS2D; 190 break; 191 case 3: 192 dofS = dofS3D; 193 break; 194 default: 195 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No layout for dimension %" PetscInt_FMT, sdim); 196 } 197 198 /* Identify cell type based on closure size only. This works for Tri/Tet/Quad/Hex meshes 199 It will not be enough to identify more exotic elements like pyramid or prisms... */ 200 PetscCall(ISGetIndices(cellIS, &cellID)); 201 PetscCall(DMPlexGetTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA)); 202 switch (closureSize) { 203 case 7: /* Tri */ 204 if (order == 1) { 205 dofU = dofUP1Tri; 206 dofA = dofAP1Tri; 207 } else { 208 dofU = dofUP2Tri; 209 dofA = dofAP2Tri; 210 } 211 break; 212 case 9: /* Quad */ 213 if (order == 1) { 214 dofU = dofUP1Quad; 215 dofA = dofAP1Quad; 216 } else { 217 dofU = dofUP2Quad; 218 dofA = dofAP2Quad; 219 } 220 break; 221 case 15: /* Tet */ 222 if (order == 1) { 223 dofU = dofUP1Tet; 224 dofA = dofAP1Tet; 225 } else { 226 dofU = dofUP2Tet; 227 dofA = dofAP2Tet; 228 } 229 break; 230 case 27: /* Hex */ 231 if (order == 1) { 232 dofU = dofUP1Hex; 233 dofA = dofAP1Hex; 234 } else { 235 dofU = dofUP2Hex; 236 dofA = dofAP2Hex; 237 } 238 break; 239 default: 240 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Unknown element with closure size %" PetscInt_FMT, closureSize); 241 } 242 PetscCall(DMPlexRestoreTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA)); 243 244 for (cell = 0; cell < numCells; cell++) { 245 PetscInt *closure = NULL; 246 247 PetscCall(DMPlexGetTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure)); 248 for (p = 0; p < closureSize; ++p) { 249 /* Find depth of p */ 250 for (d = 0; d <= sdim; ++d) { 251 if ((closure[2 * p] >= pStartDepth[d]) && (closure[2 * p] < pEndDepth[d])) { 252 PetscCall(PetscSectionSetDof(section, closure[2 * p], dofU[d] + dofA[d] + dofS[d])); 253 PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldU, dofU[d])); 254 PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldA, dofA[d])); 255 PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldS, dofS[d])); 256 } 257 } 258 } 259 PetscCall(DMPlexRestoreTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure)); 260 } 261 PetscCall(ISRestoreIndices(cellIS, &cellID)); 262 PetscCall(ISDestroy(&cellIS)); 263 } 264 } 265 if (csIS) PetscCall(ISRestoreIndices(csIS, &csID)); 266 PetscCall(ISDestroy(&csIS)); 267 PetscCall(PetscSectionSetUp(section)); 268 PetscCall(DMSetLocalSection(dm, section)); 269 PetscCall(PetscObjectViewFromOptions((PetscObject)section, NULL, "-dm_section_view")); 270 PetscCall(PetscSectionDestroy(§ion)); 271 272 { 273 PetscSF migrationSF; 274 PetscInt ovlp = 0; 275 PetscPartitioner part; 276 277 PetscCall(DMSetUseNatural(dm, PETSC_TRUE)); 278 PetscCall(DMPlexGetPartitioner(dm, &part)); 279 PetscCall(PetscPartitionerSetFromOptions(part)); 280 PetscCall(DMPlexDistribute(dm, ovlp, &migrationSF, &pdm)); 281 if (!pdm) pdm = dm; 282 /* Set the migrationSF is mandatory since useNatural = PETSC_TRUE */ 283 if (migrationSF) { 284 PetscCall(DMPlexSetMigrationSF(pdm, migrationSF)); 285 PetscCall(PetscSFDestroy(&migrationSF)); 286 } 287 PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view")); 288 } 289 290 /* Get DM and IS for each field of dm */ 291 PetscCall(DMCreateSubDM(pdm, 1, &fieldU, &isU, &dmU)); 292 PetscCall(DMCreateSubDM(pdm, 1, &fieldA, &isA, &dmA)); 293 PetscCall(DMCreateSubDM(pdm, 1, &fieldS, &isS, &dmS)); 294 PetscCall(DMCreateSubDM(pdm, 2, fieldUA, &isUA, &dmUA)); 295 296 PetscCall(PetscMalloc1(2, &dmList)); 297 dmList[0] = dmU; 298 dmList[1] = dmA; 299 300 PetscCall(DMCreateSuperDM(dmList, 2, NULL, &dmUA2)); 301 PetscCall(PetscFree(dmList)); 302 303 PetscCall(DMGetGlobalVector(pdm, &X)); 304 PetscCall(DMGetGlobalVector(dmU, &U)); 305 PetscCall(DMGetGlobalVector(dmA, &A)); 306 PetscCall(DMGetGlobalVector(dmS, &S)); 307 PetscCall(DMGetGlobalVector(dmUA, &UA)); 308 PetscCall(DMGetGlobalVector(dmUA2, &UA2)); 309 310 PetscCall(PetscObjectSetName((PetscObject)U, "U")); 311 PetscCall(PetscObjectSetName((PetscObject)A, "Alpha")); 312 PetscCall(PetscObjectSetName((PetscObject)S, "Sigma")); 313 PetscCall(PetscObjectSetName((PetscObject)UA, "UAlpha")); 314 PetscCall(PetscObjectSetName((PetscObject)UA2, "UAlpha2")); 315 PetscCall(VecSet(X, -111.)); 316 317 /* 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 */ 318 { 319 PetscSection sectionUA; 320 Vec UALoc; 321 PetscSection coordSection; 322 Vec coord; 323 PetscScalar *cval, *xyz; 324 PetscInt clSize, i, j; 325 326 PetscCall(DMGetLocalSection(dmUA, §ionUA)); 327 PetscCall(DMGetLocalVector(dmUA, &UALoc)); 328 PetscCall(VecGetArray(UALoc, &cval)); 329 PetscCall(DMGetCoordinateSection(dmUA, &coordSection)); 330 PetscCall(DMGetCoordinatesLocal(dmUA, &coord)); 331 PetscCall(DMPlexGetChart(dmUA, &pStart, &pEnd)); 332 333 for (p = pStart; p < pEnd; ++p) { 334 PetscInt dofUA, offUA; 335 336 PetscCall(PetscSectionGetDof(sectionUA, p, &dofUA)); 337 if (dofUA > 0) { 338 xyz = NULL; 339 PetscCall(PetscSectionGetOffset(sectionUA, p, &offUA)); 340 PetscCall(DMPlexVecGetClosure(dmUA, coordSection, coord, p, &clSize, &xyz)); 341 cval[offUA + sdim] = 0.; 342 for (i = 0; i < sdim; ++i) { 343 cval[offUA + i] = 0; 344 for (j = 0; j < clSize / sdim; ++j) cval[offUA + i] += xyz[j * sdim + i]; 345 cval[offUA + i] = cval[offUA + i] * sdim / clSize; 346 cval[offUA + sdim] += PetscSqr(cval[offUA + i]); 347 } 348 PetscCall(DMPlexVecRestoreClosure(dmUA, coordSection, coord, p, &clSize, &xyz)); 349 } 350 } 351 PetscCall(VecRestoreArray(UALoc, &cval)); 352 PetscCall(DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA)); 353 PetscCall(DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA)); 354 PetscCall(DMRestoreLocalVector(dmUA, &UALoc)); 355 356 /* Update X */ 357 PetscCall(VecISCopy(X, isUA, SCATTER_FORWARD, UA)); 358 PetscCall(VecViewFromOptions(UA, NULL, "-ua_vec_view")); 359 360 /* Restrict to U and Alpha */ 361 PetscCall(VecISCopy(X, isU, SCATTER_REVERSE, U)); 362 PetscCall(VecISCopy(X, isA, SCATTER_REVERSE, A)); 363 364 /* restrict to UA2 */ 365 PetscCall(VecISCopy(X, isUA, SCATTER_REVERSE, UA2)); 366 PetscCall(VecViewFromOptions(UA2, NULL, "-ua2_vec_view")); 367 } 368 369 { 370 Vec tmpVec; 371 PetscSection coordSection; 372 Vec coord; 373 PetscReal norm; 374 PetscReal time = 1.234; 375 376 /* Writing nodal variables to ExodusII file */ 377 PetscCall(DMSetOutputSequenceNumber(dmU, 0, time)); 378 PetscCall(DMSetOutputSequenceNumber(dmA, 0, time)); 379 380 PetscCall(VecView(U, viewer)); 381 PetscCall(VecView(A, viewer)); 382 383 /* Saving U and Alpha in one shot. 384 For this, we need to cheat and change the Vec's name 385 Note that in the end we write variables one component at a time, 386 so that there is no real values in doing this 387 */ 388 389 PetscCall(DMSetOutputSequenceNumber(dmUA, 1, time)); 390 PetscCall(DMGetGlobalVector(dmUA, &tmpVec)); 391 PetscCall(VecCopy(UA, tmpVec)); 392 PetscCall(PetscObjectSetName((PetscObject)tmpVec, "U")); 393 PetscCall(VecView(tmpVec, viewer)); 394 /* Reading nodal variables in Exodus file */ 395 PetscCall(VecSet(tmpVec, -1000.0)); 396 PetscCall(VecLoad(tmpVec, viewer)); 397 PetscCall(VecAXPY(UA, -1.0, tmpVec)); 398 PetscCall(VecNorm(UA, NORM_INFINITY, &norm)); 399 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "UAlpha ||Vin - Vout|| = %g", (double)norm); 400 PetscCall(DMRestoreGlobalVector(dmUA, &tmpVec)); 401 402 /* same thing with the UA2 Vec obtained from the superDM */ 403 PetscCall(DMGetGlobalVector(dmUA2, &tmpVec)); 404 PetscCall(VecCopy(UA2, tmpVec)); 405 PetscCall(PetscObjectSetName((PetscObject)tmpVec, "U")); 406 PetscCall(DMSetOutputSequenceNumber(dmUA2, 2, time)); 407 PetscCall(VecView(tmpVec, viewer)); 408 /* Reading nodal variables in Exodus file */ 409 PetscCall(VecSet(tmpVec, -1000.0)); 410 PetscCall(VecLoad(tmpVec, viewer)); 411 PetscCall(VecAXPY(UA2, -1.0, tmpVec)); 412 PetscCall(VecNorm(UA2, NORM_INFINITY, &norm)); 413 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "UAlpha2 ||Vin - Vout|| = %g", (double)norm); 414 PetscCall(DMRestoreGlobalVector(dmUA2, &tmpVec)); 415 416 /* Building and saving Sigma 417 We set sigma_0 = rank (to see partitioning) 418 sigma_1 = cell set ID 419 sigma_2 = x_coordinate of the cell center of mass 420 */ 421 PetscCall(DMGetCoordinateSection(dmS, &coordSection)); 422 PetscCall(DMGetCoordinatesLocal(dmS, &coord)); 423 PetscCall(DMGetLabelIdIS(dmS, "Cell Sets", &csIS)); 424 PetscCall(DMGetLabelSize(dmS, "Cell Sets", &numCS)); 425 PetscCall(ISGetIndices(csIS, &csID)); 426 for (set = 0; set < numCS; ++set) { 427 /* We know that all cells in a cell set have the same type, so we can dimension cval and xyz once for each cell set */ 428 IS cellIS; 429 const PetscInt *cellID; 430 PetscInt numCells, cell; 431 PetscScalar *cval = NULL, *xyz = NULL; 432 PetscInt clSize, cdimCoord, c; 433 434 PetscCall(DMGetStratumIS(dmS, "Cell Sets", csID[set], &cellIS)); 435 PetscCall(ISGetIndices(cellIS, &cellID)); 436 PetscCall(ISGetSize(cellIS, &numCells)); 437 for (cell = 0; cell < numCells; cell++) { 438 PetscCall(DMPlexVecGetClosure(dmS, NULL, S, cellID[cell], &clSize, &cval)); 439 PetscCall(DMPlexVecGetClosure(dmS, coordSection, coord, cellID[cell], &cdimCoord, &xyz)); 440 cval[0] = rank; 441 cval[1] = csID[set]; 442 cval[2] = 0.; 443 for (c = 0; c < cdimCoord / sdim; c++) cval[2] += xyz[c * sdim]; 444 cval[2] = cval[2] * sdim / cdimCoord; 445 PetscCall(DMPlexVecSetClosure(dmS, NULL, S, cellID[cell], cval, INSERT_ALL_VALUES)); 446 } 447 PetscCall(DMPlexVecRestoreClosure(dmS, NULL, S, cellID[0], &clSize, &cval)); 448 PetscCall(DMPlexVecRestoreClosure(dmS, coordSection, coord, cellID[0], NULL, &xyz)); 449 PetscCall(ISRestoreIndices(cellIS, &cellID)); 450 PetscCall(ISDestroy(&cellIS)); 451 } 452 PetscCall(ISRestoreIndices(csIS, &csID)); 453 PetscCall(ISDestroy(&csIS)); 454 PetscCall(VecViewFromOptions(S, NULL, "-s_vec_view")); 455 456 /* Writing zonal variables in Exodus file */ 457 PetscCall(DMSetOutputSequenceNumber(dmS, 0, time)); 458 PetscCall(VecView(S, viewer)); 459 460 /* Reading zonal variables in Exodus file */ 461 PetscCall(DMGetGlobalVector(dmS, &tmpVec)); 462 PetscCall(VecSet(tmpVec, -1000.0)); 463 PetscCall(PetscObjectSetName((PetscObject)tmpVec, "Sigma")); 464 PetscCall(VecLoad(tmpVec, viewer)); 465 PetscCall(VecAXPY(S, -1.0, tmpVec)); 466 PetscCall(VecNorm(S, NORM_INFINITY, &norm)); 467 PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Sigma ||Vin - Vout|| = %g", (double)norm); 468 PetscCall(DMRestoreGlobalVector(dmS, &tmpVec)); 469 } 470 PetscCall(PetscViewerDestroy(&viewer)); 471 472 PetscCall(DMRestoreGlobalVector(dmUA2, &UA2)); 473 PetscCall(DMRestoreGlobalVector(dmUA, &UA)); 474 PetscCall(DMRestoreGlobalVector(dmS, &S)); 475 PetscCall(DMRestoreGlobalVector(dmA, &A)); 476 PetscCall(DMRestoreGlobalVector(dmU, &U)); 477 PetscCall(DMRestoreGlobalVector(pdm, &X)); 478 PetscCall(DMDestroy(&dmU)); 479 PetscCall(ISDestroy(&isU)); 480 PetscCall(DMDestroy(&dmA)); 481 PetscCall(ISDestroy(&isA)); 482 PetscCall(DMDestroy(&dmS)); 483 PetscCall(ISDestroy(&isS)); 484 PetscCall(DMDestroy(&dmUA)); 485 PetscCall(ISDestroy(&isUA)); 486 PetscCall(DMDestroy(&dmUA2)); 487 PetscCall(DMDestroy(&pdm)); 488 if (size > 1) PetscCall(DMDestroy(&dm)); 489 PetscCall(PetscFree2(pStartDepth, pEndDepth)); 490 PetscCall(PetscFinalize()); 491 return 0; 492 } 493 494 /*TEST 495 496 build: 497 requires: exodusii pnetcdf !complex 498 # 2D seq 499 test: 500 suffix: 0 501 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 1 502 #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 503 test: 504 suffix: 1 505 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 1 506 test: 507 suffix: 2 508 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 1 509 #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 510 test: 511 suffix: 3 512 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 2 513 test: 514 suffix: 4 515 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 2 516 test: 517 suffix: 5 518 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 2 519 520 # 2D par 521 test: 522 suffix: 6 523 nsize: 2 524 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 1 525 #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 526 test: 527 suffix: 7 528 nsize: 2 529 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 1 530 test: 531 suffix: 8 532 nsize: 2 533 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 1 534 #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: invalid dimension ID or name 535 test: 536 suffix: 9 537 nsize: 2 538 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 2 539 test: 540 suffix: 10 541 nsize: 2 542 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 2 543 test: 544 # Something is now broken with parallel read/write for wahtever shape H is 545 TODO: broken 546 suffix: 11 547 nsize: 2 548 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 2 549 550 #3d seq 551 test: 552 suffix: 12 553 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 1 554 test: 555 suffix: 13 556 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 1 557 #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 558 test: 559 suffix: 14 560 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 2 561 test: 562 suffix: 15 563 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 2 564 #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 565 #3d par 566 test: 567 suffix: 16 568 nsize: 2 569 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 1 570 test: 571 suffix: 17 572 nsize: 2 573 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 1 574 #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 575 test: 576 suffix: 18 577 nsize: 2 578 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 2 579 test: 580 suffix: 19 581 nsize: 2 582 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 2 583 #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 584 585 TEST*/ 586