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