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