1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #include <petsc/private/petscfeimpl.h> 4 5 static PetscErrorCode DMProjectPoint_Func_Private(DM dm, PetscReal time, PetscFECellGeom *fegeom, PetscFVCellGeom *fvgeom, PetscBool isFE[], PetscDualSpace sp[], 6 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, 7 PetscScalar values[]) 8 { 9 PetscDS prob; 10 PetscInt Nf, *Nc, f, totDim, spDim, d, v; 11 PetscErrorCode ierr; 12 13 PetscFunctionBeginHot; 14 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 15 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 16 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 17 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 18 /* Get values for closure */ 19 for (f = 0, v = 0; f < Nf; ++f) { 20 void * const ctx = ctxs ? ctxs[f] : NULL; 21 22 if (!sp[f]) continue; 23 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 24 for (d = 0; d < spDim; ++d, ++v) { 25 if (funcs[f]) { 26 if (isFE[f]) {ierr = PetscDualSpaceApply(sp[f], d, time, fegeom, Nc[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);} 27 else {ierr = PetscDualSpaceApplyFVM(sp[f], d, time, fvgeom, Nc[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);} 28 } else { 29 values[v] = 0.0; 30 } 31 } 32 } 33 PetscFunctionReturn(0); 34 } 35 36 static PetscErrorCode DMProjectPoint_Field_Private(DM dm, DM dmAux, PetscReal time, Vec localU, Vec localA, PetscFECellGeom *fegeom, PetscDualSpace sp[], PetscInt p, PetscInt Ncc, const PetscInt comps[], 37 PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux, 38 void (**funcs)(PetscInt, PetscInt, PetscInt, 39 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 40 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 41 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void **ctxs, 42 PetscScalar values[]) 43 { 44 PetscDS prob, probAux = NULL; 45 PetscSection section, sectionAux = NULL; 46 PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL, *refSpaceDer, *refSpaceDerAux = NULL, *bc; 47 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 48 PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL; 49 const PetscScalar *constants; 50 PetscReal *x; 51 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 52 PetscInt dim, dimAux = 0, numConstants, Nf, NfAux = 0, f, spDim, d, c, v, tp = 0; 53 PetscErrorCode ierr; 54 55 PetscFunctionBeginHot; 56 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 57 ierr = PetscDSGetSpatialDimension(prob, &dim);CHKERRQ(ierr); 58 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 59 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 60 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 61 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 62 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 63 ierr = PetscDSGetEvaluationArrays(prob, &u, &bc /*&u_t*/, &u_x);CHKERRQ(ierr); 64 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 65 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 66 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 67 ierr = DMPlexVecGetClosure(dm, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 68 if (dmAux) { 69 DMLabel spmap; 70 PetscInt subp; 71 72 ierr = DMPlexGetSubpointMap(dmAux, &spmap);CHKERRQ(ierr); 73 if (spmap) { 74 IS subpointIS; 75 const PetscInt *subpoints; 76 PetscInt numSubpoints; 77 78 ierr = DMPlexCreateSubpointIS(dmAux, &subpointIS);CHKERRQ(ierr); 79 ierr = ISGetLocalSize(subpointIS, &numSubpoints);CHKERRQ(ierr); 80 ierr = ISGetIndices(subpointIS, &subpoints);CHKERRQ(ierr); 81 ierr = PetscFindInt(p, numSubpoints, subpoints, &subp);CHKERRQ(ierr); 82 if (subp < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", p); 83 ierr = ISRestoreIndices(subpointIS, &subpoints);CHKERRQ(ierr); 84 ierr = ISDestroy(&subpointIS);CHKERRQ(ierr); 85 } else subp = p; 86 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 87 ierr = PetscDSGetSpatialDimension(probAux, &dimAux);CHKERRQ(ierr); 88 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 89 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 90 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 91 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 92 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 93 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 94 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL /*&a_t*/, &a_x);CHKERRQ(ierr); 95 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 96 ierr = DMPlexVecGetClosure(dmAux, sectionAux, localA, subp, NULL, &coefficientsAux);CHKERRQ(ierr); 97 } 98 /* Get values for closure */ 99 for (f = 0, v = 0; f < Nf; ++f) { 100 if (!sp[f]) continue; 101 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 102 for (d = 0; d < spDim; ++d, ++v) { 103 if (funcs[f]) { 104 PetscQuadrature quad; 105 const PetscReal *points, *weights; 106 PetscInt numPoints, q; 107 108 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 109 ierr = PetscQuadratureGetData(quad, NULL, NULL, &numPoints, &points, &weights);CHKERRQ(ierr); 110 for (q = 0; q < numPoints; ++q, ++tp) { 111 CoordinatesRefToReal(dim, dim, fegeom->v0, fegeom->J, &points[q*dim], x); 112 EvaluateFieldJets(dim, Nf, Nb, Nc, tp, basisTab, basisDerTab, refSpaceDer, fegeom->invJ, coefficients, coefficients_t, u, u_x, u_t); 113 if (probAux) {EvaluateFieldJets(dimAux, NfAux, NbAux, NcAux, tp, basisTabAux, basisDerTabAux, refSpaceDerAux, fegeom->invJ, coefficientsAux, coefficientsAux_t, a, a_x, a_t);} 114 (*funcs[f])(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, time, x, numConstants, constants, bc); 115 if (Ncc) { 116 for (c = 0; c < Ncc; ++c) values[v] += bc[comps[c]]*weights[comps[c]]; 117 } else { 118 for (c = 0; c < Nc[f]; ++c) values[v] += bc[c]*weights[c]; 119 } 120 } 121 } else { 122 values[v] = 0.0; 123 } 124 } 125 } 126 ierr = DMPlexVecRestoreClosure(dm, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 127 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);CHKERRQ(ierr);} 128 PetscFunctionReturn(0); 129 } 130 131 static PetscErrorCode DMProjectPoint_Private(DM dm, DM dmAux, PetscInt h, PetscReal time, Vec localU, Vec localA, PetscBool hasFE, PetscBool hasFV, PetscBool isFE[], 132 PetscDualSpace sp[], PetscInt p, PetscInt Ncc, const PetscInt comps[], 133 PetscReal **basisTab, PetscReal **basisDerTab, PetscReal **basisTabAux, PetscReal **basisDerTabAux, 134 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[]) 135 { 136 PetscFECellGeom fegeom; 137 PetscFVCellGeom fvgeom; 138 PetscInt dim, dimEmbed; 139 PetscErrorCode ierr; 140 141 PetscFunctionBeginHot; 142 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 143 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 144 if (hasFE) { 145 ierr = DMPlexComputeCellGeometryFEM(dm, p, NULL, fegeom.v0, fegeom.J, fegeom.invJ, &fegeom.detJ);CHKERRQ(ierr); 146 fegeom.dim = dim - h; 147 fegeom.dimEmbed = dimEmbed; 148 } 149 if (hasFV) {ierr = DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);CHKERRQ(ierr);} 150 switch (type) { 151 case DM_BC_ESSENTIAL: 152 case DM_BC_NATURAL: 153 ierr = DMProjectPoint_Func_Private(dm, time, &fegeom, &fvgeom, isFE, sp, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);CHKERRQ(ierr);break; 154 case DM_BC_ESSENTIAL_FIELD: 155 case DM_BC_NATURAL_FIELD: 156 ierr = DMProjectPoint_Field_Private(dm, dmAux, time, localU, localA, &fegeom, sp, p, Ncc, comps, 157 basisTab, basisDerTab, basisTabAux, basisDerTabAux, 158 (void (**)(PetscInt, PetscInt, PetscInt, 159 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 160 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 161 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);CHKERRQ(ierr);break; 162 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type); 163 } 164 PetscFunctionReturn(0); 165 } 166 167 static PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU, 168 PetscInt Ncc, const PetscInt comps[], DMLabel label, PetscInt numIds, const PetscInt ids[], 169 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, 170 InsertMode mode, Vec localX) 171 { 172 DM dmAux = NULL; 173 PetscDS prob, probAux = NULL; 174 Vec localA = NULL; 175 PetscSection section; 176 PetscDualSpace *sp, *cellsp; 177 PetscReal **basisTab = NULL, **basisDerTab = NULL, **basisTabAux = NULL, **basisDerTabAux = NULL; 178 PetscInt *Nc; 179 PetscInt dim, dimEmbed, minHeight, maxHeight, h, Nf, NfAux = 0, f; 180 PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE, auxBd = PETSC_FALSE; 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 185 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);CHKERRQ(ierr); 186 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 187 ierr = DMPlexGetVTKCellHeight(dm, &minHeight);CHKERRQ(ierr); 188 if (!minHeight && dmAux) { 189 DMLabel spmap; 190 191 /* If dmAux is a surface, then force the projection to take place over a surface */ 192 ierr = DMPlexGetSubpointMap(dmAux, &spmap);CHKERRQ(ierr); 193 if (spmap) {minHeight = 1; auxBd = PETSC_TRUE;} 194 } 195 ierr = DMPlexGetMaxProjectionHeight(dm, &maxHeight);CHKERRQ(ierr); 196 maxHeight = PetscMax(maxHeight, minHeight); 197 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %d not in [0, %d)\n", maxHeight, dim);} 198 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 199 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 200 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 201 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 202 if (dmAux) { 203 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 204 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 205 } 206 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 207 ierr = PetscMalloc2(Nf, &isFE, Nf, &sp);CHKERRQ(ierr); 208 if (maxHeight > 0) {ierr = PetscMalloc1(Nf, &cellsp);CHKERRQ(ierr);} 209 else {cellsp = sp;} 210 if (localU && localU != localX) {ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, time, NULL, NULL, NULL);CHKERRQ(ierr);} 211 /* Get cell dual spaces */ 212 for (f = 0; f < Nf; ++f) { 213 PetscObject obj; 214 PetscClassId id; 215 216 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 217 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 218 if (id == PETSCFE_CLASSID) { 219 PetscFE fe = (PetscFE) obj; 220 221 hasFE = PETSC_TRUE; 222 isFE[f] = PETSC_TRUE; 223 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 224 } else if (id == PETSCFV_CLASSID) { 225 PetscFV fv = (PetscFV) obj; 226 227 hasFV = PETSC_TRUE; 228 isFE[f] = PETSC_FALSE; 229 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 230 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 231 } 232 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 233 PetscInt effectiveHeight = auxBd ? minHeight : 0; 234 PetscFE fem, subfem; 235 PetscReal *points; 236 PetscInt numPoints, spDim, qdim = 0, d; 237 238 if (maxHeight > minHeight) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Field projection not supported for face interpolation"); 239 numPoints = 0; 240 for (f = 0; f < Nf; ++f) { 241 if (!effectiveHeight) {sp[f] = cellsp[f];} 242 else {ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr);} 243 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 244 for (d = 0; d < spDim; ++d) { 245 if (funcs[f]) { 246 PetscQuadrature quad; 247 PetscInt Nq; 248 249 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 250 ierr = PetscQuadratureGetData(quad, &qdim, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 251 numPoints += Nq; 252 } 253 } 254 } 255 ierr = PetscMalloc1(numPoints*qdim, &points);CHKERRQ(ierr); 256 numPoints = 0; 257 for (f = 0; f < Nf; ++f) { 258 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 259 for (d = 0; d < spDim; ++d) { 260 if (funcs[f]) { 261 PetscQuadrature quad; 262 const PetscReal *qpoints; 263 PetscInt Nq, q; 264 265 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 266 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, NULL);CHKERRQ(ierr); 267 for (q = 0; q < Nq*qdim; ++q) points[numPoints*qdim+q] = qpoints[q]; 268 numPoints += Nq; 269 } 270 } 271 } 272 ierr = PetscMalloc4(Nf, &basisTab, Nf, &basisDerTab, NfAux, &basisTabAux, NfAux, &basisDerTabAux);CHKERRQ(ierr); 273 for (f = 0; f < Nf; ++f) { 274 if (!isFE[f]) continue; 275 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 276 if (!effectiveHeight) {subfem = fem;} 277 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 278 ierr = PetscFEGetTabulation(subfem, numPoints, points, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 279 } 280 for (f = 0; f < NfAux; ++f) { 281 ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 282 if (!effectiveHeight || auxBd) {subfem = fem;} 283 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 284 ierr = PetscFEGetTabulation(subfem, numPoints, points, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 285 } 286 ierr = PetscFree(points);CHKERRQ(ierr); 287 } 288 /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */ 289 for (h = minHeight; h <= maxHeight; h++) { 290 PetscInt effectiveHeight = h - (auxBd ? 0 : minHeight); 291 PetscScalar *values; 292 PetscBool *fieldActive; 293 PetscInt pStart, pEnd, p, spDim, totDim, numValues; 294 295 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 296 if (!h) { 297 PetscInt cEndInterior; 298 299 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 300 pEnd = cEndInterior < 0 ? pEnd : cEndInterior; 301 } 302 if (pEnd <= pStart) continue; 303 /* Compute totDim, the number of dofs in the closure of a point at this height */ 304 totDim = 0; 305 for (f = 0; f < Nf; ++f) { 306 if (!effectiveHeight) { 307 sp[f] = cellsp[f]; 308 } else { 309 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], effectiveHeight, &sp[f]);CHKERRQ(ierr); 310 if (!sp[f]) continue; 311 } 312 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 313 totDim += spDim; 314 } 315 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 316 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim); 317 if (!totDim) continue; 318 /* Loop over points at this height */ 319 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 320 ierr = DMGetWorkArray(dm, Nf, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 321 for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 322 if (label) { 323 PetscInt i; 324 325 for (i = 0; i < numIds; ++i) { 326 IS pointIS; 327 const PetscInt *points; 328 PetscInt n; 329 330 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 331 if (!pointIS) continue; /* No points with that id on this process */ 332 ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 333 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 334 for (p = 0; p < n; ++p) { 335 const PetscInt point = points[p]; 336 337 if ((point < pStart) || (point >= pEnd)) continue; 338 ierr = PetscMemzero(values, numValues * sizeof(PetscScalar));CHKERRQ(ierr); 339 ierr = DMProjectPoint_Private(dm, dmAux, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, point, Ncc, comps, basisTab, basisDerTab, basisTabAux, basisDerTabAux, type, funcs, ctxs, fieldActive, values); 340 if (ierr) { 341 PetscErrorCode ierr2; 342 ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2); 343 ierr2 = DMRestoreWorkArray(dm, Nf, PETSC_BOOL, &fieldActive);CHKERRQ(ierr2); 344 CHKERRQ(ierr); 345 } 346 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);CHKERRQ(ierr); 347 } 348 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 349 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 350 } 351 } else { 352 for (p = pStart; p < pEnd; ++p) { 353 ierr = PetscMemzero(values, numValues * sizeof(PetscScalar));CHKERRQ(ierr); 354 ierr = DMProjectPoint_Private(dm, dmAux, effectiveHeight, time, localU, localA, hasFE, hasFV, isFE, sp, p, Ncc, comps, basisTab, basisDerTab, basisTabAux, basisDerTabAux, type, funcs, ctxs, fieldActive, values); 355 if (ierr) { 356 PetscErrorCode ierr2; 357 ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2); 358 ierr2 = DMRestoreWorkArray(dm, Nf, PETSC_BOOL, &fieldActive);CHKERRQ(ierr2); 359 CHKERRQ(ierr); 360 } 361 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, p, values, mode);CHKERRQ(ierr); 362 } 363 } 364 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 365 ierr = DMRestoreWorkArray(dm, Nf, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 366 } 367 /* Cleanup */ 368 if (type == DM_BC_ESSENTIAL_FIELD || type == DM_BC_NATURAL_FIELD) { 369 PetscInt effectiveHeight = auxBd ? minHeight : 0; 370 PetscFE fem, subfem; 371 372 for (f = 0; f < Nf; ++f) { 373 if (!isFE[f]) continue; 374 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fem);CHKERRQ(ierr); 375 if (!effectiveHeight) {subfem = fem;} 376 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 377 ierr = PetscFERestoreTabulation(subfem, 0, NULL, &basisTab[f], &basisDerTab[f], NULL);CHKERRQ(ierr); 378 } 379 for (f = 0; f < NfAux; ++f) { 380 ierr = PetscDSGetDiscretization(probAux, f, (PetscObject *) &fem);CHKERRQ(ierr); 381 if (!effectiveHeight || auxBd) {subfem = fem;} 382 else {ierr = PetscFEGetHeightSubspace(fem, effectiveHeight, &subfem);CHKERRQ(ierr);} 383 ierr = PetscFERestoreTabulation(subfem, 0, NULL, &basisTabAux[f], &basisDerTabAux[f], NULL);CHKERRQ(ierr); 384 } 385 ierr = PetscFree4(basisTab, basisDerTab, basisTabAux, basisDerTabAux);CHKERRQ(ierr); 386 } 387 ierr = PetscFree2(isFE, sp);CHKERRQ(ierr); 388 if (maxHeight > 0) {ierr = PetscFree(cellsp);CHKERRQ(ierr);} 389 PetscFunctionReturn(0); 390 } 391 392 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 393 { 394 PetscErrorCode ierr; 395 396 PetscFunctionBegin; 397 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 398 PetscFunctionReturn(0); 399 } 400 401 PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 402 { 403 PetscErrorCode ierr; 404 405 PetscFunctionBegin; 406 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 407 PetscFunctionReturn(0); 408 } 409 410 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU, 411 void (**funcs)(PetscInt, PetscInt, PetscInt, 412 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 413 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 414 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 415 InsertMode mode, Vec localX) 416 { 417 PetscErrorCode ierr; 418 419 PetscFunctionBegin; 420 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, 0, NULL, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 421 PetscFunctionReturn(0); 422 } 423 424 PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Ncc, const PetscInt comps[], Vec localU, 425 void (**funcs)(PetscInt, PetscInt, PetscInt, 426 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 427 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 428 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 429 InsertMode mode, Vec localX) 430 { 431 PetscErrorCode ierr; 432 433 PetscFunctionBegin; 434 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, Ncc, comps, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 435 PetscFunctionReturn(0); 436 } 437