1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petscdmplex.h> 3 4 static PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp) 5 { 6 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 7 PetscInt i; 8 PetscErrorCode ierr; 9 10 PetscFunctionBegin; 11 if (lag->symmetries) { 12 PetscInt **selfSyms = lag->symmetries[0]; 13 14 if (selfSyms) { 15 PetscInt i, **allocated = &selfSyms[-lag->selfSymOff]; 16 17 for (i = 0; i < lag->numSelfSym; i++) { 18 ierr = PetscFree(allocated[i]);CHKERRQ(ierr); 19 } 20 ierr = PetscFree(allocated);CHKERRQ(ierr); 21 } 22 ierr = PetscFree(lag->symmetries);CHKERRQ(ierr); 23 } 24 for (i = 0; i < lag->height; i++) { 25 ierr = PetscDualSpaceDestroy(&lag->subspaces[i]);CHKERRQ(ierr); 26 } 27 ierr = PetscFree(lag->subspaces);CHKERRQ(ierr); 28 ierr = PetscFree(lag);CHKERRQ(ierr); 29 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);CHKERRQ(ierr); 30 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);CHKERRQ(ierr); 31 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);CHKERRQ(ierr); 32 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 36 static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer) 37 { 38 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 39 PetscErrorCode ierr; 40 41 PetscFunctionBegin; 42 ierr = PetscViewerASCIIPrintf(viewer, "%s %sLagrange dual space\n", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "tensor " : "");CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 46 static PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer) 47 { 48 PetscBool iascii; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 53 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 54 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 55 if (iascii) {ierr = PetscDualSpaceLagrangeView_Ascii(sp, viewer);CHKERRQ(ierr);} 56 PetscFunctionReturn(0); 57 } 58 59 static PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp) 60 { 61 PetscBool continuous, tensor, flg; 62 PetscErrorCode ierr; 63 64 PetscFunctionBegin; 65 ierr = PetscDualSpaceLagrangeGetContinuity(sp, &continuous);CHKERRQ(ierr); 66 ierr = PetscDualSpaceLagrangeGetTensor(sp, &tensor);CHKERRQ(ierr); 67 ierr = PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");CHKERRQ(ierr); 68 ierr = PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);CHKERRQ(ierr); 69 if (flg) {ierr = PetscDualSpaceLagrangeSetContinuity(sp, continuous);CHKERRQ(ierr);} 70 ierr = PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetContinuity", tensor, &tensor, &flg);CHKERRQ(ierr); 71 if (flg) {ierr = PetscDualSpaceLagrangeSetTensor(sp, tensor);CHKERRQ(ierr);} 72 ierr = PetscOptionsTail();CHKERRQ(ierr); 73 PetscFunctionReturn(0); 74 } 75 76 static PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace spNew) 77 { 78 PetscBool cont, tensor; 79 PetscErrorCode ierr; 80 81 PetscFunctionBegin; 82 ierr = PetscDualSpaceLagrangeGetContinuity(sp, &cont);CHKERRQ(ierr); 83 ierr = PetscDualSpaceLagrangeSetContinuity(spNew, cont);CHKERRQ(ierr); 84 ierr = PetscDualSpaceLagrangeGetTensor(sp, &tensor);CHKERRQ(ierr); 85 ierr = PetscDualSpaceLagrangeSetTensor(spNew, tensor);CHKERRQ(ierr); 86 PetscFunctionReturn(0); 87 } 88 89 static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt *dim) 90 { 91 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 92 PetscReal D = 1.0; 93 PetscInt n, d; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 *dim = -1; 98 ierr = DMGetDimension(sp->dm, &n);CHKERRQ(ierr); 99 if (!lag->tensorSpace) { 100 for (d = 1; d <= n; ++d) { 101 D *= ((PetscReal) (order+d))/d; 102 } 103 *dim = (PetscInt) (D + 0.5); 104 } else { 105 *dim = 1; 106 for (d = 0; d < n; ++d) *dim *= (order+1); 107 } 108 *dim *= sp->Nc; 109 PetscFunctionReturn(0); 110 } 111 112 static PetscErrorCode PetscDualSpaceCreateHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp) 113 { 114 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 115 PetscBool continuous, tensor; 116 PetscInt order; 117 PetscErrorCode ierr; 118 119 PetscFunctionBegin; 120 ierr = PetscDualSpaceLagrangeGetContinuity(sp,&continuous);CHKERRQ(ierr); 121 ierr = PetscDualSpaceGetOrder(sp,&order);CHKERRQ(ierr); 122 if (height == 0) { 123 ierr = PetscObjectReference((PetscObject)sp);CHKERRQ(ierr); 124 *bdsp = sp; 125 } else if (continuous == PETSC_FALSE || !order) { 126 *bdsp = NULL; 127 } else { 128 DM dm, K; 129 PetscInt dim; 130 131 ierr = PetscDualSpaceGetDM(sp,&dm);CHKERRQ(ierr); 132 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 133 if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Asked for dual space at height %d for dimension %d reference element\n",height,dim); 134 ierr = PetscDualSpaceDuplicate(sp,bdsp);CHKERRQ(ierr); 135 ierr = PetscDualSpaceCreateReferenceCell(*bdsp, dim-height, lag->simplexCell, &K);CHKERRQ(ierr); 136 ierr = PetscDualSpaceSetDM(*bdsp, K);CHKERRQ(ierr); 137 ierr = DMDestroy(&K);CHKERRQ(ierr); 138 ierr = PetscDualSpaceLagrangeGetTensor(sp,&tensor);CHKERRQ(ierr); 139 ierr = PetscDualSpaceLagrangeSetTensor(*bdsp,tensor);CHKERRQ(ierr); 140 ierr = PetscDualSpaceSetUp(*bdsp);CHKERRQ(ierr); 141 } 142 PetscFunctionReturn(0); 143 } 144 145 static PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp) 146 { 147 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 148 DM dm = sp->dm; 149 PetscInt order = sp->order; 150 PetscInt Nc = sp->Nc; 151 MPI_Comm comm; 152 PetscBool continuous; 153 PetscSection csection; 154 Vec coordinates; 155 PetscReal *qpoints, *qweights; 156 PetscInt depth, dim, pdimMax, pStart, pEnd, p, *pStratStart, *pStratEnd, coneSize, d, f = 0, c; 157 PetscBool simplex, tensorSpace; 158 PetscErrorCode ierr; 159 160 PetscFunctionBegin; 161 ierr = PetscObjectGetComm((PetscObject) sp, &comm);CHKERRQ(ierr); 162 if (!order) lag->continuous = PETSC_FALSE; 163 continuous = lag->continuous; 164 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 165 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 166 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 167 ierr = PetscMalloc2(depth+1,&pStratStart,depth+1,&pStratEnd);CHKERRQ(ierr); 168 for (d = 0; d <= depth; ++d) {ierr = DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);CHKERRQ(ierr);} 169 ierr = DMPlexGetConeSize(dm, pStratStart[depth], &coneSize);CHKERRQ(ierr); 170 ierr = DMGetCoordinateSection(dm, &csection);CHKERRQ(ierr); 171 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 172 if (depth == 1) { 173 if (coneSize == dim+1) simplex = PETSC_TRUE; 174 else if (coneSize == 1 << dim) simplex = PETSC_FALSE; 175 else SETERRQ(comm, PETSC_ERR_SUP, "Only support simplices and tensor product cells"); 176 } else if (depth == dim) { 177 if (coneSize == dim+1) simplex = PETSC_TRUE; 178 else if (coneSize == 2 * dim) simplex = PETSC_FALSE; 179 else SETERRQ(comm, PETSC_ERR_SUP, "Only support simplices and tensor product cells"); 180 } else SETERRQ(comm, PETSC_ERR_SUP, "Only support cell-vertex meshes or fully interpolated meshes"); 181 lag->simplexCell = simplex; 182 if (dim > 1 && continuous && lag->simplexCell == lag->tensorSpace) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Mismatching simplex/tensor cells and spaces only allowed for discontinuous elements"); 183 tensorSpace = lag->tensorSpace; 184 lag->height = 0; 185 lag->subspaces = NULL; 186 if (continuous && order > 0 && dim > 0) { 187 PetscInt i; 188 189 lag->height = dim; 190 ierr = PetscMalloc1(dim,&lag->subspaces);CHKERRQ(ierr); 191 ierr = PetscDualSpaceCreateHeightSubspace_Lagrange(sp,1,&lag->subspaces[0]);CHKERRQ(ierr); 192 ierr = PetscDualSpaceSetUp(lag->subspaces[0]);CHKERRQ(ierr); 193 for (i = 1; i < dim; i++) { 194 ierr = PetscDualSpaceGetHeightSubspace(lag->subspaces[i-1],1,&lag->subspaces[i]);CHKERRQ(ierr); 195 ierr = PetscObjectReference((PetscObject)(lag->subspaces[i]));CHKERRQ(ierr); 196 } 197 } 198 ierr = PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);CHKERRQ(ierr); 199 pdimMax *= (pStratEnd[depth] - pStratStart[depth]); 200 ierr = PetscMalloc1(pdimMax, &sp->functional);CHKERRQ(ierr); 201 if (!dim) { 202 for (c = 0; c < Nc; ++c) { 203 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);CHKERRQ(ierr); 204 ierr = PetscCalloc1(Nc, &qweights);CHKERRQ(ierr); 205 ierr = PetscQuadratureSetOrder(sp->functional[f], 0);CHKERRQ(ierr); 206 ierr = PetscQuadratureSetData(sp->functional[f], 0, Nc, 1, NULL, qweights);CHKERRQ(ierr); 207 qweights[c] = 1.0; 208 ++f; 209 } 210 } else { 211 PetscSection section; 212 PetscReal *v0, *hv0, *J, *invJ, detJ, hdetJ; 213 PetscInt *tup; 214 215 ierr = PetscSectionCreate(PETSC_COMM_SELF,§ion);CHKERRQ(ierr); 216 ierr = PetscSectionSetChart(section,pStart,pEnd);CHKERRQ(ierr); 217 ierr = PetscCalloc5(dim+1,&tup,dim,&v0,dim,&hv0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 218 for (p = pStart; p < pEnd; p++) { 219 PetscInt pointDim, d, nFunc = 0; 220 PetscDualSpace hsp; 221 222 ierr = DMPlexComputeCellGeometryFEM(dm, p, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 223 for (d = 0; d < depth; d++) {if (p >= pStratStart[d] && p < pStratEnd[d]) break;} 224 pointDim = (depth == 1 && d == 1) ? dim : d; 225 hsp = ((pointDim < dim) && lag->subspaces) ? lag->subspaces[dim - pointDim - 1] : NULL; 226 if (hsp) { 227 DM hdm; 228 229 ierr = PetscDualSpaceGetDM(hsp,&hdm);CHKERRQ(ierr); 230 ierr = DMPlexComputeCellGeometryFEM(hdm, 0, NULL, hv0, NULL, NULL, &hdetJ);CHKERRQ(ierr); 231 nFunc = sp->numDof[pointDim]; 232 } 233 if (pointDim == dim) { 234 /* Cells, create for self */ 235 PetscInt orderEff = continuous ? (!tensorSpace ? order-1-dim : order-2) : order; 236 PetscReal denom = continuous ? order : (!tensorSpace ? order+1+dim : order+2); 237 PetscReal numer = (!simplex || !tensorSpace) ? 2. : (2./dim); 238 PetscReal dx = numer/denom; 239 PetscInt cdim, d, d2; 240 241 if (orderEff < 0) continue; 242 ierr = PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, orderEff, &cdim);CHKERRQ(ierr); 243 ierr = PetscArrayzero(tup,dim+1);CHKERRQ(ierr); 244 if (!tensorSpace) { 245 while (!tup[dim]) { 246 for (c = 0; c < Nc; ++c) { 247 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);CHKERRQ(ierr); 248 ierr = PetscMalloc1(dim, &qpoints);CHKERRQ(ierr); 249 ierr = PetscCalloc1(Nc, &qweights);CHKERRQ(ierr); 250 ierr = PetscQuadratureSetOrder(sp->functional[f], 0);CHKERRQ(ierr); 251 ierr = PetscQuadratureSetData(sp->functional[f], dim, Nc, 1, qpoints, qweights);CHKERRQ(ierr); 252 for (d = 0; d < dim; ++d) { 253 qpoints[d] = v0[d]; 254 for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx); 255 } 256 qweights[c] = 1.0; 257 ++f; 258 } 259 ierr = PetscDualSpaceLatticePointLexicographic_Internal(dim, orderEff, tup);CHKERRQ(ierr); 260 } 261 } else { 262 while (!tup[dim]) { 263 for (c = 0; c < Nc; ++c) { 264 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);CHKERRQ(ierr); 265 ierr = PetscMalloc1(dim, &qpoints);CHKERRQ(ierr); 266 ierr = PetscCalloc1(Nc, &qweights);CHKERRQ(ierr); 267 ierr = PetscQuadratureSetOrder(sp->functional[f], 0);CHKERRQ(ierr); 268 ierr = PetscQuadratureSetData(sp->functional[f], dim, Nc, 1, qpoints, qweights);CHKERRQ(ierr); 269 for (d = 0; d < dim; ++d) { 270 qpoints[d] = v0[d]; 271 for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx); 272 } 273 qweights[c] = 1.0; 274 ++f; 275 } 276 ierr = PetscDualSpaceTensorPointLexicographic_Internal(dim, orderEff, tup);CHKERRQ(ierr); 277 } 278 } 279 } else { /* transform functionals from subspaces */ 280 PetscInt q; 281 282 for (q = 0; q < nFunc; q++, f++) { 283 PetscQuadrature fn; 284 PetscInt fdim, Nc, c, nPoints, i; 285 const PetscReal *points; 286 const PetscReal *weights; 287 PetscReal *qpoints; 288 PetscReal *qweights; 289 290 ierr = PetscDualSpaceGetFunctional(hsp, q, &fn);CHKERRQ(ierr); 291 ierr = PetscQuadratureGetData(fn,&fdim,&Nc,&nPoints,&points,&weights);CHKERRQ(ierr); 292 if (fdim != pointDim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected height dual space dim %D, got %D",pointDim,fdim); 293 ierr = PetscMalloc1(nPoints * dim, &qpoints);CHKERRQ(ierr); 294 ierr = PetscCalloc1(nPoints * Nc, &qweights);CHKERRQ(ierr); 295 for (i = 0; i < nPoints; i++) { 296 PetscInt j, k; 297 PetscReal *qp = &qpoints[i * dim]; 298 299 for (c = 0; c < Nc; ++c) qweights[i*Nc+c] = weights[i*Nc+c]; 300 for (j = 0; j < dim; ++j) qp[j] = v0[j]; 301 for (j = 0; j < dim; ++j) { 302 for (k = 0; k < pointDim; k++) qp[j] += J[dim * j + k] * (points[pointDim * i + k] - hv0[k]); 303 } 304 } 305 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);CHKERRQ(ierr); 306 ierr = PetscQuadratureSetOrder(sp->functional[f],0);CHKERRQ(ierr); 307 ierr = PetscQuadratureSetData(sp->functional[f],dim,Nc,nPoints,qpoints,qweights);CHKERRQ(ierr); 308 } 309 } 310 } 311 ierr = PetscFree5(tup,v0,hv0,J,invJ);CHKERRQ(ierr); 312 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 313 { /* reorder to closure order */ 314 PetscInt *key, count; 315 PetscQuadrature *reorder = NULL; 316 317 ierr = PetscCalloc1(f,&key);CHKERRQ(ierr); 318 ierr = PetscMalloc1(f*sp->Nc,&reorder);CHKERRQ(ierr); 319 320 for (p = pStratStart[depth], count = 0; p < pStratEnd[depth]; p++) { 321 PetscInt *closure = NULL, closureSize, c; 322 323 ierr = DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 324 for (c = 0; c < closureSize; c++) { 325 PetscInt point = closure[2 * c], dof, off, i; 326 327 ierr = PetscSectionGetDof(section,point,&dof);CHKERRQ(ierr); 328 ierr = PetscSectionGetOffset(section,point,&off);CHKERRQ(ierr); 329 for (i = 0; i < dof; i++) { 330 PetscInt fi = i + off; 331 if (!key[fi]) { 332 key[fi] = 1; 333 reorder[count++] = sp->functional[fi]; 334 } 335 } 336 } 337 ierr = DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 338 } 339 ierr = PetscFree(sp->functional);CHKERRQ(ierr); 340 sp->functional = reorder; 341 ierr = PetscFree(key);CHKERRQ(ierr); 342 } 343 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 344 } 345 if (pStratEnd[depth] == 1 && f != pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %D not equal to dimension %D", f, pdimMax); 346 if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %D is greater than max size %D", f, pdimMax); 347 ierr = PetscFree2(pStratStart, pStratEnd);CHKERRQ(ierr); 348 PetscFunctionReturn(0); 349 } 350 351 static PetscErrorCode PetscDualSpaceGetHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp) 352 { 353 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 354 PetscErrorCode ierr; 355 356 PetscFunctionBegin; 357 if (height == 0) { 358 *bdsp = sp; 359 } else { 360 DM dm; 361 PetscInt dim; 362 363 ierr = PetscDualSpaceGetDM(sp,&dm);CHKERRQ(ierr); 364 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 365 if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Asked for dual space at height %D for dimension %D reference element\n",height,dim); 366 if (height <= lag->height) {*bdsp = lag->subspaces[height-1];} 367 else {*bdsp = NULL;} 368 } 369 PetscFunctionReturn(0); 370 } 371 372 #define BaryIndex(perEdge,a,b,c) (((b)*(2*perEdge+1-(b)))/2)+(c) 373 374 #define CartIndex(perEdge,a,b) (perEdge*(a)+b) 375 376 static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips) 377 { 378 379 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 380 PetscInt dim, order, p, Nc; 381 PetscErrorCode ierr; 382 383 PetscFunctionBegin; 384 ierr = PetscDualSpaceGetOrder(sp,&order);CHKERRQ(ierr); 385 ierr = PetscDualSpaceGetNumComponents(sp,&Nc);CHKERRQ(ierr); 386 ierr = DMGetDimension(sp->dm,&dim);CHKERRQ(ierr); 387 if (!dim || !lag->continuous || order < 3) PetscFunctionReturn(0); 388 if (dim > 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Lagrange symmetries not implemented for dim = %D > 3",dim); 389 if (!lag->symmetries) { /* store symmetries */ 390 PetscDualSpace hsp; 391 DM K; 392 PetscInt numPoints = 1, d; 393 PetscInt numFaces; 394 PetscInt ***symmetries; 395 const PetscInt ***hsymmetries; 396 397 if (lag->simplexCell) { 398 numFaces = 1 + dim; 399 for (d = 0; d < dim; d++) numPoints = numPoints * 2 + 1; 400 } else { 401 numPoints = PetscPowInt(3,dim); 402 numFaces = 2 * dim; 403 } 404 ierr = PetscCalloc1(numPoints,&symmetries);CHKERRQ(ierr); 405 if (0 < dim && dim < 3) { /* compute self symmetries */ 406 PetscInt **cellSymmetries; 407 408 lag->numSelfSym = 2 * numFaces; 409 lag->selfSymOff = numFaces; 410 ierr = PetscCalloc1(2*numFaces,&cellSymmetries);CHKERRQ(ierr); 411 /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */ 412 symmetries[0] = &cellSymmetries[numFaces]; 413 if (dim == 1) { 414 PetscInt dofPerEdge = order - 1; 415 416 if (dofPerEdge > 1) { 417 PetscInt i, j, *reverse; 418 419 ierr = PetscMalloc1(dofPerEdge*Nc,&reverse);CHKERRQ(ierr); 420 for (i = 0; i < dofPerEdge; i++) { 421 for (j = 0; j < Nc; j++) { 422 reverse[i*Nc + j] = Nc * (dofPerEdge - 1 - i) + j; 423 } 424 } 425 symmetries[0][-2] = reverse; 426 427 /* yes, this is redundant, but it makes it easier to cleanup if I don't have to worry about what not to free */ 428 ierr = PetscMalloc1(dofPerEdge*Nc,&reverse);CHKERRQ(ierr); 429 for (i = 0; i < dofPerEdge; i++) { 430 for (j = 0; j < Nc; j++) { 431 reverse[i*Nc + j] = Nc * (dofPerEdge - 1 - i) + j; 432 } 433 } 434 symmetries[0][1] = reverse; 435 } 436 } else { 437 PetscInt dofPerEdge = lag->simplexCell ? (order - 2) : (order - 1), s; 438 PetscInt dofPerFace; 439 440 if (dofPerEdge > 1) { 441 for (s = -numFaces; s < numFaces; s++) { 442 PetscInt *sym, i, j, k, l; 443 444 if (!s) continue; 445 if (lag->simplexCell) { 446 dofPerFace = (dofPerEdge * (dofPerEdge + 1))/2; 447 ierr = PetscMalloc1(Nc*dofPerFace,&sym);CHKERRQ(ierr); 448 for (j = 0, l = 0; j < dofPerEdge; j++) { 449 for (k = 0; k < dofPerEdge - j; k++, l++) { 450 i = dofPerEdge - 1 - j - k; 451 switch (s) { 452 case -3: 453 sym[Nc*l] = BaryIndex(dofPerEdge,i,k,j); 454 break; 455 case -2: 456 sym[Nc*l] = BaryIndex(dofPerEdge,j,i,k); 457 break; 458 case -1: 459 sym[Nc*l] = BaryIndex(dofPerEdge,k,j,i); 460 break; 461 case 1: 462 sym[Nc*l] = BaryIndex(dofPerEdge,k,i,j); 463 break; 464 case 2: 465 sym[Nc*l] = BaryIndex(dofPerEdge,j,k,i); 466 break; 467 } 468 } 469 } 470 } else { 471 dofPerFace = dofPerEdge * dofPerEdge; 472 ierr = PetscMalloc1(Nc*dofPerFace,&sym);CHKERRQ(ierr); 473 for (j = 0, l = 0; j < dofPerEdge; j++) { 474 for (k = 0; k < dofPerEdge; k++, l++) { 475 switch (s) { 476 case -4: 477 sym[Nc*l] = CartIndex(dofPerEdge,k,j); 478 break; 479 case -3: 480 sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),k); 481 break; 482 case -2: 483 sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),(dofPerEdge - 1 - j)); 484 break; 485 case -1: 486 sym[Nc*l] = CartIndex(dofPerEdge,j,(dofPerEdge - 1 - k)); 487 break; 488 case 1: 489 sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),j); 490 break; 491 case 2: 492 sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),(dofPerEdge - 1 - k)); 493 break; 494 case 3: 495 sym[Nc*l] = CartIndex(dofPerEdge,k,(dofPerEdge - 1 - j)); 496 break; 497 } 498 } 499 } 500 } 501 for (i = 0; i < dofPerFace; i++) { 502 sym[Nc*i] *= Nc; 503 for (j = 1; j < Nc; j++) { 504 sym[Nc*i+j] = sym[Nc*i] + j; 505 } 506 } 507 symmetries[0][s] = sym; 508 } 509 } 510 } 511 } 512 ierr = PetscDualSpaceGetHeightSubspace(sp,1,&hsp);CHKERRQ(ierr); 513 ierr = PetscDualSpaceGetSymmetries(hsp,&hsymmetries,NULL);CHKERRQ(ierr); 514 if (hsymmetries) { 515 PetscBool *seen; 516 const PetscInt *cone; 517 PetscInt KclosureSize, *Kclosure = NULL; 518 519 ierr = PetscDualSpaceGetDM(sp,&K);CHKERRQ(ierr); 520 ierr = PetscCalloc1(numPoints,&seen);CHKERRQ(ierr); 521 ierr = DMPlexGetCone(K,0,&cone);CHKERRQ(ierr); 522 ierr = DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&KclosureSize,&Kclosure);CHKERRQ(ierr); 523 for (p = 0; p < numFaces; p++) { 524 PetscInt closureSize, *closure = NULL, q; 525 526 ierr = DMPlexGetTransitiveClosure(K,cone[p],PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 527 for (q = 0; q < closureSize; q++) { 528 PetscInt point = closure[2*q], r; 529 530 if(!seen[point]) { 531 for (r = 0; r < KclosureSize; r++) { 532 if (Kclosure[2 * r] == point) break; 533 } 534 seen[point] = PETSC_TRUE; 535 symmetries[r] = (PetscInt **) hsymmetries[q]; 536 } 537 } 538 ierr = DMPlexRestoreTransitiveClosure(K,cone[p],PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 539 } 540 ierr = DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&KclosureSize,&Kclosure);CHKERRQ(ierr); 541 ierr = PetscFree(seen);CHKERRQ(ierr); 542 } 543 lag->symmetries = symmetries; 544 } 545 if (perms) *perms = (const PetscInt ***) lag->symmetries; 546 PetscFunctionReturn(0); 547 } 548 549 static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous) 550 { 551 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 552 553 PetscFunctionBegin; 554 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 555 PetscValidPointer(continuous, 2); 556 *continuous = lag->continuous; 557 PetscFunctionReturn(0); 558 } 559 560 static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous) 561 { 562 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 563 564 PetscFunctionBegin; 565 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 566 lag->continuous = continuous; 567 PetscFunctionReturn(0); 568 } 569 570 /*@ 571 PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity 572 573 Not Collective 574 575 Input Parameter: 576 . sp - the PetscDualSpace 577 578 Output Parameter: 579 . continuous - flag for element continuity 580 581 Level: intermediate 582 583 .seealso: PetscDualSpaceLagrangeSetContinuity() 584 @*/ 585 PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous) 586 { 587 PetscErrorCode ierr; 588 589 PetscFunctionBegin; 590 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 591 PetscValidPointer(continuous, 2); 592 ierr = PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));CHKERRQ(ierr); 593 PetscFunctionReturn(0); 594 } 595 596 /*@ 597 PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous 598 599 Logically Collective on sp 600 601 Input Parameters: 602 + sp - the PetscDualSpace 603 - continuous - flag for element continuity 604 605 Options Database: 606 . -petscdualspace_lagrange_continuity <bool> 607 608 Level: intermediate 609 610 .seealso: PetscDualSpaceLagrangeGetContinuity() 611 @*/ 612 PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous) 613 { 614 PetscErrorCode ierr; 615 616 PetscFunctionBegin; 617 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 618 PetscValidLogicalCollectiveBool(sp, continuous, 2); 619 ierr = PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));CHKERRQ(ierr); 620 PetscFunctionReturn(0); 621 } 622 623 static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor) 624 { 625 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 626 627 PetscFunctionBegin; 628 *tensor = lag->tensorSpace; 629 PetscFunctionReturn(0); 630 } 631 632 static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor) 633 { 634 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 635 636 PetscFunctionBegin; 637 lag->tensorSpace = tensor; 638 PetscFunctionReturn(0); 639 } 640 641 /*@ 642 PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space 643 644 Not collective 645 646 Input Parameter: 647 . sp - The PetscDualSpace 648 649 Output Parameter: 650 . tensor - Whether the dual space has tensor layout (vs. simplicial) 651 652 Level: intermediate 653 654 .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate() 655 @*/ 656 PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor) 657 { 658 PetscErrorCode ierr; 659 660 PetscFunctionBegin; 661 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 662 PetscValidPointer(tensor, 2); 663 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));CHKERRQ(ierr); 664 PetscFunctionReturn(0); 665 } 666 667 /*@ 668 PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space 669 670 Not collective 671 672 Input Parameters: 673 + sp - The PetscDualSpace 674 - tensor - Whether the dual space has tensor layout (vs. simplicial) 675 676 Level: intermediate 677 678 .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate() 679 @*/ 680 PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor) 681 { 682 PetscErrorCode ierr; 683 684 PetscFunctionBegin; 685 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 686 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));CHKERRQ(ierr); 687 PetscFunctionReturn(0); 688 } 689 690 static PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp) 691 { 692 PetscFunctionBegin; 693 sp->ops->destroy = PetscDualSpaceDestroy_Lagrange; 694 sp->ops->view = PetscDualSpaceView_Lagrange; 695 sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange; 696 sp->ops->duplicate = PetscDualSpaceDuplicate_Lagrange; 697 sp->ops->setup = PetscDualSpaceSetUp_Lagrange; 698 sp->ops->createheightsubspace = PetscDualSpaceGetHeightSubspace_Lagrange; 699 sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_Lagrange; 700 sp->ops->apply = PetscDualSpaceApplyDefault; 701 sp->ops->applyall = PetscDualSpaceApplyAllDefault; 702 sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault; 703 sp->ops->applyint = PetscDualSpaceApplyInteriorDefault; 704 sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault; 705 PetscFunctionReturn(0); 706 } 707 708 /*MC 709 PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals 710 711 Level: intermediate 712 713 .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType() 714 M*/ 715 716 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp) 717 { 718 PetscDualSpace_Lag *lag; 719 PetscErrorCode ierr; 720 721 PetscFunctionBegin; 722 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 723 ierr = PetscNewLog(sp,&lag);CHKERRQ(ierr); 724 sp->data = lag; 725 726 lag->simplexCell = PETSC_TRUE; 727 lag->tensorSpace = PETSC_FALSE; 728 lag->continuous = PETSC_TRUE; 729 730 ierr = PetscDualSpaceInitialize_Lagrange(sp);CHKERRQ(ierr); 731 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);CHKERRQ(ierr); 732 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);CHKERRQ(ierr); 733 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);CHKERRQ(ierr); 734 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);CHKERRQ(ierr); 735 PetscFunctionReturn(0); 736 } 737 738