1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #include <petsc/private/petscfeimpl.h> 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "DMProjectPoint_Func_Private" 7 static PetscErrorCode DMProjectPoint_Func_Private(DM dm, PetscReal time, PetscFECellGeom *fegeom, PetscFVCellGeom *fvgeom, PetscBool isFE[], PetscDualSpace sp[], PetscInt Nc[], 8 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, 9 PetscScalar values[]) 10 { 11 PetscDS prob; 12 PetscInt Nf, f, spDim, d, c, v; 13 PetscErrorCode ierr; 14 15 PetscFunctionBeginHot; 16 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 17 ierr = PetscDSGetNumFields(prob, &Nf);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) { 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 for (c = 0; c < Nc[f]; ++c) values[v+c] = 0.0; 30 } 31 v += Nc[f]; 32 } 33 } 34 PetscFunctionReturn(0); 35 } 36 37 #undef __FUNCT__ 38 #define __FUNCT__ "DMProjectPoint_Field_Private" 39 static PetscErrorCode DMProjectPoint_Field_Private(DM dm, DM dmAux, PetscReal time, Vec localU, Vec localA, PetscFECellGeom *fegeom, PetscDualSpace sp[], PetscInt Nc[], PetscInt p, 40 void (**funcs)(PetscInt, PetscInt, PetscInt, 41 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 42 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 43 PetscReal, const PetscReal[], PetscScalar[]), void **ctxs, 44 PetscScalar values[]) 45 { 46 PetscDS prob, probAux = NULL; 47 PetscSection section, sectionAux = NULL; 48 PetscScalar *u, *u_t = NULL, *u_x, *a = NULL, *a_t = NULL, *a_x = NULL; 49 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 50 PetscScalar *coefficients_t = NULL, *coefficientsAux_t = NULL; 51 PetscReal *x; 52 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 53 PetscInt dim, Nf, NfAux = 0, f, spDim, d, c, v; 54 PetscErrorCode ierr; 55 56 PetscFunctionBeginHot; 57 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 58 ierr = PetscDSGetSpatialDimension(prob, &dim);CHKERRQ(ierr); 59 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 60 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 61 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 62 ierr = PetscDSGetEvaluationArrays(prob, &u, NULL /*&u_t*/, &u_x);CHKERRQ(ierr); 63 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 64 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 65 ierr = DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 66 if (dmAux) { 67 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 68 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 69 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 70 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 71 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 72 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL /*&a_t*/, &a_x);CHKERRQ(ierr); 73 ierr = DMPlexVecGetClosure(dmAux, sectionAux, localA, c, NULL, &coefficientsAux);CHKERRQ(ierr); 74 } 75 /* Get values for closure */ 76 for (f = 0, v = 0; f < Nf; ++f) { 77 void * const ctx = ctxs ? ctxs[f] : NULL; 78 79 if (!sp[f]) continue; 80 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 81 for (d = 0; d < spDim; ++d) { 82 if (funcs[f]) { 83 PetscQuadrature quad; 84 const PetscReal *points, *weights; 85 PetscInt numPoints, q; 86 87 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 88 ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);CHKERRQ(ierr); 89 for (q = 0; q < numPoints; ++q) { 90 CoordinatesRefToReal(dim, dim, fegeom->v0, fegeom->J, &points[q*dim], x); 91 ierr = EvaluateFieldJets(prob, PETSC_FALSE, q, fegeom->invJ, coefficients, coefficients_t, u, u_x, u_t);CHKERRQ(ierr); 92 ierr = EvaluateFieldJets(probAux, PETSC_FALSE, q, fegeom->invJ, coefficientsAux, coefficientsAux_t, a, a_x, a_t);CHKERRQ(ierr); 93 (*funcs[f])(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, time, x, &values[v]); 94 } 95 } else { 96 for (c = 0; c < Nc[f]; ++c) values[v+c] = 0.0; 97 } 98 v += Nc[f]; 99 } 100 } 101 ierr = DMPlexVecRestoreClosure(dm, section, localU, p, NULL, &coefficients);CHKERRQ(ierr); 102 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, localA, p, NULL, &coefficientsAux);CHKERRQ(ierr);} 103 PetscFunctionReturn(0); 104 } 105 106 #undef __FUNCT__ 107 #define __FUNCT__ "DMProjectPoint_Private" 108 static PetscErrorCode DMProjectPoint_Private(DM dm, DM dmAux, PetscInt h, PetscReal time, Vec localU, Vec localA, PetscBool hasFE, PetscBool hasFV, PetscBool isFE[], PetscDualSpace sp[], PetscInt Nc[], PetscInt p, 109 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, PetscBool fieldActive[], PetscScalar values[]) 110 { 111 PetscFECellGeom fegeom; 112 PetscFVCellGeom fvgeom; 113 PetscInt dim, dimEmbed; 114 PetscErrorCode ierr; 115 116 PetscFunctionBeginHot; 117 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 118 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 119 if (hasFE) { 120 ierr = DMPlexComputeCellGeometryFEM(dm, p, NULL, fegeom.v0, fegeom.J, NULL, &fegeom.detJ);CHKERRQ(ierr); 121 fegeom.dim = dim - h; 122 fegeom.dimEmbed = dimEmbed; 123 } 124 if (hasFV) {ierr = DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);CHKERRQ(ierr);} 125 switch (type) { 126 case DM_BC_ESSENTIAL: 127 case DM_BC_NATURAL: 128 ierr = DMProjectPoint_Func_Private(dm, time, &fegeom, &fvgeom, isFE, sp, Nc, (PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *)) funcs, ctxs, values);CHKERRQ(ierr);break; 129 case DM_BC_ESSENTIAL_FIELD: 130 case DM_BC_NATURAL_FIELD: 131 ierr = DMProjectPoint_Field_Private(dm, dmAux, time, localU, localA, &fegeom, sp, Nc, p, (void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscScalar[])) funcs, ctxs, values);CHKERRQ(ierr);break; 132 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown boundary condition type: %d", (int) type); 133 } 134 PetscFunctionReturn(0); 135 } 136 137 #undef __FUNCT__ 138 #define __FUNCT__ "DMProjectLocal_Generic_Plex" 139 PetscErrorCode DMProjectLocal_Generic_Plex(DM dm, PetscReal time, Vec localU, 140 DMLabel label, PetscInt numIds, const PetscInt ids[], 141 DMBoundaryConditionType type, void (**funcs)(void), void **ctxs, 142 InsertMode mode, Vec localX) 143 { 144 DM dmAux = NULL; 145 Vec localA = NULL; 146 PetscSection section; 147 PetscDualSpace *sp, *cellsp; 148 PetscInt *Nc; 149 PetscInt dim, dimEmbed, maxHeight, h, Nf, f; 150 PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE; 151 PetscErrorCode ierr; 152 153 PetscFunctionBegin; 154 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 155 ierr = DMPlexGetMaxProjectionHeight(dm, &maxHeight);CHKERRQ(ierr); 156 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Maximum projection height %d not in [0, %d)\n", maxHeight, dim);} 157 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 158 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 159 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 160 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 161 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &localA);CHKERRQ(ierr); 162 ierr = PetscMalloc3(Nf, &isFE, Nf, &sp, Nf, &Nc);CHKERRQ(ierr); 163 if (maxHeight > 0) {ierr = PetscMalloc1(Nf, &cellsp);CHKERRQ(ierr);} 164 else {cellsp = sp;} 165 if (localU && localU != localX) {ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, time, NULL, NULL, NULL);CHKERRQ(ierr);} 166 /* Get cell dual spaces */ 167 for (f = 0; f < Nf; ++f) { 168 PetscObject obj; 169 PetscClassId id; 170 171 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 172 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 173 if (id == PETSCFE_CLASSID) { 174 PetscFE fe = (PetscFE) obj; 175 176 hasFE = PETSC_TRUE; 177 isFE[f] = PETSC_TRUE; 178 ierr = PetscFEGetNumComponents(fe, &Nc[f]);CHKERRQ(ierr); 179 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 180 } else if (id == PETSCFV_CLASSID) { 181 PetscFV fv = (PetscFV) obj; 182 183 hasFV = PETSC_TRUE; 184 isFE[f] = PETSC_FALSE; 185 ierr = PetscFVGetNumComponents(fv, &Nc[f]);CHKERRQ(ierr); 186 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 187 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 188 } 189 /* Note: We make no attempt to optimize for height. Higher height things just overwrite the lower height results. */ 190 for (h = 0; h <= maxHeight; h++) { 191 PetscScalar *values; 192 PetscBool *fieldActive; 193 PetscInt pStart, pEnd, p, spDim, totDim, d, numValues, v; 194 195 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 196 if (!h) { 197 PetscInt cEndInterior; 198 199 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 200 pEnd = cEndInterior < 0 ? pEnd : cEndInterior; 201 } 202 if (pEnd <= pStart) continue; 203 /* Compute totDim, the number of dofs in the closure of a point at this height */ 204 totDim = 0; 205 for (f = 0; f < Nf; ++f) { 206 if (!h) { 207 sp[f] = cellsp[f]; 208 } else { 209 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr); 210 if (!sp[f]) continue; 211 } 212 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 213 totDim += spDim*Nc[f]; 214 } 215 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 216 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The section point closure size %d != dual space dimension %d", numValues, totDim); 217 if (!totDim) continue; 218 /* Loop over points at this height */ 219 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 220 ierr = DMGetWorkArray(dm, Nf, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 221 for (f = 0; f < Nf; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 222 if (label) { 223 PetscInt i; 224 225 for (i = 0; i < numIds; ++i) { 226 IS pointIS; 227 const PetscInt *points; 228 PetscInt n; 229 230 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 231 if (!pointIS) continue; /* No points with that id on this process */ 232 ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 233 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 234 for (p = 0; p < n; ++p) { 235 const PetscInt point = points[p]; 236 237 if ((point < pStart) || (point >= pEnd)) continue; 238 ierr = DMProjectPoint_Private(dm, dmAux, h, time, localU, localA, hasFE, hasFV, isFE, sp, Nc, point, type, funcs, ctxs, fieldActive, values); 239 if (ierr) { 240 PetscErrorCode ierr2; 241 ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2); 242 ierr2 = DMRestoreWorkArray(dm, Nf, PETSC_BOOL, &fieldActive);CHKERRQ(ierr2); 243 CHKERRQ(ierr); 244 } 245 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);CHKERRQ(ierr); 246 } 247 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 248 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 249 } 250 } else { 251 for (p = pStart; p < pEnd; ++p) { 252 ierr = DMProjectPoint_Private(dm, dmAux, h, time, localU, localA, hasFE, hasFV, isFE, sp, Nc, p, type, funcs, ctxs, fieldActive, values); 253 if (ierr) { 254 PetscErrorCode ierr2; 255 ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2); 256 ierr2 = DMRestoreWorkArray(dm, Nf, PETSC_BOOL, &fieldActive);CHKERRQ(ierr2); 257 CHKERRQ(ierr); 258 } 259 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, p, values, mode);CHKERRQ(ierr); 260 } 261 } 262 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 263 ierr = DMRestoreWorkArray(dm, Nf, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 264 } 265 /* Cleanup */ 266 ierr = PetscFree3(isFE, sp, Nc);CHKERRQ(ierr); 267 if (maxHeight > 0) {ierr = PetscFree(cellsp);CHKERRQ(ierr);} 268 PetscFunctionReturn(0); 269 } 270 271 #undef __FUNCT__ 272 #define __FUNCT__ "DMProjectFunctionLocal_Plex" 273 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 274 { 275 PetscErrorCode ierr; 276 277 PetscFunctionBegin; 278 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, NULL, 0, NULL, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 279 PetscFunctionReturn(0); 280 } 281 282 #undef __FUNCT__ 283 #define __FUNCT__ "DMProjectFunctionLabelLocal_Plex" 284 PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 285 { 286 PetscErrorCode ierr; 287 288 PetscFunctionBegin; 289 ierr = DMProjectLocal_Generic_Plex(dm, time, localX, label, numIds, ids, DM_BC_ESSENTIAL, (void (**)(void)) funcs, ctxs, mode, localX);CHKERRQ(ierr); 290 PetscFunctionReturn(0); 291 } 292 293 #undef __FUNCT__ 294 #define __FUNCT__ "DMProjectFieldLocal_Plex" 295 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, PetscReal time, Vec localU, 296 void (**funcs)(PetscInt, PetscInt, PetscInt, 297 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 298 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 299 PetscReal, const PetscReal[], PetscScalar[]), 300 InsertMode mode, Vec localX) 301 { 302 PetscErrorCode ierr; 303 304 PetscFunctionBegin; 305 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, NULL, 0, NULL, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 306 PetscFunctionReturn(0); 307 } 308 309 #undef __FUNCT__ 310 #define __FUNCT__ "DMProjectFieldLabelLocal_Plex" 311 PetscErrorCode DMProjectFieldLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], Vec localU, 312 void (**funcs)(PetscInt, PetscInt, PetscInt, 313 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 314 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 315 PetscReal, const PetscReal[], PetscScalar[]), 316 InsertMode mode, Vec localX) 317 { 318 PetscErrorCode ierr; 319 320 PetscFunctionBegin; 321 ierr = DMProjectLocal_Generic_Plex(dm, time, localU, label, numIds, ids, DM_BC_ESSENTIAL_FIELD, (void (**)(void)) funcs, NULL, mode, localX);CHKERRQ(ierr); 322 PetscFunctionReturn(0); 323 } 324