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 PetscErrorCode ierr; 32 33 ierr = PetscInitialize(&argc, &argv,NULL, help);if (ierr) return ierr; 34 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr); 35 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr); 36 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FEM Layout Options", "ex26");CHKERRQ(ierr); 37 ierr = PetscOptionsString("-i", "Filename to read", "ex26", ifilename, ifilename, sizeof(ifilename), NULL);CHKERRQ(ierr); 38 ierr = PetscOptionsString("-o", "Filename to write", "ex26", ofilename, ofilename, sizeof(ofilename), NULL);CHKERRQ(ierr); 39 ierr = PetscOptionsBoundedInt("-order", "FEM polynomial order", "ex26", order, &order, NULL,1);CHKERRQ(ierr); 40 ierr = PetscOptionsEnd();CHKERRQ(ierr); 41 if ((order > 2) || (order < 1)) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported polynomial order %D not in [1, 2]", order); 42 43 ex_opts(EX_VERBOSE+EX_DEBUG); 44 ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename, PETSC_TRUE, &dm);CHKERRQ(ierr); 45 ierr = DMSetFromOptions(dm);CHKERRQ(ierr); 46 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr); 47 48 /* Create the main section containning all fields */ 49 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);CHKERRQ(ierr); 50 ierr = PetscSectionSetNumFields(section, 3);CHKERRQ(ierr); 51 ierr = PetscSectionSetFieldName(section, fieldU, "U");CHKERRQ(ierr); 52 ierr = PetscSectionSetFieldName(section, fieldA, "Alpha");CHKERRQ(ierr); 53 ierr = PetscSectionSetFieldName(section, fieldS, "Sigma");CHKERRQ(ierr); 54 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 55 ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr); 56 ierr = DMGetDimension(dm, &sdim);CHKERRQ(ierr); 57 ierr = PetscMalloc2(sdim+1, &pStartDepth, sdim+1, &pEndDepth);CHKERRQ(ierr); 58 for (d = 0; d <= sdim; ++d) {ierr = DMPlexGetDepthStratum(dm, d, &pStartDepth[d], &pEndDepth[d]);CHKERRQ(ierr);} 59 /* Vector field U, Scalar field Alpha, Tensor field Sigma */ 60 ierr = PetscSectionSetFieldComponents(section, fieldU, sdim);CHKERRQ(ierr); 61 ierr = PetscSectionSetFieldComponents(section, fieldA, 1);CHKERRQ(ierr); 62 ierr = PetscSectionSetFieldComponents(section, fieldS, sdim*(sdim+1)/2);CHKERRQ(ierr); 63 64 /* Going through cell sets then cells, and setting up storage for the sections */ 65 ierr = DMGetLabelSize(dm, "Cell Sets", &numCS); 66 ierr = DMGetLabelIdIS(dm, "Cell Sets", &csIS);CHKERRQ(ierr); 67 if (csIS) {ierr = ISGetIndices(csIS, &csID);CHKERRQ(ierr);} 68 for (set = 0; set < numCS; set++) { 69 IS cellIS; 70 const PetscInt *cellID; 71 PetscInt numCells, cell, closureSize, *closureA = NULL; 72 73 ierr = DMGetStratumSize(dm, "Cell Sets", csID[set], &numCells);CHKERRQ(ierr); 74 ierr = DMGetStratumIS(dm, "Cell Sets", csID[set], &cellIS);CHKERRQ(ierr); 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: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No layout for dimension %D", 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 ierr = ISGetIndices(cellIS, &cellID);CHKERRQ(ierr); 106 ierr = DMPlexGetTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA);CHKERRQ(ierr); 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: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Unknown element with closure size %D\n", closureSize); 145 } 146 ierr = DMPlexRestoreTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA);CHKERRQ(ierr); 147 148 for (cell = 0; cell < numCells; cell++) { 149 PetscInt *closure = NULL; 150 151 ierr = DMPlexGetTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 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 ierr = PetscSectionSetDof(section, closure[2*p], dofU[d]+dofA[d]+dofS[d]);CHKERRQ(ierr); 157 ierr = PetscSectionSetFieldDof(section, closure[2*p], fieldU, dofU[d]);CHKERRQ(ierr); 158 ierr = PetscSectionSetFieldDof(section, closure[2*p], fieldA, dofA[d]);CHKERRQ(ierr); 159 ierr = PetscSectionSetFieldDof(section, closure[2*p], fieldS, dofS[d]);CHKERRQ(ierr); 160 } 161 } 162 } 163 ierr = DMPlexRestoreTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 164 } 165 ierr = ISRestoreIndices(cellIS, &cellID);CHKERRQ(ierr); 166 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 167 } 168 } 169 if (csIS) {ierr = ISRestoreIndices(csIS, &csID);CHKERRQ(ierr);} 170 ierr = ISDestroy(&csIS);CHKERRQ(ierr); 171 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 172 ierr = DMSetLocalSection(dm, section);CHKERRQ(ierr); 173 ierr = PetscObjectViewFromOptions((PetscObject) section, NULL, "-dm_section_view");CHKERRQ(ierr); 174 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 175 176 { 177 /* TODO: Replace with ExodusII viewer */ 178 /* Create the exodus result file */ 179 PetscInt numstep = 3, step; 180 char *nodalVarName[4]; 181 char *zonalVarName[6]; 182 int *truthtable; 183 PetscInt numNodalVar, numZonalVar, i; 184 int CPU_word_size, IO_word_size, EXO_mode; 185 186 ex_opts(EX_VERBOSE+EX_DEBUG); 187 if (!rank) { 188 CPU_word_size = sizeof(PetscReal); 189 IO_word_size = sizeof(PetscReal); 190 EXO_mode = EX_CLOBBER; 191 #if defined(PETSC_USE_64BIT_INDICES) 192 EXO_mode += EX_ALL_INT64_API; 193 #endif 194 exoid = ex_create(ofilename, EXO_mode, &CPU_word_size, &IO_word_size); 195 if (exoid < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to open exodus file %\n", ofilename); 196 } 197 ierr = DMPlexView_ExodusII_Internal(dm, exoid, order);CHKERRQ(ierr); 198 199 if (!rank) { 200 /* "Format" the exodus result file, i.e. allocate space for nodal and zonal variables */ 201 switch (sdim) { 202 case 2: 203 numNodalVar = 3; 204 nodalVarName[0] = (char *) "U_x"; 205 nodalVarName[1] = (char *) "U_y"; 206 nodalVarName[2] = (char *) "Alpha"; 207 numZonalVar = 3; 208 zonalVarName[0] = (char *) "Sigma_11"; 209 zonalVarName[1] = (char *) "Sigma_22"; 210 zonalVarName[2] = (char *) "Sigma_12"; 211 break; 212 case 3: 213 numNodalVar = 4; 214 nodalVarName[0] = (char *) "U_x"; 215 nodalVarName[1] = (char *) "U_y"; 216 nodalVarName[2] = (char *) "U_z"; 217 nodalVarName[3] = (char *) "Alpha"; 218 numZonalVar = 6; 219 zonalVarName[0] = (char *) "Sigma_11"; 220 zonalVarName[1] = (char *) "Sigma_22"; 221 zonalVarName[2] = (char *) "Sigma_33"; 222 zonalVarName[3] = (char *) "Sigma_23"; 223 zonalVarName[4] = (char *) "Sigma_13"; 224 zonalVarName[5] = (char *) "Sigma_12"; 225 break; 226 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No layout for dimension %D", sdim); 227 } 228 ierr = ex_put_variable_param(exoid, EX_ELEM_BLOCK, numZonalVar);CHKERRQ(ierr); 229 ierr = ex_put_variable_names(exoid, EX_ELEM_BLOCK, numZonalVar, zonalVarName);CHKERRQ(ierr); 230 ierr = ex_put_variable_param(exoid, EX_NODAL, numNodalVar);CHKERRQ(ierr); 231 ierr = ex_put_variable_names(exoid, EX_NODAL, numNodalVar, nodalVarName);CHKERRQ(ierr); 232 numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK); 233 ierr = PetscMalloc1(numZonalVar * numCS, &truthtable);CHKERRQ(ierr); 234 for (i = 0; i < numZonalVar * numCS; ++i) truthtable[i] = 1; 235 ierr = ex_put_truth_table(exoid, EX_ELEM_BLOCK, numCS, numZonalVar, truthtable);CHKERRQ(ierr); 236 ierr = PetscFree(truthtable);CHKERRQ(ierr); 237 /* Writing time step information in the file. Note that this is currently broken in the exodus library for netcdf4 (HDF5-based) files */ 238 for (step = 0; step < numstep; ++step) { 239 PetscReal time = step; 240 ierr = ex_put_time(exoid, step+1, &time);CHKERRQ(ierr); 241 } 242 ierr = ex_close(exoid);CHKERRQ(ierr); 243 } 244 } 245 246 { 247 DM pdm; 248 PetscSF migrationSF; 249 PetscInt ovlp = 0; 250 PetscPartitioner part; 251 252 ierr = DMSetUseNatural(dm,PETSC_TRUE);CHKERRQ(ierr); 253 ierr = DMPlexGetPartitioner(dm,&part);CHKERRQ(ierr); 254 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 255 ierr = DMPlexDistribute(dm,ovlp,&migrationSF,&pdm);CHKERRQ(ierr); 256 if (pdm) { 257 ierr = DMPlexSetMigrationSF(pdm,migrationSF);CHKERRQ(ierr); 258 ierr = PetscSFDestroy(&migrationSF);CHKERRQ(ierr); 259 ierr = DMDestroy(&dm);CHKERRQ(ierr); 260 dm = pdm; 261 ierr = DMViewFromOptions(dm,NULL,"-dm_view");CHKERRQ(ierr); 262 } 263 } 264 265 { 266 /* TODO Replace with ExodusII viewer */ 267 /* Reopen the exodus result file on all processors */ 268 MPI_Info mpi_info = MPI_INFO_NULL; 269 int CPU_word_size, IO_word_size, EXO_mode; 270 float EXO_version; 271 272 EXO_mode = EX_WRITE; 273 CPU_word_size = sizeof(PetscReal); 274 IO_word_size = sizeof(PetscReal); 275 #if defined(PETSC_USE_64BIT_INDICES) 276 EXO_mode += EX_ALL_INT64_API; 277 #endif 278 exoid = ex_open_par(ofilename, EXO_mode, &CPU_word_size, &IO_word_size, &EXO_version, PetscObjectComm((PetscObject) dm), mpi_info); 279 } 280 281 /* Get DM and IS for each field of dm */ 282 ierr = DMCreateSubDM(dm, 1, &fieldU, &isU, &dmU);CHKERRQ(ierr); 283 ierr = DMCreateSubDM(dm, 1, &fieldA, &isA, &dmA);CHKERRQ(ierr); 284 ierr = DMCreateSubDM(dm, 1, &fieldS, &isS, &dmS);CHKERRQ(ierr); 285 ierr = DMCreateSubDM(dm, 2, fieldUA, &isUA, &dmUA);CHKERRQ(ierr); 286 287 ierr = PetscMalloc1(2,&dmList);CHKERRQ(ierr); 288 dmList[0] = dmU; 289 dmList[1] = dmA; 290 /* We temporarily disable dmU->useNatural to test that we can reconstruct the 291 NaturaltoGlobal SF from any of the dm in dms 292 */ 293 dmU->useNatural = PETSC_FALSE; 294 ierr = DMCreateSuperDM(dmList,2,NULL,&dmUA2);CHKERRQ(ierr); 295 dmU->useNatural = PETSC_TRUE; 296 ierr = PetscFree(dmList);CHKERRQ(ierr); 297 298 ierr = DMGetGlobalVector(dm, &X);CHKERRQ(ierr); 299 ierr = DMGetGlobalVector(dmU, &U);CHKERRQ(ierr); 300 ierr = DMGetGlobalVector(dmA, &A);CHKERRQ(ierr); 301 ierr = DMGetGlobalVector(dmS, &S);CHKERRQ(ierr); 302 ierr = DMGetGlobalVector(dmUA, &UA);CHKERRQ(ierr); 303 ierr = DMGetGlobalVector(dmUA2, &UA2);CHKERRQ(ierr); 304 305 ierr = PetscObjectSetName((PetscObject) U, "U");CHKERRQ(ierr); 306 ierr = PetscObjectSetName((PetscObject) A, "Alpha");CHKERRQ(ierr); 307 ierr = PetscObjectSetName((PetscObject) S, "Sigma");CHKERRQ(ierr); 308 ierr = PetscObjectSetName((PetscObject) UA, "UAlpha");CHKERRQ(ierr); 309 ierr = PetscObjectSetName((PetscObject) UA2, "UAlpha2");CHKERRQ(ierr); 310 ierr = VecSet(X, -111.);CHKERRQ(ierr); 311 312 /* 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 */ 313 { 314 PetscSection sectionUA; 315 Vec UALoc; 316 PetscSection coordSection; 317 Vec coord; 318 PetscScalar *cval, *xyz; 319 PetscInt cdimCoord = 24; 320 321 ierr = DMGetLocalSection(dmUA, §ionUA);CHKERRQ(ierr); 322 ierr = DMGetLocalVector(dmUA, &UALoc);CHKERRQ(ierr); 323 ierr = VecGetArray(UALoc, &cval);CHKERRQ(ierr); 324 ierr = DMGetCoordinateSection(dmUA, &coordSection);CHKERRQ(ierr); 325 ierr = DMGetCoordinatesLocal(dmUA, &coord);CHKERRQ(ierr); 326 ierr = DMPlexGetChart(dmUA, &pStart, &pEnd);CHKERRQ(ierr); 327 ierr = PetscMalloc1(cdimCoord, &xyz);CHKERRQ(ierr); 328 /* I know that UA has 0 or sdim+1 dof at each point, since both U and Alpha use the same discretization 329 The maximum possible size for the closure of the coordinate section is 8*3 at the cell center */ 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 PetscInt clSize = cdimCoord, i, j; 336 337 ierr = PetscSectionGetOffset(sectionUA, p, &offUA);CHKERRQ(ierr); 338 ierr = DMPlexVecGetClosure(dmUA, coordSection, coord, p, &clSize, &xyz);CHKERRQ(ierr); 339 cval[offUA+sdim] = 0.; 340 for (i = 0; i < sdim; ++i) { 341 cval[offUA+i] = 0; 342 for (j = 0; j < clSize/sdim; ++j) { 343 cval[offUA+i] += xyz[j*sdim+i]; 344 } 345 cval[offUA+i] = cval[offUA+i] * sdim / clSize; 346 cval[offUA+sdim] += PetscSqr(cval[offUA+i]); 347 } 348 } 349 } 350 ierr = PetscFree(xyz);CHKERRQ(ierr); 351 ierr = VecRestoreArray(UALoc, &cval);CHKERRQ(ierr); 352 ierr = DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA);CHKERRQ(ierr); 353 ierr = DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA);CHKERRQ(ierr); 354 ierr = DMRestoreLocalVector(dmUA, &UALoc);CHKERRQ(ierr); 355 /* Update X */ 356 ierr = VecISCopy(X, isUA, SCATTER_FORWARD, UA);CHKERRQ(ierr); 357 ierr = VecViewFromOptions(UA, NULL, "-ua_vec_view");CHKERRQ(ierr); 358 /* Restrict to U and Alpha */ 359 ierr = VecISCopy(X, isU, SCATTER_REVERSE, U);CHKERRQ(ierr); 360 ierr = VecISCopy(X, isA, SCATTER_REVERSE, A);CHKERRQ(ierr); 361 /* restrict to UA2 */ 362 ierr = VecISCopy(X, isUA, SCATTER_REVERSE, UA2);CHKERRQ(ierr); 363 ierr = VecViewFromOptions(UA2, NULL, "-ua2_vec_view");CHKERRQ(ierr); 364 } 365 366 { 367 Vec tmpVec; 368 PetscSection coordSection; 369 Vec coord; 370 PetscReal norm; 371 372 /* Writing nodal variables to ExodusII file */ 373 ierr = VecViewPlex_ExodusII_Nodal_Internal(U, exoid, 1);CHKERRQ(ierr); 374 ierr = VecViewPlex_ExodusII_Nodal_Internal(A, exoid, 1);CHKERRQ(ierr); 375 /* Saving U and Alpha in one shot. 376 For this, we need to cheat and change the Vec's name 377 Note that in the end we write variables one component at a time, so that there is no real values in doing this */ 378 ierr = DMGetGlobalVector(dmUA, &tmpVec);CHKERRQ(ierr); 379 ierr = VecCopy(UA, tmpVec);CHKERRQ(ierr); 380 ierr = PetscObjectSetName((PetscObject) tmpVec, "U");CHKERRQ(ierr); 381 ierr = VecViewPlex_ExodusII_Nodal_Internal(tmpVec, exoid, 2);CHKERRQ(ierr); 382 /* Reading nodal variables in Exodus file */ 383 ierr = VecSet(tmpVec, -1000.0);CHKERRQ(ierr); 384 ierr = VecLoadPlex_ExodusII_Nodal_Internal(tmpVec, exoid, 2);CHKERRQ(ierr); 385 ierr = VecAXPY(UA, -1.0, tmpVec);CHKERRQ(ierr); 386 ierr = VecNorm(UA, NORM_INFINITY, &norm);CHKERRQ(ierr); 387 if (norm > PETSC_SQRT_MACHINE_EPSILON) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "UAlpha ||Vin - Vout|| = %g\n", (double) norm); 388 ierr = DMRestoreGlobalVector(dmUA, &tmpVec);CHKERRQ(ierr); 389 390 /* same thing with the UA2 Vec obtained from the superDM */ 391 ierr = DMGetGlobalVector(dmUA2, &tmpVec);CHKERRQ(ierr); 392 ierr = VecCopy(UA2, tmpVec);CHKERRQ(ierr); 393 ierr = PetscObjectSetName((PetscObject) tmpVec, "U");CHKERRQ(ierr); 394 ierr = VecViewPlex_ExodusII_Nodal_Internal(tmpVec, exoid, 3);CHKERRQ(ierr); 395 /* Reading nodal variables in Exodus file */ 396 ierr = VecSet(tmpVec, -1000.0);CHKERRQ(ierr); 397 ierr = VecLoadPlex_ExodusII_Nodal_Internal(tmpVec, exoid, 3);CHKERRQ(ierr); 398 ierr = VecAXPY(UA2, -1.0, tmpVec);CHKERRQ(ierr); 399 ierr = VecNorm(UA2, NORM_INFINITY, &norm);CHKERRQ(ierr); 400 if (norm > PETSC_SQRT_MACHINE_EPSILON) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "UAlpha2 ||Vin - Vout|| = %g\n", (double) norm); 401 ierr = DMRestoreGlobalVector(dmUA2, &tmpVec);CHKERRQ(ierr); 402 403 /* Building and saving Sigma 404 We set sigma_0 = rank (to see partitioning) 405 sigma_1 = cell set ID 406 sigma_2 = x_coordinate of the cell center of mass */ 407 ierr = DMGetCoordinateSection(dmS, &coordSection);CHKERRQ(ierr); 408 ierr = DMGetCoordinatesLocal(dmS, &coord);CHKERRQ(ierr); 409 ierr = DMGetLabelIdIS(dmS, "Cell Sets", &csIS);CHKERRQ(ierr); 410 ierr = DMGetLabelSize(dmS, "Cell Sets", &numCS); 411 ierr = ISGetIndices(csIS, &csID);CHKERRQ(ierr); 412 for (set = 0; set < numCS; ++set) { 413 /* 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 */ 414 IS cellIS; 415 const PetscInt *cellID; 416 PetscInt numCells, cell; 417 PetscScalar *cval = NULL, *xyz = NULL; 418 PetscInt clSize, cdimCoord, c; 419 420 ierr = DMGetStratumIS(dmS, "Cell Sets", csID[set], &cellIS);CHKERRQ(ierr); 421 ierr = ISGetIndices(cellIS, &cellID);CHKERRQ(ierr); 422 ierr = ISGetSize(cellIS, &numCells);CHKERRQ(ierr); 423 for (cell = 0; cell < numCells; cell++) { 424 ierr = DMPlexVecGetClosure(dmS, NULL, S, cellID[cell], &clSize, &cval);CHKERRQ(ierr); 425 ierr = DMPlexVecGetClosure(dmS, coordSection, coord, cellID[cell], &cdimCoord, &xyz);CHKERRQ(ierr); 426 cval[0] = rank; 427 cval[1] = csID[set]; 428 cval[2] = 0.; 429 for (c = 0; c < cdimCoord/sdim; c++) cval[2] += xyz[c*sdim]; 430 cval[2] = cval[2] * sdim / cdimCoord; 431 ierr = DMPlexVecSetClosure(dmS, NULL, S, cellID[cell], cval, INSERT_ALL_VALUES);CHKERRQ(ierr); 432 } 433 ierr = DMPlexVecRestoreClosure(dmS, NULL, S, cellID[0], &clSize, &cval);CHKERRQ(ierr); 434 ierr = DMPlexVecRestoreClosure(dmS, coordSection, coord, cellID[0], NULL, &xyz);CHKERRQ(ierr); 435 ierr = ISRestoreIndices(cellIS, &cellID);CHKERRQ(ierr); 436 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 437 } 438 ierr = ISRestoreIndices(csIS, &csID);CHKERRQ(ierr); 439 ierr = ISDestroy(&csIS);CHKERRQ(ierr); 440 ierr = VecViewFromOptions(S, NULL, "-s_vec_view");CHKERRQ(ierr); 441 /* Writing zonal variables in Exodus file */ 442 ierr = VecViewPlex_ExodusII_Zonal_Internal(S, exoid, 1);CHKERRQ(ierr); 443 /* Reading zonal variables in Exodus file */ 444 ierr = DMGetGlobalVector(dmS, &tmpVec);CHKERRQ(ierr); 445 ierr = VecSet(tmpVec, -1000.0);CHKERRQ(ierr); 446 ierr = PetscObjectSetName((PetscObject) tmpVec, "Sigma");CHKERRQ(ierr); 447 ierr = VecLoadPlex_ExodusII_Zonal_Internal(tmpVec, exoid, 1);CHKERRQ(ierr); 448 ierr = VecAXPY(S, -1.0, tmpVec);CHKERRQ(ierr); 449 ierr = VecNorm(S, NORM_INFINITY, &norm); 450 if (norm > PETSC_SQRT_MACHINE_EPSILON) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Sigma ||Vin - Vout|| = %g\n", (double) norm); 451 ierr = DMRestoreGlobalVector(dmS, &tmpVec);CHKERRQ(ierr); 452 } 453 ierr = ex_close(exoid);CHKERRQ(ierr); 454 455 ierr = DMRestoreGlobalVector(dmUA2, &UA2);CHKERRQ(ierr); 456 ierr = DMRestoreGlobalVector(dmUA, &UA);CHKERRQ(ierr); 457 ierr = DMRestoreGlobalVector(dmS, &S);CHKERRQ(ierr); 458 ierr = DMRestoreGlobalVector(dmA, &A);CHKERRQ(ierr); 459 ierr = DMRestoreGlobalVector(dmU, &U);CHKERRQ(ierr); 460 ierr = DMRestoreGlobalVector(dm, &X);CHKERRQ(ierr); 461 ierr = DMDestroy(&dmU);CHKERRQ(ierr); ierr = ISDestroy(&isU);CHKERRQ(ierr); 462 ierr = DMDestroy(&dmA);CHKERRQ(ierr); ierr = ISDestroy(&isA);CHKERRQ(ierr); 463 ierr = DMDestroy(&dmS);CHKERRQ(ierr); ierr = ISDestroy(&isS);CHKERRQ(ierr); 464 ierr = DMDestroy(&dmUA);CHKERRQ(ierr);ierr = ISDestroy(&isUA);CHKERRQ(ierr); 465 ierr = DMDestroy(&dmUA2);CHKERRQ(ierr); 466 ierr = DMDestroy(&dm);CHKERRQ(ierr); 467 ierr = PetscFree2(pStartDepth, pEndDepth);CHKERRQ(ierr); 468 ierr = PetscFinalize(); 469 return ierr; 470 } 471 472 /*TEST 473 474 build: 475 requires: exodusii pnetcdf !complex 476 # 2D seq 477 test: 478 suffix: 0 479 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 480 #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 481 test: 482 suffix: 1 483 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 484 test: 485 suffix: 2 486 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 487 #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 488 test: 489 suffix: 3 490 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 491 test: 492 suffix: 4 493 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 494 test: 495 suffix: 5 496 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 497 498 # 2D par 499 test: 500 suffix: 6 501 nsize: 2 502 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 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: 7 506 nsize: 2 507 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 508 test: 509 suffix: 8 510 nsize: 2 511 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 512 #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: invalid dimension ID or name 513 test: 514 suffix: 9 515 nsize: 2 516 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 517 test: 518 suffix: 10 519 nsize: 2 520 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 521 test: 522 # Something is now broken with parallel read/write for wahtever shape H is 523 TODO: broken 524 suffix: 11 525 nsize: 2 526 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 527 528 #3d seq 529 test: 530 suffix: 12 531 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 532 test: 533 suffix: 13 534 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 535 #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 536 test: 537 suffix: 14 538 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 539 test: 540 suffix: 15 541 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 542 #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 543 #3d par 544 test: 545 suffix: 16 546 nsize: 2 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: 17 550 nsize: 2 551 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 552 #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 553 test: 554 suffix: 18 555 nsize: 2 556 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 557 test: 558 suffix: 19 559 nsize: 2 560 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 561 #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 562 563 TEST*/ 564