1 static char help[] = "Spectral element access patterns with Plex\n\n"; 2 3 #include <petscdmplex.h> 4 5 typedef struct { 6 PetscInt Nf; /* Number of fields */ 7 PetscInt *Nc; /* Number of components per field */ 8 PetscInt *k; /* Spectral order per field */ 9 } AppCtx; 10 11 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 12 { 13 PetscInt len; 14 PetscBool flg; 15 16 PetscFunctionBeginUser; 17 options->Nf = 0; 18 options->Nc = NULL; 19 options->k = NULL; 20 21 PetscOptionsBegin(comm, "", "SEM Problem Options", "DMPLEX"); 22 PetscCall(PetscOptionsBoundedInt("-num_fields", "The number of fields", "ex6.c", options->Nf, &options->Nf, NULL, 0)); 23 if (options->Nf) { 24 len = options->Nf; 25 PetscCall(PetscMalloc1(len, &options->Nc)); 26 PetscCall(PetscOptionsIntArray("-num_components", "The number of components per field", "ex6.c", options->Nc, &len, &flg)); 27 PetscCheck(!flg || !(len != options->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of components array is %" PetscInt_FMT " should be %" PetscInt_FMT, len, options->Nf); 28 len = options->Nf; 29 PetscCall(PetscMalloc1(len, &options->k)); 30 PetscCall(PetscOptionsIntArray("-order", "The spectral order per field", "ex6.c", options->k, &len, &flg)); 31 PetscCheck(!flg || !(len != options->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of order array is %" PetscInt_FMT " should be %" PetscInt_FMT, len, options->Nf); 32 } 33 PetscOptionsEnd(); 34 PetscFunctionReturn(PETSC_SUCCESS); 35 } 36 37 static PetscErrorCode LoadData2D(DM dm, PetscInt Ni, PetscInt Nj, PetscInt clSize, Vec u, AppCtx *user) 38 { 39 PetscInt i, j, f, c; 40 PetscScalar *closure; 41 42 PetscFunctionBeginUser; 43 PetscCall(PetscMalloc1(clSize, &closure)); 44 for (j = 0; j < Nj; ++j) { 45 for (i = 0; i < Ni; ++i) { 46 PetscInt ki, kj, o = 0; 47 PetscCall(PetscArrayzero(closure, clSize)); 48 49 for (f = 0; f < user->Nf; ++f) { 50 PetscInt ioff = i * user->k[f], joff = j * user->k[f]; 51 52 for (kj = 0; kj <= user->k[f]; ++kj) { 53 for (ki = 0; ki <= user->k[f]; ++ki) { 54 for (c = 0; c < user->Nc[f]; ++c) closure[o++] = ((kj + joff) * (Ni * user->k[f] + 1) + ki + ioff) * user->Nc[f] + c; 55 } 56 } 57 } 58 PetscCall(DMPlexVecSetClosure(dm, NULL, u, j * Ni + i, closure, INSERT_VALUES)); 59 } 60 } 61 PetscCall(PetscFree(closure)); 62 PetscFunctionReturn(PETSC_SUCCESS); 63 } 64 65 static PetscErrorCode LoadData3D(DM dm, PetscInt Ni, PetscInt Nj, PetscInt Nk, PetscInt clSize, Vec u, AppCtx *user) 66 { 67 PetscInt i, j, k, f, c; 68 PetscScalar *closure; 69 70 PetscFunctionBeginUser; 71 PetscCall(PetscMalloc1(clSize, &closure)); 72 for (k = 0; k < Nk; ++k) { 73 for (j = 0; j < Nj; ++j) { 74 for (i = 0; i < Ni; ++i) { 75 PetscInt ki, kj, kk, o = 0; 76 PetscCall(PetscArrayzero(closure, clSize)); 77 78 for (f = 0; f < user->Nf; ++f) { 79 PetscInt ioff = i * user->k[f], joff = j * user->k[f], koff = k * user->k[f]; 80 81 for (kk = 0; kk <= user->k[f]; ++kk) { 82 for (kj = 0; kj <= user->k[f]; ++kj) { 83 for (ki = 0; ki <= user->k[f]; ++ki) { 84 for (c = 0; c < user->Nc[f]; ++c) closure[o++] = (((kk + koff) * (Nj * user->k[f] + 1) + kj + joff) * (Ni * user->k[f] + 1) + ki + ioff) * user->Nc[f] + c; 85 } 86 } 87 } 88 } 89 PetscCall(DMPlexVecSetClosure(dm, NULL, u, (k * Nj + j) * Ni + i, closure, INSERT_VALUES)); 90 } 91 } 92 } 93 PetscCall(PetscFree(closure)); 94 PetscFunctionReturn(PETSC_SUCCESS); 95 } 96 97 static PetscErrorCode CheckPoint(DM dm, Vec u, PetscInt point, AppCtx *user) 98 { 99 PetscSection s; 100 PetscScalar *a; 101 const PetscScalar *array; 102 PetscInt dof, d; 103 104 PetscFunctionBeginUser; 105 PetscCall(DMGetLocalSection(dm, &s)); 106 PetscCall(VecGetArrayRead(u, &array)); 107 PetscCall(DMPlexPointLocalRead(dm, point, array, &a)); 108 PetscCall(PetscSectionGetDof(s, point, &dof)); 109 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Point %" PetscInt_FMT ": ", point)); 110 for (d = 0; d < dof; ++d) { 111 if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", ")); 112 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double)PetscRealPart(a[d]))); 113 } 114 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 115 PetscCall(VecRestoreArrayRead(u, &array)); 116 PetscFunctionReturn(PETSC_SUCCESS); 117 } 118 119 static PetscErrorCode ReadData2D(DM dm, Vec u, AppCtx *user) 120 { 121 PetscInt cStart, cEnd, cell; 122 123 PetscFunctionBeginUser; 124 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 125 for (cell = cStart; cell < cEnd; ++cell) { 126 PetscScalar *closure = NULL; 127 PetscInt closureSize, ki, kj, f, c, foff = 0; 128 129 PetscCall(DMPlexVecGetClosure(dm, NULL, u, cell, &closureSize, &closure)); 130 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT "\n", cell)); 131 for (f = 0; f < user->Nf; ++f) { 132 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Field %" PetscInt_FMT "\n", f)); 133 for (kj = user->k[f]; kj >= 0; --kj) { 134 for (ki = 0; ki <= user->k[f]; ++ki) { 135 if (ki > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 136 for (c = 0; c < user->Nc[f]; ++c) { 137 if (c > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ",")); 138 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double)PetscRealPart(closure[(kj * (user->k[f] + 1) + ki) * user->Nc[f] + c + foff]))); 139 } 140 } 141 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 142 } 143 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n\n")); 144 foff += PetscSqr(user->k[f] + 1); 145 } 146 PetscCall(DMPlexVecRestoreClosure(dm, NULL, u, cell, &closureSize, &closure)); 147 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n\n")); 148 } 149 PetscFunctionReturn(PETSC_SUCCESS); 150 } 151 152 static PetscErrorCode ReadData3D(DM dm, Vec u, AppCtx *user) 153 { 154 PetscInt cStart, cEnd, cell; 155 156 PetscFunctionBeginUser; 157 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 158 for (cell = cStart; cell < cEnd; ++cell) { 159 PetscScalar *closure = NULL; 160 PetscInt closureSize, ki, kj, kk, f, c, foff = 0; 161 162 PetscCall(DMPlexVecGetClosure(dm, NULL, u, cell, &closureSize, &closure)); 163 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT "\n", cell)); 164 for (f = 0; f < user->Nf; ++f) { 165 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Field %" PetscInt_FMT "\n", f)); 166 for (kk = user->k[f]; kk >= 0; --kk) { 167 for (kj = user->k[f]; kj >= 0; --kj) { 168 for (ki = 0; ki <= user->k[f]; ++ki) { 169 if (ki > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 170 for (c = 0; c < user->Nc[f]; ++c) { 171 if (c > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ",")); 172 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double)PetscRealPart(closure[((kk * (user->k[f] + 1) + kj) * (user->k[f] + 1) + ki) * user->Nc[f] + c + foff]))); 173 } 174 } 175 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 176 } 177 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 178 } 179 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n\n")); 180 foff += PetscSqr(user->k[f] + 1); 181 } 182 PetscCall(DMPlexVecRestoreClosure(dm, NULL, u, cell, &closureSize, &closure)); 183 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n\n")); 184 } 185 PetscFunctionReturn(PETSC_SUCCESS); 186 } 187 188 static PetscErrorCode SetSymmetries(DM dm, PetscSection s, AppCtx *user) 189 { 190 PetscInt dim, f, o, i, j, k, c, d; 191 DMLabel depthLabel; 192 193 PetscFunctionBegin; 194 PetscCall(DMGetDimension(dm, &dim)); 195 PetscCall(DMGetLabel(dm, "depth", &depthLabel)); 196 for (f = 0; f < user->Nf; f++) { 197 PetscSectionSym sym; 198 199 if (user->k[f] < 3) continue; /* No symmetries needed for order < 3, because no cell, facet, edge or vertex has more than one node */ 200 PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)s), depthLabel, &sym)); 201 202 for (d = 0; d <= dim; d++) { 203 if (d == 1) { 204 PetscInt numDof = user->k[f] - 1; 205 PetscInt numComp = user->Nc[f]; 206 PetscInt minOrnt = -1; 207 PetscInt maxOrnt = 1; 208 PetscInt **perms; 209 210 PetscCall(PetscCalloc1(maxOrnt - minOrnt, &perms)); 211 for (o = minOrnt; o < maxOrnt; o++) { 212 PetscInt *perm; 213 214 if (!o) { /* identity */ 215 perms[o - minOrnt] = NULL; 216 } else { 217 PetscCall(PetscMalloc1(numDof * numComp, &perm)); 218 for (i = numDof - 1, k = 0; i >= 0; i--) { 219 for (j = 0; j < numComp; j++, k++) perm[k] = i * numComp + j; 220 } 221 perms[o - minOrnt] = perm; 222 } 223 } 224 PetscCall(PetscSectionSymLabelSetStratum(sym, d, numDof * numComp, minOrnt, maxOrnt, PETSC_OWN_POINTER, (const PetscInt **)perms, NULL)); 225 } else if (d == 2) { 226 PetscInt perEdge = user->k[f] - 1; 227 PetscInt numDof = perEdge * perEdge; 228 PetscInt numComp = user->Nc[f]; 229 PetscInt minOrnt = -4; 230 PetscInt maxOrnt = 4; 231 PetscInt **perms; 232 233 PetscCall(PetscCalloc1(maxOrnt - minOrnt, &perms)); 234 for (o = minOrnt; o < maxOrnt; o++) { 235 PetscInt *perm; 236 237 if (!o) continue; /* identity */ 238 PetscCall(PetscMalloc1(numDof * numComp, &perm)); 239 /* We want to perm[k] to list which *localArray* position the *sectionArray* position k should go to for the given orientation*/ 240 switch (o) { 241 case 0: 242 break; /* identity */ 243 case -2: /* flip along (-1,-1)--( 1, 1), which swaps edges 0 and 3 and edges 1 and 2. This swaps the i and j variables */ 244 for (i = 0, k = 0; i < perEdge; i++) { 245 for (j = 0; j < perEdge; j++, k++) { 246 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * j + i) * numComp + c; 247 } 248 } 249 break; 250 case -1: /* flip along (-1, 0)--( 1, 0), which swaps edges 0 and 2. This reverses the i variable */ 251 for (i = 0, k = 0; i < perEdge; i++) { 252 for (j = 0; j < perEdge; j++, k++) { 253 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - i) + j) * numComp + c; 254 } 255 } 256 break; 257 case -4: /* flip along ( 1,-1)--(-1, 1), which swaps edges 0 and 1 and edges 2 and 3. This swaps the i and j variables and reverse both */ 258 for (i = 0, k = 0; i < perEdge; i++) { 259 for (j = 0; j < perEdge; j++, k++) { 260 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - j) + (perEdge - 1 - i)) * numComp + c; 261 } 262 } 263 break; 264 case -3: /* flip along ( 0,-1)--( 0, 1), which swaps edges 3 and 1. This reverses the j variable */ 265 for (i = 0, k = 0; i < perEdge; i++) { 266 for (j = 0; j < perEdge; j++, k++) { 267 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * i + (perEdge - 1 - j)) * numComp + c; 268 } 269 } 270 break; 271 case 1: /* rotate section edge 1 to local edge 0. This swaps the i and j variables and then reverses the j variable */ 272 for (i = 0, k = 0; i < perEdge; i++) { 273 for (j = 0; j < perEdge; j++, k++) { 274 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - j) + i) * numComp + c; 275 } 276 } 277 break; 278 case 2: /* rotate section edge 2 to local edge 0. This reverse both i and j variables */ 279 for (i = 0, k = 0; i < perEdge; i++) { 280 for (j = 0; j < perEdge; j++, k++) { 281 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - i) + (perEdge - 1 - j)) * numComp + c; 282 } 283 } 284 break; 285 case 3: /* rotate section edge 3 to local edge 0. This swaps the i and j variables and then reverses the i variable */ 286 for (i = 0, k = 0; i < perEdge; i++) { 287 for (j = 0; j < perEdge; j++, k++) { 288 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * j + (perEdge - 1 - i)) * numComp + c; 289 } 290 } 291 break; 292 default: 293 break; 294 } 295 perms[o - minOrnt] = perm; 296 } 297 PetscCall(PetscSectionSymLabelSetStratum(sym, d, numDof * numComp, minOrnt, maxOrnt, PETSC_OWN_POINTER, (const PetscInt **)perms, NULL)); 298 } 299 } 300 PetscCall(PetscSectionSetFieldSym(s, f, sym)); 301 PetscCall(PetscSectionSymDestroy(&sym)); 302 } 303 PetscCall(PetscSectionViewFromOptions(s, NULL, "-section_with_sym_view")); 304 PetscFunctionReturn(PETSC_SUCCESS); 305 } 306 307 int main(int argc, char **argv) 308 { 309 DM dm; 310 PetscSection s; 311 Vec u; 312 AppCtx user; 313 PetscInt dim, size = 0, f; 314 315 PetscFunctionBeginUser; 316 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 317 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 318 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 319 PetscCall(DMSetType(dm, DMPLEX)); 320 PetscCall(DMSetFromOptions(dm)); 321 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 322 PetscCall(DMGetDimension(dm, &dim)); 323 /* Create a section for SEM order k */ 324 { 325 PetscInt *numDof, d; 326 327 PetscCall(PetscMalloc1(user.Nf * (dim + 1), &numDof)); 328 for (f = 0; f < user.Nf; ++f) { 329 for (d = 0; d <= dim; ++d) numDof[f * (dim + 1) + d] = PetscPowInt(user.k[f] - 1, d) * user.Nc[f]; 330 size += PetscPowInt(user.k[f] + 1, d) * user.Nc[f]; 331 } 332 PetscCall(DMSetNumFields(dm, user.Nf)); 333 PetscCall(DMPlexCreateSection(dm, NULL, user.Nc, numDof, 0, NULL, NULL, NULL, NULL, &s)); 334 PetscCall(SetSymmetries(dm, s, &user)); 335 PetscCall(PetscFree(numDof)); 336 } 337 PetscCall(DMSetLocalSection(dm, s)); 338 /* Create spectral ordering and load in data */ 339 PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); 340 PetscCall(DMGetLocalVector(dm, &u)); 341 switch (dim) { 342 case 2: 343 PetscCall(LoadData2D(dm, 2, 2, size, u, &user)); 344 break; 345 case 3: 346 PetscCall(LoadData3D(dm, 2, 2, 2, size, u, &user)); 347 break; 348 } 349 /* Remove ordering and check some values */ 350 PetscCall(PetscSectionSetClosurePermutation(s, (PetscObject)dm, dim, NULL)); 351 switch (dim) { 352 case 2: 353 PetscCall(CheckPoint(dm, u, 0, &user)); 354 PetscCall(CheckPoint(dm, u, 13, &user)); 355 PetscCall(CheckPoint(dm, u, 15, &user)); 356 PetscCall(CheckPoint(dm, u, 19, &user)); 357 break; 358 case 3: 359 PetscCall(CheckPoint(dm, u, 0, &user)); 360 PetscCall(CheckPoint(dm, u, 13, &user)); 361 PetscCall(CheckPoint(dm, u, 15, &user)); 362 PetscCall(CheckPoint(dm, u, 19, &user)); 363 break; 364 } 365 /* Recreate spectral ordering and read out data */ 366 PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, s)); 367 switch (dim) { 368 case 2: 369 PetscCall(ReadData2D(dm, u, &user)); 370 break; 371 case 3: 372 PetscCall(ReadData3D(dm, u, &user)); 373 break; 374 } 375 PetscCall(DMRestoreLocalVector(dm, &u)); 376 PetscCall(PetscSectionDestroy(&s)); 377 PetscCall(DMDestroy(&dm)); 378 PetscCall(PetscFree(user.Nc)); 379 PetscCall(PetscFree(user.k)); 380 PetscCall(PetscFinalize()); 381 return 0; 382 } 383 384 /*TEST 385 386 # Spectral ordering 2D 0-5 387 testset: 388 args: -dm_plex_simplex 0 -dm_plex_box_faces 2,2 389 390 test: 391 suffix: 0 392 args: -num_fields 1 -num_components 1 -order 2 393 test: 394 suffix: 1 395 args: -num_fields 1 -num_components 1 -order 3 396 test: 397 suffix: 2 398 args: -num_fields 1 -num_components 1 -order 5 399 test: 400 suffix: 3 401 args: -num_fields 1 -num_components 2 -order 2 402 test: 403 suffix: 4 404 args: -num_fields 2 -num_components 1,1 -order 2,2 405 test: 406 suffix: 5 407 args: -num_fields 2 -num_components 1,2 -order 2,3 408 409 # Spectral ordering 3D 6-11 410 testset: 411 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 2,2,2 412 413 test: 414 suffix: 6 415 args: -num_fields 1 -num_components 1 -order 2 416 test: 417 suffix: 7 418 args: -num_fields 1 -num_components 1 -order 3 419 test: 420 suffix: 8 421 args: -num_fields 1 -num_components 1 -order 5 422 test: 423 suffix: 9 424 args: -num_fields 1 -num_components 2 -order 2 425 test: 426 suffix: 10 427 args: -num_fields 2 -num_components 1,1 -order 2,2 428 test: 429 suffix: 11 430 args: -num_fields 2 -num_components 1,2 -order 2,3 431 432 TEST*/ 433