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