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