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