1 #include <petsc/private/dmfieldimpl.h> /*I "petscdmfield.h" I*/ 2 #include <petsc/private/petscfeimpl.h> /*I "petscdmfield.h" I*/ 3 #include <petscfe.h> 4 #include <petscdmplex.h> 5 #include <petscds.h> 6 7 typedef struct _n_DMField_DS 8 { 9 PetscInt fieldNum; 10 Vec vec; 11 PetscInt height; 12 PetscObject *disc; 13 PetscBool multifieldVec; 14 } 15 DMField_DS; 16 17 static PetscErrorCode DMFieldDestroy_DS(DMField field) 18 { 19 DMField_DS *dsfield; 20 PetscInt i; 21 PetscErrorCode ierr; 22 23 PetscFunctionBegin; 24 dsfield = (DMField_DS *) field->data; 25 ierr = VecDestroy(&dsfield->vec);CHKERRQ(ierr); 26 for (i = 0; i < dsfield->height; i++) { 27 ierr = PetscObjectDereference(dsfield->disc[i]);CHKERRQ(ierr); 28 } 29 ierr = PetscFree(dsfield->disc);CHKERRQ(ierr); 30 ierr = PetscFree(dsfield);CHKERRQ(ierr); 31 PetscFunctionReturn(0); 32 } 33 34 static PetscErrorCode DMFieldView_DS(DMField field,PetscViewer viewer) 35 { 36 DMField_DS *dsfield = (DMField_DS *) field->data; 37 PetscBool iascii; 38 PetscObject disc; 39 PetscErrorCode ierr; 40 41 PetscFunctionBegin; 42 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 43 disc = dsfield->disc[0]; 44 if (iascii) { 45 PetscViewerASCIIPrintf(viewer, "PetscDS field %D\n",dsfield->fieldNum);CHKERRQ(ierr); 46 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 47 ierr = PetscObjectView(disc,viewer);CHKERRQ(ierr); 48 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 49 } 50 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 51 if (dsfield->multifieldVec) { 52 SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"View of subfield not implemented yet"); 53 } else { 54 ierr = VecView(dsfield->vec,viewer);CHKERRQ(ierr); 55 } 56 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 57 PetscFunctionReturn(0); 58 } 59 60 static PetscErrorCode DMFieldEvaluate_DS(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H) 61 { 62 PetscFunctionBegin; 63 SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented yet"); 64 PetscFunctionReturn(0); 65 } 66 67 static PetscErrorCode DMFieldDSGetHeightDisc(DMField field, PetscInt height, PetscObject *disc) 68 { 69 DMField_DS *dsfield = (DMField_DS *) field->data; 70 PetscErrorCode ierr; 71 72 PetscFunctionBegin; 73 if (!dsfield->disc[height]) { 74 PetscClassId id; 75 76 ierr = PetscObjectGetClassId(dsfield->disc[0],&id);CHKERRQ(ierr); 77 if (id == PETSCFE_CLASSID) { 78 PetscFE fe = (PetscFE) dsfield->disc[0]; 79 80 ierr = PetscFECreateHeightTrace(fe,height,(PetscFE *)&dsfield->disc[height]);CHKERRQ(ierr); 81 } 82 } 83 *disc = dsfield->disc[height]; 84 PetscFunctionReturn(0); 85 } 86 87 #define DMFieldDSdot(y,A,b,m,n,c,cast) \ 88 do { \ 89 PetscInt _i, _j, _k; \ 90 for (_i = 0; _i < (m); _i++) { \ 91 for (_k = 0; _k < (c); _k++) { \ 92 (y)[_i * (c) + _k] = 0.; \ 93 } \ 94 for (_j = 0; _j < (n); _j++) { \ 95 for (_k = 0; _k < (c); _k++) { \ 96 (y)[_i * (c) + _k] += (A)[(_i * (n) + _j) * (c) + _k] * cast((b)[_j]); \ 97 } \ 98 } \ 99 } \ 100 } while (0) 101 102 static PetscErrorCode DMFieldEvaluateFE_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H) 103 { 104 DMField_DS *dsfield = (DMField_DS *) field->data; 105 DM dm; 106 PetscObject disc; 107 PetscClassId classid; 108 PetscInt nq, nc, dim, meshDim, numCells; 109 PetscSection section; 110 const PetscReal *qpoints; 111 PetscBool isStride; 112 const PetscInt *points = NULL; 113 PetscInt sfirst = -1, stride = -1; 114 PetscErrorCode ierr; 115 116 PetscFunctionBeginHot; 117 dm = field->dm; 118 nc = field->numComponents; 119 ierr = PetscQuadratureGetData(quad,&dim,NULL,&nq,&qpoints,NULL);CHKERRQ(ierr); 120 ierr = DMFieldDSGetHeightDisc(field,dsfield->height - 1 - dim,&disc);CHKERRQ(ierr); 121 ierr = DMGetDimension(dm,&meshDim);CHKERRQ(ierr); 122 ierr = DMGetDefaultSection(dm,§ion);CHKERRQ(ierr); 123 ierr = PetscSectionGetField(section,dsfield->fieldNum,§ion);CHKERRQ(ierr); 124 ierr = PetscObjectGetClassId(disc,&classid);CHKERRQ(ierr); 125 /* TODO: batch */ 126 ierr = PetscObjectTypeCompare((PetscObject)pointIS,ISSTRIDE,&isStride);CHKERRQ(ierr); 127 ierr = ISGetLocalSize(pointIS,&numCells);CHKERRQ(ierr); 128 if (isStride) { 129 ierr = ISStrideGetInfo(pointIS,&sfirst,&stride);CHKERRQ(ierr); 130 } else { 131 ierr = ISGetIndices(pointIS,&points);CHKERRQ(ierr); 132 } 133 if (classid == PETSCFE_CLASSID) { 134 PetscFE fe = (PetscFE) disc; 135 PetscInt feDim, i; 136 PetscReal *fB = NULL, *fD = NULL, *fH = NULL; 137 138 if (dim == meshDim - 1) { 139 /* TODO */ 140 } 141 ierr = PetscFEGetDimension(fe,&feDim);CHKERRQ(ierr); 142 ierr = PetscFEGetTabulation(fe,nq,qpoints,B ? &fB : NULL,D ? &fD : NULL,H ? &fH : NULL);CHKERRQ(ierr); 143 for (i = 0; i < numCells; i++) { 144 PetscInt c = isStride ? (sfirst + i * stride) : points[i]; 145 PetscInt closureSize; 146 PetscScalar *elem = NULL; 147 148 ierr = DMPlexVecGetClosure(dm,section,dsfield->vec,c,&closureSize,&elem);CHKERRQ(ierr); 149 if (B) { 150 if (type == PETSC_SCALAR) { 151 PetscScalar *cB = &((PetscScalar *) B)[nc * nq * i]; 152 153 DMFieldDSdot(cB,fB,elem,nq,feDim,nc,(PetscScalar)); 154 } else { 155 PetscReal *cB = &((PetscReal *) B)[nc * nq * i]; 156 157 DMFieldDSdot(cB,fB,elem,nq,feDim,nc,PetscRealPart); 158 } 159 } 160 if (D) { 161 if (type == PETSC_SCALAR) { 162 PetscScalar *cD = &((PetscScalar *) D)[nc * nq * dim * i]; 163 164 DMFieldDSdot(cD,fD,elem,nq,feDim,(nc * dim),(PetscScalar)); 165 } else { 166 PetscReal *cD = &((PetscReal *) D)[nc * nq * dim * i]; 167 168 DMFieldDSdot(cD,fD,elem,nq,feDim,(nc * dim),PetscRealPart); 169 } 170 } 171 if (H) { 172 if (type == PETSC_SCALAR) { 173 PetscScalar *cH = &((PetscScalar *) H)[nc * nq * dim * dim * i]; 174 175 DMFieldDSdot(cH,fH,elem,nq,feDim,(nc * dim * dim),(PetscScalar)); 176 } else { 177 PetscReal *cH = &((PetscReal *) H)[nc * nq * dim * dim * i]; 178 179 DMFieldDSdot(cH,fH,elem,nq,feDim,(nc * dim * dim),PetscRealPart); 180 } 181 } 182 ierr = DMPlexVecRestoreClosure(dm,section,dsfield->vec,c,&closureSize,&elem);CHKERRQ(ierr); 183 } 184 ierr = PetscFERestoreTabulation(fe,nq,qpoints,B ? &fB : NULL,D ? &fD : NULL,H ? &fH : NULL);CHKERRQ(ierr); 185 } else {SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented");} 186 if (!isStride) { 187 ierr = ISRestoreIndices(pointIS,&points);CHKERRQ(ierr); 188 } 189 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 190 PetscFunctionReturn(0); 191 } 192 193 static PetscErrorCode DMFieldGetFEInvariance_DS(DMField field, IS pointIS, PetscBool *isConstant, PetscBool *isAffine, PetscBool *isQuadratic) 194 { 195 DMField_DS *dsfield; 196 PetscObject disc; 197 PetscInt h, imin; 198 PetscClassId id; 199 PetscErrorCode ierr; 200 201 PetscFunctionBegin; 202 dsfield = (DMField_DS *) field->data; 203 ierr = ISGetMinMax(pointIS,&imin,NULL);CHKERRQ(ierr); 204 for (h = 0; h < dsfield->height; h++) { 205 PetscInt hEnd; 206 207 ierr = DMPlexGetHeightStratum(field->dm,h,NULL,&hEnd);CHKERRQ(ierr); 208 if (imin < hEnd) break; 209 } 210 ierr = DMFieldDSGetHeightDisc(field,h,&disc);CHKERRQ(ierr); 211 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 212 if (id == PETSCFE_CLASSID) { 213 PetscFE fe = (PetscFE) disc; 214 PetscInt order, maxOrder; 215 PetscBool tensor = PETSC_FALSE; 216 PetscSpace sp; 217 218 ierr = PetscFEGetBasisSpace(fe, &sp);CHKERRQ(ierr); 219 ierr = PetscSpaceGetOrder(sp,&order);CHKERRQ(ierr); 220 ierr = PetscSpacePolynomialGetTensor(sp,&tensor);CHKERRQ(ierr); 221 if (tensor) { 222 PetscInt dim; 223 224 ierr = DMGetDimension(field->dm,&dim);CHKERRQ(ierr); 225 maxOrder = order * dim; 226 } else { 227 maxOrder = order; 228 } 229 if (isConstant) *isConstant = (maxOrder < 1) ? PETSC_TRUE : PETSC_FALSE; 230 if (isAffine) *isAffine = (maxOrder < 2) ? PETSC_TRUE : PETSC_FALSE; 231 if (isQuadratic) *isQuadratic = (maxOrder < 3) ? PETSC_TRUE : PETSC_FALSE; 232 } 233 PetscFunctionReturn(0); 234 } 235 236 static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad) 237 { 238 PetscInt h, dim, imax, imin; 239 DM dm; 240 DMField_DS *dsfield; 241 PetscObject disc; 242 PetscFE fe; 243 PetscClassId id; 244 PetscErrorCode ierr; 245 246 247 PetscFunctionBegin; 248 dm = field->dm; 249 dsfield = (DMField_DS *) field->data; 250 ierr = ISGetMinMax(pointIS,&imax,&imin);CHKERRQ(ierr); 251 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 252 for (h = 0; h <= dim; h++) { 253 PetscInt hStart, hEnd; 254 255 ierr = DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);CHKERRQ(ierr); 256 if (imin >= hStart && imax < hEnd) break; 257 } 258 *quad = NULL; 259 if (h < dsfield->height) { 260 ierr = DMFieldDSGetHeightDisc(field,h,&disc);CHKERRQ(ierr); 261 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 262 if (id != PETSCFE_CLASSID) PetscFunctionReturn(0); 263 fe = (PetscFE) disc; 264 ierr = PetscFEGetQuadrature(fe,quad);CHKERRQ(ierr); 265 ierr = PetscObjectReference((PetscObject)*quad);CHKERRQ(ierr); 266 } 267 PetscFunctionReturn(0); 268 } 269 270 static PetscErrorCode DMFieldComputeFaceData_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscFEGeom *geom) 271 { 272 const PetscInt *points; 273 PetscInt p, dim, dE, numFaces, Nq; 274 PetscBool affineCells; 275 DMLabel depthLabel; 276 IS cellIS; 277 DM dm = field->dm; 278 PetscErrorCode ierr; 279 280 PetscFunctionBegin; 281 dim = geom->dim; 282 dE = geom->dimEmbed; 283 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 284 ierr = DMLabelGetStratumIS(depthLabel, dim + 1, &cellIS);CHKERRQ(ierr); 285 ierr = DMFieldGetFEInvariance(field,cellIS,NULL,&affineCells,NULL);CHKERRQ(ierr); 286 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 287 numFaces = geom->numCells; 288 Nq = geom->numPoints; 289 if (affineCells) { 290 PetscInt numCells, offset, *cells; 291 PetscFEGeom *cellGeom; 292 IS suppIS; 293 PetscQuadrature cellQuad = NULL; 294 295 ierr = DMFieldCreateDefaultQuadrature(field,cellIS,&cellQuad);CHKERRQ(ierr); 296 for (p = 0, numCells = 0; p < numFaces; p++) { 297 PetscInt point = points[p]; 298 PetscInt numSupp, numChildren; 299 300 ierr = DMPlexGetTreeChildren(dm, point, &numChildren, NULL); CHKERRQ(ierr); 301 if (numChildren) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children"); 302 ierr = DMPlexGetSupportSize(dm, point,&numSupp);CHKERRQ(ierr); 303 if (numSupp > 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D has %D support, expected at most 2\n", point, numSupp); 304 numCells += numSupp; 305 } 306 ierr = PetscMalloc1(numCells, &cells);CHKERRQ(ierr); 307 for (p = 0, offset = 0; p < numFaces; p++) { 308 PetscInt point = points[p]; 309 PetscInt numSupp, s; 310 const PetscInt *supp; 311 312 ierr = DMPlexGetSupportSize(dm, point,&numSupp);CHKERRQ(ierr); 313 ierr = DMPlexGetSupport(dm, point, &supp);CHKERRQ(ierr); 314 for (s = 0; s < numSupp; s++, offset++) { 315 cells[offset] = supp[s]; 316 } 317 } 318 ierr = ISCreateGeneral(PETSC_COMM_SELF,numCells,cells,PETSC_USE_POINTER, &suppIS);CHKERRQ(ierr); 319 ierr = DMFieldCreateFEGeom(field,suppIS,cellQuad,PETSC_FALSE,&cellGeom);CHKERRQ(ierr); 320 for (p = 0, offset = 0; p < numFaces; p++) { 321 PetscInt point = points[p]; 322 PetscInt numSupp, s, q; 323 const PetscInt *supp; 324 325 ierr = DMPlexGetSupportSize(dm, point,&numSupp);CHKERRQ(ierr); 326 ierr = DMPlexGetSupport(dm, point, &supp);CHKERRQ(ierr); 327 for (s = 0; s < numSupp; s++, offset++) { 328 for (q = 0; q < Nq * dE * dE; q++) { 329 geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q]; 330 } 331 } 332 } 333 ierr = PetscFEGeomDestroy(&cellGeom);CHKERRQ(ierr); 334 ierr = ISDestroy(&suppIS);CHKERRQ(ierr); 335 ierr = PetscFree(cells);CHKERRQ(ierr); 336 ierr = PetscQuadratureDestroy(&cellQuad);CHKERRQ(ierr); 337 } else { 338 PetscObject faceDisc, cellDisc; 339 PetscClassId faceId, cellId; 340 PetscDualSpace dsp; 341 DM K; 342 PetscInt (*co)[2][3]; 343 PetscInt coneSize; 344 PetscInt **counts; 345 PetscInt f, i, o, q, s; 346 const PetscInt *coneK; 347 PetscInt minOrient, maxOrient, numOrient; 348 PetscInt *orients; 349 PetscReal **orientPoints; 350 PetscReal *cellPoints; 351 PetscReal *dummyWeights; 352 PetscQuadrature cellQuad = NULL; 353 354 ierr = DMFieldDSGetHeightDisc(field, 1, &faceDisc);CHKERRQ(ierr); 355 ierr = DMFieldDSGetHeightDisc(field, 0, &cellDisc);CHKERRQ(ierr); 356 ierr = PetscObjectGetClassId(faceDisc,&faceId);CHKERRQ(ierr); 357 ierr = PetscObjectGetClassId(cellDisc,&cellId);CHKERRQ(ierr); 358 if (faceId != PETSCFE_CLASSID || cellId != PETSCFE_CLASSID) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not supported\n"); 359 ierr = PetscFEGetDualSpace((PetscFE)cellDisc, &dsp);CHKERRQ(ierr); 360 ierr = PetscDualSpaceGetDM(dsp, &K); CHKERRQ(ierr); 361 ierr = DMPlexGetConeSize(K,0,&coneSize);CHKERRQ(ierr); 362 ierr = DMPlexGetCone(K,0,&coneK);CHKERRQ(ierr); 363 ierr = PetscMalloc4(numFaces,&co,dE*Nq,&cellPoints,coneSize,&counts,Nq,&dummyWeights);CHKERRQ(ierr); 364 ierr = PetscQuadratureCreate(PetscObjectComm((PetscObject)field), &cellQuad);CHKERRQ(ierr); 365 ierr = PetscQuadratureSetData(cellQuad, dE, 1, Nq, cellPoints, dummyWeights);CHKERRQ(ierr); 366 minOrient = PETSC_MAX_INT; 367 maxOrient = PETSC_MIN_INT; 368 for (p = 0; p < numFaces; p++) { /* record the orientation of the facet wrt the support cells */ 369 PetscInt point = points[p]; 370 PetscInt numSupp, numChildren; 371 const PetscInt *supp; 372 373 ierr = DMPlexGetTreeChildren(dm, point, &numChildren, NULL); CHKERRQ(ierr); 374 if (numChildren) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children"); 375 ierr = DMPlexGetSupportSize(dm, point,&numSupp);CHKERRQ(ierr); 376 if (numSupp > 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D has %D support, expected at most 2\n", point, numSupp); 377 ierr = DMPlexGetSupport(dm, point, &supp);CHKERRQ(ierr); 378 for (s = 0; s < numSupp; s++) { 379 PetscInt cell = supp[s]; 380 PetscInt numCone; 381 const PetscInt *cone, *orient; 382 383 ierr = DMPlexGetConeSize(dm, cell, &numCone);CHKERRQ(ierr); 384 if (numCone != coneSize) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support point does not match reference element"); 385 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 386 ierr = DMPlexGetConeOrientation(dm, cell, &orient);CHKERRQ(ierr); 387 for (f = 0; f < coneSize; f++) { 388 if (cone[f] == point) break; 389 } 390 co[p][s][0] = f; 391 co[p][s][1] = orient[f]; 392 co[p][s][2] = cell; 393 minOrient = PetscMin(minOrient, orient[f]); 394 maxOrient = PetscMin(maxOrient, orient[f]); 395 } 396 for (; s < 2; s++) { 397 co[p][s][0] = -1; 398 co[p][s][1] = -1; 399 co[p][s][2] = -1; 400 } 401 } 402 numOrient = maxOrient + 1 - minOrient; 403 ierr = DMPlexGetCone(K,0,&coneK);CHKERRQ(ierr); 404 /* count all (face,orientation) doubles that appear */ 405 ierr = PetscCalloc2(numOrient,&orients,numOrient,&orientPoints);CHKERRQ(ierr); 406 for (f = 0; f < coneSize; f++) {ierr = PetscCalloc1(numOrient, &counts[f]);CHKERRQ(ierr);} 407 for (p = 0; p < numFaces; p++) { 408 for (s = 0; s < 2; s++) { 409 if (co[p][s][0] >= 0) { 410 counts[co[p][s][0]][co[p][s][1] - minOrient]++; 411 orients[co[p][s][1] - minOrient]++; 412 } 413 } 414 } 415 for (o = 0; o < numOrient; o++) { 416 if (orients[o]) { 417 PetscInt orient = o + minOrient; 418 PetscInt q; 419 420 ierr = PetscMalloc1(Nq * dim, &orientPoints[o]);CHKERRQ(ierr); 421 /* rotate the quadrature points appropriately */ 422 switch (dim) { 423 case 0: 424 break; 425 case 1: 426 if (orient == -2 || orient == 1) { 427 for (q = 0; q < Nq; q++) { 428 orientPoints[o][q] = -geom->xi[q]; 429 } 430 } else { 431 for (q = 0; q < Nq; q++) { 432 orientPoints[o][q] = geom->xi[q]; 433 } 434 } 435 break; 436 case 2: 437 switch (coneSize) { 438 case 3: 439 for (q = 0; q < Nq; q++) { 440 PetscReal lambda[3]; 441 PetscReal lambdao[3]; 442 443 /* convert to barycentric */ 444 lambda[0] = - (geom->xi[2 * q] + geom->xi[2 * q + 1]) / 2.; 445 lambda[1] = (geom->xi[2 * q] + 1.) / 2.; 446 lambda[2] = (geom->xi[2 * q + 1] + 1.) / 2.; 447 if (orient >= 0) { 448 for (i = 0; i < 3; i++) { 449 lambdao[i] = lambda[(orient + i) % 3]; 450 } 451 } else { 452 for (i = 0; i < 3; i++) { 453 lambdao[i] = lambda[(-(orient + i) + 3) % 3]; 454 } 455 } 456 /* convert to coordinates */ 457 orientPoints[o][2 * q + 0] = -(lambdao[0] + lambdao[2]) + lambdao[1]; 458 orientPoints[o][2 * q + 1] = -(lambdao[0] + lambdao[1]) + lambdao[2]; 459 } 460 break; 461 case 4: 462 for (q = 0; q < Nq; q++) { 463 PetscReal xi[2], xio[2]; 464 PetscInt oabs = (orient >= 0) ? orient : -(orient + 1); 465 466 xi[0] = geom->xi[2 * q]; 467 xi[1] = geom->xi[2 * q + 1]; 468 switch (oabs) { 469 case 0: 470 xio[0] = xi[0]; 471 xio[1] = xi[1]; 472 break; 473 case 1: 474 xio[0] = xi[1]; 475 xio[1] = -xi[0]; 476 break; 477 case 2: 478 xio[0] = -xi[0]; 479 xio[1] = -xi[1]; 480 case 3: 481 xio[0] = -xi[1]; 482 xio[1] = xi[0]; 483 } 484 if (orient < 0) { 485 xio[0] = -xio[0]; 486 } 487 orientPoints[o][2 * q + 0] = xio[0]; 488 orientPoints[o][2 * q + 1] = xio[1]; 489 } 490 break; 491 default: 492 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not implemented yet\n"); 493 } 494 default: 495 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not implemented yet\n"); 496 } 497 } 498 } 499 for (f = 0; f < coneSize; f++) { 500 PetscInt face = coneK[f]; 501 PetscScalar v0[3]; 502 PetscScalar J[9]; 503 PetscInt numCells, offset; 504 PetscInt *cells; 505 IS suppIS; 506 507 ierr = DMPlexComputeCellGeometryFEM(K, face, NULL, v0, J, NULL, NULL);CHKERRQ(ierr); 508 for (o = 0; o <= numOrient; o++) { 509 PetscFEGeom *cellGeom; 510 511 if (!counts[f][o]) continue; 512 /* If this (face,orientation) double appears, 513 * convert the face quadrature points into volume quadrature points */ 514 for (q = 0; q < Nq; q++) { 515 PetscReal xi0[3] = {-1., -1., -1.}; 516 517 CoordinatesRefToReal(dE, dim, xi0, v0, J, &orientPoints[o][dim * q + 0], &cellPoints[dE * q + 0]); 518 } 519 for (p = 0, numCells = 0; p < numFaces; p++) { 520 for (s = 0; s < 2; s++) { 521 if (co[p][s][0] == f && co[p][s][1] == o + minOrient) numCells++; 522 } 523 } 524 ierr = PetscMalloc1(numCells, &cells);CHKERRQ(ierr); 525 for (p = 0, offset = 0; p < numFaces; p++) { 526 for (s = 0; s < 2; s++) { 527 if (co[p][s][0] == f && co[p][s][1] == o + minOrient) { 528 cells[offset++] = co[p][s][2]; 529 } 530 } 531 } 532 ierr = ISCreateGeneral(PETSC_COMM_SELF,numCells,cells,PETSC_USE_POINTER, &suppIS);CHKERRQ(ierr); 533 ierr = DMFieldCreateFEGeom(field,suppIS,cellQuad,PETSC_FALSE,&cellGeom);CHKERRQ(ierr); 534 for (p = 0, offset = 0; p < numFaces; p++) { 535 for (s = 0; s < 2; s++) { 536 if (co[p][s][0] == f && co[p][s][1] == o + minOrient) { 537 for (q = 0; q < Nq * dE * dE; q++) { 538 geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q]; 539 } 540 offset++; 541 } 542 } 543 } 544 ierr = PetscFEGeomDestroy(&cellGeom);CHKERRQ(ierr); 545 ierr = ISDestroy(&suppIS);CHKERRQ(ierr); 546 ierr = PetscFree(cells);CHKERRQ(ierr); 547 } 548 } 549 for (o = 0; o < numOrient; o++) { 550 if (orients[o]) { 551 ierr = PetscFree(orientPoints[o]);CHKERRQ(ierr); 552 } 553 } 554 ierr = PetscFree2(orients,orientPoints);CHKERRQ(ierr); 555 ierr = PetscQuadratureDestroy(&cellQuad);CHKERRQ(ierr); 556 ierr = PetscFree4(co,cellPoints,counts,dummyWeights);CHKERRQ(ierr); 557 } 558 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 559 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 560 PetscFunctionReturn(0); 561 } 562 563 static PetscErrorCode DMFieldInitialize_DS(DMField field) 564 { 565 PetscFunctionBegin; 566 field->ops->destroy = DMFieldDestroy_DS; 567 field->ops->evaluate = DMFieldEvaluate_DS; 568 field->ops->evaluateFE = DMFieldEvaluateFE_DS; 569 field->ops->getFEInvariance = DMFieldGetFEInvariance_DS; 570 field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS; 571 field->ops->view = DMFieldView_DS; 572 field->ops->computeFaceData = DMFieldComputeFaceData_DS; 573 PetscFunctionReturn(0); 574 } 575 576 PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field) 577 { 578 DMField_DS *dsfield; 579 PetscErrorCode ierr; 580 581 PetscFunctionBegin; 582 ierr = PetscNewLog(field,&dsfield);CHKERRQ(ierr); 583 field->data = dsfield; 584 ierr = DMFieldInitialize_DS(field);CHKERRQ(ierr); 585 PetscFunctionReturn(0); 586 } 587 588 PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec,DMField *field) 589 { 590 DMField b; 591 DMField_DS *dsfield; 592 PetscObject disc = NULL; 593 PetscBool isContainer = PETSC_FALSE; 594 PetscClassId id = -1; 595 PetscInt numComponents = -1, dsNumFields; 596 PetscSection section; 597 PetscErrorCode ierr; 598 599 PetscFunctionBegin; 600 ierr = DMGetDefaultSection(dm,§ion);CHKERRQ(ierr); 601 ierr = PetscSectionGetFieldComponents(section,fieldNum,&numComponents);CHKERRQ(ierr); 602 ierr = DMGetNumFields(dm,&dsNumFields);CHKERRQ(ierr); 603 if (dsNumFields) {ierr = DMGetField(dm,fieldNum,&disc);CHKERRQ(ierr);} 604 if (disc) { 605 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 606 isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE; 607 } 608 if (!disc || isContainer) { 609 MPI_Comm comm = PetscObjectComm((PetscObject) dm); 610 PetscInt cStart, cEnd, dim; 611 PetscInt localConeSize = 0, coneSize; 612 PetscFE fe; 613 PetscDualSpace Q; 614 PetscSpace P; 615 DM K; 616 PetscQuadrature quad, fquad; 617 PetscBool isSimplex; 618 619 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 620 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 621 if (cEnd > cStart) { 622 ierr = DMPlexGetConeSize(dm, cStart, &localConeSize);CHKERRQ(ierr); 623 } 624 ierr = MPI_Allreduce(&localConeSize,&coneSize,1,MPIU_INT,MPI_MAX,comm);CHKERRQ(ierr); 625 isSimplex = (coneSize == (dim + 1)) ? PETSC_TRUE : PETSC_FALSE; 626 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 627 ierr = PetscSpaceSetOrder(P, 1);CHKERRQ(ierr); 628 ierr = PetscSpaceSetNumComponents(P, numComponents);CHKERRQ(ierr); 629 ierr = PetscSpaceSetType(P,PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 630 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 631 ierr = PetscSpacePolynomialSetTensor(P, isSimplex ? PETSC_FALSE : PETSC_TRUE);CHKERRQ(ierr); 632 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 633 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 634 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 635 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 636 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 637 ierr = DMDestroy(&K);CHKERRQ(ierr); 638 ierr = PetscDualSpaceSetNumComponents(Q, numComponents);CHKERRQ(ierr); 639 ierr = PetscDualSpaceSetOrder(Q, 1);CHKERRQ(ierr); 640 ierr = PetscDualSpaceLagrangeSetTensor(Q, isSimplex ? PETSC_FALSE : PETSC_TRUE);CHKERRQ(ierr); 641 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 642 ierr = PetscFECreate(comm, &fe);CHKERRQ(ierr); 643 ierr = PetscFESetType(fe,PETSCFEBASIC);CHKERRQ(ierr); 644 ierr = PetscFESetBasisSpace(fe, P);CHKERRQ(ierr); 645 ierr = PetscFESetDualSpace(fe, Q);CHKERRQ(ierr); 646 ierr = PetscFESetNumComponents(fe, numComponents);CHKERRQ(ierr); 647 ierr = PetscFESetUp(fe);CHKERRQ(ierr); 648 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 649 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 650 if (isSimplex) { 651 ierr = PetscDTGaussJacobiQuadrature(dim, 1, 1, -1.0, 1.0, &quad);CHKERRQ(ierr); 652 ierr = PetscDTGaussJacobiQuadrature(dim-1, 1, 1, -1.0, 1.0, &fquad);CHKERRQ(ierr); 653 } 654 else { 655 ierr = PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, &quad);CHKERRQ(ierr); 656 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, 1, -1.0, 1.0, &fquad);CHKERRQ(ierr); 657 } 658 ierr = PetscFESetQuadrature(fe, quad);CHKERRQ(ierr); 659 ierr = PetscFESetFaceQuadrature(fe, fquad);CHKERRQ(ierr); 660 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 661 ierr = PetscQuadratureDestroy(&fquad);CHKERRQ(ierr); 662 disc = (PetscObject) fe; 663 } else { 664 ierr = PetscObjectReference(disc);CHKERRQ(ierr); 665 } 666 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 667 if (id == PETSCFE_CLASSID) { 668 PetscFE fe = (PetscFE) disc; 669 670 ierr = PetscFEGetNumComponents(fe,&numComponents);CHKERRQ(ierr); 671 } else {SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented");} 672 ierr = DMFieldCreate(dm,numComponents,DMFIELD_VERTEX,&b);CHKERRQ(ierr); 673 ierr = DMFieldSetType(b,DMFIELDDS);CHKERRQ(ierr); 674 dsfield = (DMField_DS *) b->data; 675 dsfield->fieldNum = fieldNum; 676 ierr = DMGetDimension(dm,&dsfield->height);CHKERRQ(ierr); 677 dsfield->height++; 678 ierr = PetscCalloc1(dsfield->height,&dsfield->disc);CHKERRQ(ierr); 679 dsfield->disc[0] = disc; 680 ierr = PetscObjectReference((PetscObject)vec);CHKERRQ(ierr); 681 dsfield->vec = vec; 682 *field = b; 683 684 PetscFunctionReturn(0); 685 } 686