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