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