1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 4 #include <petsc/private/petscfeimpl.h> 5 #include <petsc/private/petscfvimpl.h> 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "DMPlexGetScale" 9 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 10 { 11 DM_Plex *mesh = (DM_Plex*) dm->data; 12 13 PetscFunctionBegin; 14 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 15 PetscValidPointer(scale, 3); 16 *scale = mesh->scale[unit]; 17 PetscFunctionReturn(0); 18 } 19 20 #undef __FUNCT__ 21 #define __FUNCT__ "DMPlexSetScale" 22 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 23 { 24 DM_Plex *mesh = (DM_Plex*) dm->data; 25 26 PetscFunctionBegin; 27 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 28 mesh->scale[unit] = scale; 29 PetscFunctionReturn(0); 30 } 31 32 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k) 33 { 34 switch (i) { 35 case 0: 36 switch (j) { 37 case 0: return 0; 38 case 1: 39 switch (k) { 40 case 0: return 0; 41 case 1: return 0; 42 case 2: return 1; 43 } 44 case 2: 45 switch (k) { 46 case 0: return 0; 47 case 1: return -1; 48 case 2: return 0; 49 } 50 } 51 case 1: 52 switch (j) { 53 case 0: 54 switch (k) { 55 case 0: return 0; 56 case 1: return 0; 57 case 2: return -1; 58 } 59 case 1: return 0; 60 case 2: 61 switch (k) { 62 case 0: return 1; 63 case 1: return 0; 64 case 2: return 0; 65 } 66 } 67 case 2: 68 switch (j) { 69 case 0: 70 switch (k) { 71 case 0: return 0; 72 case 1: return 1; 73 case 2: return 0; 74 } 75 case 1: 76 switch (k) { 77 case 0: return -1; 78 case 1: return 0; 79 case 2: return 0; 80 } 81 case 2: return 0; 82 } 83 } 84 return 0; 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "DMPlexProjectRigidBody" 89 static PetscErrorCode DMPlexProjectRigidBody(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx) 90 { 91 PetscInt *ctxInt = (PetscInt *) ctx; 92 PetscInt dim2 = ctxInt[0]; 93 PetscInt d = ctxInt[1]; 94 PetscInt i, j, k = dim > 2 ? d - dim : d; 95 96 PetscFunctionBegin; 97 if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2); 98 for (i = 0; i < dim; i++) mode[i] = 0.; 99 if (d < dim) { 100 mode[d] = 1.; 101 } else { 102 for (i = 0; i < dim; i++) { 103 for (j = 0; j < dim; j++) { 104 mode[j] += epsilon(i, j, k)*X[i]; 105 } 106 } 107 } 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "DMPlexCreateRigidBody" 113 /*@C 114 DMPlexCreateRigidBody - for the default global section, create rigid body modes from coordinates 115 116 Collective on DM 117 118 Input Arguments: 119 . dm - the DM 120 121 Output Argument: 122 . sp - the null space 123 124 Note: This is necessary to take account of Dirichlet conditions on the displacements 125 126 Level: advanced 127 128 .seealso: MatNullSpaceCreate() 129 @*/ 130 PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp) 131 { 132 MPI_Comm comm; 133 Vec mode[6]; 134 PetscSection section, globalSection; 135 PetscInt dim, dimEmbed, n, m, d, i, j; 136 PetscErrorCode ierr; 137 138 PetscFunctionBegin; 139 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 140 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 141 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 142 if (dim == 1) { 143 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 144 PetscFunctionReturn(0); 145 } 146 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 147 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 148 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 149 m = (dim*(dim+1))/2; 150 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 151 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 152 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 153 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 154 for (d = 0; d < m; d++) { 155 PetscInt ctx[2]; 156 void *voidctx = (void *) (&ctx[0]); 157 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody; 158 159 ctx[0] = dimEmbed; 160 ctx[1] = d; 161 ierr = DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 162 } 163 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 164 /* Orthonormalize system */ 165 for (i = dim; i < m; ++i) { 166 PetscScalar dots[6]; 167 168 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 169 for (j = 0; j < i; ++j) dots[j] *= -1.0; 170 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 171 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 172 } 173 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 174 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "DMPlexSetMaxProjectionHeight" 180 /*@ 181 DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs 182 are computed by associating the basis function with one of the mesh points in its transitively-closed support, and 183 evaluating the dual space basis of that point. A basis function is associated with the point in its 184 transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum 185 projection height, which is set with this function. By default, the maximum projection height is zero, which means 186 that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior 187 basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face. 188 189 Input Parameters: 190 + dm - the DMPlex object 191 - height - the maximum projection height >= 0 192 193 Level: advanced 194 195 .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 196 @*/ 197 PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height) 198 { 199 DM_Plex *plex = (DM_Plex *) dm->data; 200 201 PetscFunctionBegin; 202 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 203 plex->maxProjectionHeight = height; 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "DMPlexGetMaxProjectionHeight" 209 /*@ 210 DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in 211 DMPlexProjectXXXLocal() functions. 212 213 Input Parameters: 214 . dm - the DMPlex object 215 216 Output Parameters: 217 . height - the maximum projection height 218 219 Level: intermediate 220 221 .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 222 @*/ 223 PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height) 224 { 225 DM_Plex *plex = (DM_Plex *) dm->data; 226 227 PetscFunctionBegin; 228 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 229 *height = plex->maxProjectionHeight; 230 PetscFunctionReturn(0); 231 } 232 233 #undef __FUNCT__ 234 #define __FUNCT__ "DMProjectFunctionLabelLocal_Plex" 235 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) 236 { 237 PetscDualSpace *sp, *cellsp; 238 PetscInt *numComp; 239 PetscSection section; 240 PetscScalar *values; 241 PetscBool *fieldActive; 242 PetscInt numFields, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, cStart, cEnd, cEndInterior, f, d, v, i, comp, maxHeight, h; 243 PetscErrorCode ierr; 244 245 PetscFunctionBegin; 246 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 247 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 248 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 249 if (cEnd <= cStart) PetscFunctionReturn(0); 250 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 251 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 252 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 253 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 254 ierr = PetscMalloc2(numFields,&sp,numFields,&numComp);CHKERRQ(ierr); 255 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 256 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);} 257 if (maxHeight > 0) {ierr = PetscMalloc1(numFields,&cellsp);CHKERRQ(ierr);} 258 else {cellsp = sp;} 259 for (h = 0; h <= maxHeight; h++) { 260 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 261 if (!h) {pStart = cStart; pEnd = cEnd;} 262 if (pEnd <= pStart) continue; 263 totDim = 0; 264 for (f = 0; f < numFields; ++f) { 265 PetscObject obj; 266 PetscClassId id; 267 268 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 269 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 270 if (id == PETSCFE_CLASSID) { 271 PetscFE fe = (PetscFE) obj; 272 273 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 274 if (!h) { 275 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 276 sp[f] = cellsp[f]; 277 } else { 278 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr); 279 if (!sp[f]) continue; 280 } 281 } else if (id == PETSCFV_CLASSID) { 282 PetscFV fv = (PetscFV) obj; 283 284 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 285 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 286 sp[f] = cellsp[f]; 287 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 288 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 289 totDim += spDim*numComp[f]; 290 } 291 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 292 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim); 293 if (!totDim) continue; 294 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 295 ierr = DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 296 for (f = 0; f < numFields; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 297 for (i = 0; i < numIds; ++i) { 298 IS pointIS; 299 const PetscInt *points; 300 PetscInt n, p; 301 302 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 303 if (!pointIS) continue; /* No points with that id on this process */ 304 ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 305 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 306 for (p = 0; p < n; ++p) { 307 const PetscInt point = points[p]; 308 PetscFECellGeom geom; 309 310 if ((point < pStart) || (point >= pEnd)) continue; 311 ierr = DMPlexComputeCellGeometryFEM(dm, point, NULL, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr); 312 geom.dim = dim - h; 313 geom.dimEmbed = dimEmbed; 314 for (f = 0, v = 0; f < numFields; ++f) { 315 void * const ctx = ctxs ? ctxs[f] : NULL; 316 317 if (!sp[f]) continue; 318 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 319 for (d = 0; d < spDim; ++d) { 320 if (funcs[f]) { 321 ierr = PetscDualSpaceApply(sp[f], d, time, &geom, numComp[f], funcs[f], ctx, &values[v]); 322 if (ierr) { 323 PetscErrorCode ierr2; 324 ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2); 325 ierr2 = DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr2); 326 CHKERRQ(ierr); 327 } 328 } else { 329 for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0; 330 } 331 v += numComp[f]; 332 } 333 } 334 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);CHKERRQ(ierr); 335 } 336 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 337 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 338 } 339 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 340 ierr = DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 341 } 342 ierr = PetscFree2(sp, numComp);CHKERRQ(ierr); 343 if (maxHeight > 0) { 344 ierr = PetscFree(cellsp);CHKERRQ(ierr); 345 } 346 PetscFunctionReturn(0); 347 } 348 349 #undef __FUNCT__ 350 #define __FUNCT__ "DMProjectFunctionLocal_Plex" 351 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 352 { 353 PetscDualSpace *sp, *cellsp; 354 PetscInt *numComp; 355 PetscSection section; 356 PetscScalar *values; 357 PetscInt Nf, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, p, cStart, cEnd, cEndInterior, f, d, v, comp, h, maxHeight; 358 PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE; 359 PetscErrorCode ierr; 360 361 PetscFunctionBegin; 362 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 363 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 364 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 365 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 366 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 367 ierr = PetscMalloc3(Nf, &isFE, Nf, &sp, Nf, &numComp);CHKERRQ(ierr); 368 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 369 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 370 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 371 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);} 372 if (maxHeight > 0) { 373 ierr = PetscMalloc1(Nf,&cellsp);CHKERRQ(ierr); 374 } 375 else { 376 cellsp = sp; 377 } 378 for (f = 0; f < Nf; ++f) { 379 PetscObject obj; 380 PetscClassId id; 381 382 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 383 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 384 if (id == PETSCFE_CLASSID) { 385 PetscFE fe = (PetscFE) obj; 386 387 hasFE = PETSC_TRUE; 388 isFE[f] = PETSC_TRUE; 389 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 390 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 391 } else if (id == PETSCFV_CLASSID) { 392 PetscFV fv = (PetscFV) obj; 393 394 hasFV = PETSC_TRUE; 395 isFE[f] = PETSC_FALSE; 396 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 397 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 398 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 399 } 400 for (h = 0; h <= maxHeight; h++) { 401 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 402 if (!h) {pStart = cStart; pEnd = cEnd;} 403 if (pEnd <= pStart) continue; 404 totDim = 0; 405 for (f = 0; f < Nf; ++f) { 406 if (!h) { 407 sp[f] = cellsp[f]; 408 } 409 else { 410 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr); 411 if (!sp[f]) { 412 continue; 413 } 414 } 415 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 416 totDim += spDim*numComp[f]; 417 } 418 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 419 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim); 420 if (!totDim) continue; 421 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 422 for (p = pStart; p < pEnd; ++p) { 423 PetscFECellGeom fegeom; 424 PetscFVCellGeom fvgeom; 425 426 if (hasFE) { 427 ierr = DMPlexComputeCellGeometryFEM(dm, p, NULL, fegeom.v0, fegeom.J, NULL, &fegeom.detJ);CHKERRQ(ierr); 428 fegeom.dim = dim - h; 429 fegeom.dimEmbed = dimEmbed; 430 } 431 if (hasFV) {ierr = DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);CHKERRQ(ierr);} 432 for (f = 0, v = 0; f < Nf; ++f) { 433 void * const ctx = ctxs ? ctxs[f] : NULL; 434 435 if (!sp[f]) continue; 436 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 437 for (d = 0; d < spDim; ++d) { 438 if (funcs[f]) { 439 if (isFE[f]) {ierr = PetscDualSpaceApply(sp[f], d, time, &fegeom, numComp[f], funcs[f], ctx, &values[v]);} 440 else {ierr = PetscDualSpaceApplyFVM(sp[f], d, time, &fvgeom, numComp[f], funcs[f], ctx, &values[v]);} 441 if (ierr) { 442 PetscErrorCode ierr2; 443 ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2); 444 CHKERRQ(ierr); 445 } 446 } else { 447 for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0; 448 } 449 v += numComp[f]; 450 } 451 } 452 ierr = DMPlexVecSetClosure(dm, section, localX, p, values, mode);CHKERRQ(ierr); 453 } 454 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 455 } 456 ierr = PetscFree3(isFE, sp, numComp);CHKERRQ(ierr); 457 if (maxHeight > 0) { 458 ierr = PetscFree(cellsp);CHKERRQ(ierr); 459 } 460 PetscFunctionReturn(0); 461 } 462 463 #undef __FUNCT__ 464 #define __FUNCT__ "DMProjectFieldLocal_Plex" 465 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, Vec localU, 466 void (**funcs)(PetscInt, PetscInt, PetscInt, 467 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 468 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 469 PetscReal, const PetscReal[], PetscScalar[]), 470 InsertMode mode, Vec localX) 471 { 472 DM dmAux; 473 PetscDS prob, probAux = NULL; 474 Vec A; 475 PetscSection section, sectionAux = NULL; 476 PetscDualSpace *sp; 477 PetscFECellGeom cgeom; 478 PetscReal ****B, ****D, ****BAux, ****DAux; 479 PetscScalar *values, *u, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux; 480 PetscReal *x; 481 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 482 PetscInt Nf, NfAux = 0, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight; 483 PetscErrorCode ierr; 484 485 PetscFunctionBegin; 486 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 487 if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");} 488 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 489 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 490 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 491 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 492 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 493 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 494 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 495 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 496 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 497 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 498 ierr = PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);CHKERRQ(ierr); 499 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 500 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 501 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 502 if (dmAux) { 503 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 504 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 505 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 506 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 507 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 508 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 509 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 510 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 511 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 512 } 513 ierr = PetscMalloc1(Nf, &sp);CHKERRQ(ierr); 514 for (f = 0, v = 0; f < Nf; ++f) { 515 PetscObject obj; 516 PetscClassId id; 517 518 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 519 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 520 if (id == PETSCFE_CLASSID) { 521 PetscFE fe = (PetscFE) obj; 522 523 ierr = PetscFEGetDualSpace(fe, &sp[f]);CHKERRQ(ierr); 524 } else if (id == PETSCFV_CLASSID) { 525 PetscFV fv = (PetscFV) obj; 526 527 ierr = PetscFVGetDualSpace(fv, &sp[f]);CHKERRQ(ierr); 528 } 529 } 530 ierr = PetscMalloc4(Nf, &B, Nf, &D, Nf, &BAux, Nf, &DAux);CHKERRQ(ierr); 531 for (f = 0, v = 0; f < Nf; ++f) { 532 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 533 ierr = PetscMalloc4(spDim, &B[f], spDim, &D[f], spDim, &BAux[f], spDim, &DAux[f]);CHKERRQ(ierr); 534 for (d = 0; d < spDim; ++d) { 535 PetscQuadrature quad; 536 const PetscReal *points; 537 PetscInt numPoints, f2; 538 539 ierr = PetscMalloc4(Nf, &B[f][d], Nf, &D[f][d], NfAux, &BAux[f][d], NfAux, &DAux[f][d]);CHKERRQ(ierr); 540 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 541 ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, NULL);CHKERRQ(ierr); 542 for (f2 = 0; f2 < Nf; ++f2) { 543 PetscObject obj; 544 PetscClassId id; 545 546 ierr = PetscDSGetDiscretization(prob, f2, &obj);CHKERRQ(ierr); 547 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 548 if (id == PETSCFE_CLASSID) { 549 PetscFE fe = (PetscFE) obj; 550 551 ierr = PetscFEGetTabulation(fe, numPoints, points, &B[f][d][f2], &D[f][d][f2], NULL);CHKERRQ(ierr); 552 } else if (id == PETSCFV_CLASSID) { 553 PetscFV fv = (PetscFV) obj; 554 555 ierr = PetscFVGetTabulation(fv, numPoints, points, &B[f][d][f2], &D[f][d][f2], NULL);CHKERRQ(ierr); 556 } 557 } 558 for (f2 = 0; f2 < NfAux; ++f2) { 559 PetscObject obj; 560 PetscClassId id; 561 562 ierr = PetscDSGetDiscretization(probAux, f2, &obj);CHKERRQ(ierr); 563 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 564 if (id == PETSCFE_CLASSID) { 565 PetscFE fe = (PetscFE) obj; 566 567 ierr = PetscFEGetTabulation(fe, numPoints, points, &BAux[f][d][f2], &DAux[f][d][f2], NULL);CHKERRQ(ierr); 568 } else if (id == PETSCFV_CLASSID) { 569 PetscFV fv = (PetscFV) obj; 570 571 ierr = PetscFVGetTabulation(fv, numPoints, points, &BAux[f][d][f2], &DAux[f][d][f2], NULL);CHKERRQ(ierr); 572 } 573 } 574 } 575 } 576 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 577 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 578 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 579 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 580 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 581 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 582 for (c = cStart; c < cEnd; ++c) { 583 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 584 585 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom.v0, cgeom.J, cgeom.invJ, &cgeom.detJ);CHKERRQ(ierr); 586 ierr = DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 587 if (dmAux) {ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 588 for (f = 0, v = 0; f < Nf; ++f) { 589 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 590 for (d = 0; d < spDim; ++d) { 591 PetscQuadrature quad; 592 const PetscReal *points, *weights; 593 PetscInt numPoints, q; 594 595 if (funcs[f]) { 596 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 597 ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);CHKERRQ(ierr); 598 for (q = 0; q < numPoints; ++q) { 599 CoordinatesRefToReal(dim, dim, cgeom.v0, cgeom.J, &points[q*dim], x); 600 EvaluateFieldJets(dim, Nf, Nb, Nc, q, B[f][d], D[f][d], refSpaceDer, cgeom.invJ, coefficients, NULL, u, u_x, NULL); 601 if (probAux) EvaluateFieldJets(dim, NfAux, NbAux, NcAux, q, BAux[f][d], DAux[f][d], refSpaceDerAux, cgeom.invJ, coefficientsAux, NULL, a, a_x, NULL); 602 (*funcs[f])(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &values[v]); 603 } 604 } else { 605 for (comp = 0; comp < Nc[f]; ++comp) values[v+comp] = 0.0; 606 } 607 v += Nc[f]; 608 } 609 } 610 ierr = DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 611 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 612 ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 613 } 614 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 615 for (f = 0, v = 0; f < Nf; ++f) { 616 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 617 for (d = 0; d < spDim; ++d) { 618 PetscInt f2; 619 620 for (f2 = 0; f2 < Nf; ++f2) { 621 PetscObject obj; 622 PetscClassId id; 623 624 ierr = PetscDSGetDiscretization(prob, f2, &obj);CHKERRQ(ierr); 625 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 626 if (id == PETSCFE_CLASSID) { 627 PetscFE fe = (PetscFE) obj; 628 629 ierr = PetscFERestoreTabulation(fe, 0, NULL, &B[f][d][f2], &D[f][d][f2], NULL);CHKERRQ(ierr); 630 } else if (id == PETSCFV_CLASSID) { 631 PetscFV fv = (PetscFV) obj; 632 633 ierr = PetscFVGetTabulation(fv, 0, NULL, &B[f][d][f2], &D[f][d][f2], NULL);CHKERRQ(ierr); 634 } 635 } 636 for (f2 = 0; f2 < NfAux; ++f2) { 637 PetscObject obj; 638 PetscClassId id; 639 640 ierr = PetscDSGetDiscretization(probAux, f2, &obj);CHKERRQ(ierr); 641 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 642 if (id == PETSCFE_CLASSID) { 643 PetscFE fe = (PetscFE) obj; 644 645 ierr = PetscFERestoreTabulation(fe, 0, NULL, &BAux[f][d][f2], &DAux[f][d][f2], NULL);CHKERRQ(ierr); 646 } else if (id == PETSCFV_CLASSID) { 647 PetscFV fv = (PetscFV) obj; 648 649 ierr = PetscFVGetTabulation(fv, 0, NULL, &BAux[f][d][f2], &DAux[f][d][f2], NULL);CHKERRQ(ierr); 650 } 651 } 652 ierr = PetscFree4(B[f][d], D[f][d], BAux[f][d], DAux[f][d]);CHKERRQ(ierr); 653 } 654 ierr = PetscFree4(B[f], D[f], BAux[f], DAux[f]);CHKERRQ(ierr); 655 } 656 ierr = PetscFree4(B, D, BAux, DAux);CHKERRQ(ierr); 657 ierr = PetscFree(sp);CHKERRQ(ierr); 658 PetscFunctionReturn(0); 659 } 660 661 #undef __FUNCT__ 662 #define __FUNCT__ "DMPlexInsertBoundaryValues_FEM_Internal" 663 static PetscErrorCode DMPlexInsertBoundaryValues_FEM_Internal(DM dm, PetscReal time, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX) 664 { 665 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 666 void **ctxs; 667 PetscInt numFields; 668 PetscErrorCode ierr; 669 670 PetscFunctionBegin; 671 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 672 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 673 funcs[field] = func; 674 ctxs[field] = ctx; 675 ierr = DMProjectFunctionLabelLocal(dm, time, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 676 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 677 PetscFunctionReturn(0); 678 } 679 680 #undef __FUNCT__ 681 #define __FUNCT__ "DMPlexInsertBoundaryValues_FVM_Internal" 682 /* This ignores numcomps/comps */ 683 static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, 684 PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 685 { 686 PetscDS prob; 687 PetscSF sf; 688 DM dmFace, dmCell, dmGrad; 689 const PetscScalar *facegeom, *cellgeom = NULL, *grad; 690 const PetscInt *leaves; 691 PetscScalar *x, *fx; 692 PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i; 693 PetscErrorCode ierr, ierru = 0; 694 695 PetscFunctionBegin; 696 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 697 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 698 nleaves = PetscMax(0, nleaves); 699 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 700 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 701 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 702 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 703 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 704 if (cellGeometry) { 705 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 706 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 707 } 708 if (Grad) { 709 PetscFV fv; 710 711 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr); 712 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 713 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 714 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 715 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 716 } 717 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 718 for (i = 0; i < numids; ++i) { 719 IS faceIS; 720 const PetscInt *faces; 721 PetscInt numFaces, f; 722 723 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 724 if (!faceIS) continue; /* No points with that id on this process */ 725 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 726 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 727 for (f = 0; f < numFaces; ++f) { 728 const PetscInt face = faces[f], *cells; 729 PetscFVFaceGeom *fg; 730 731 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 732 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 733 if (loc >= 0) continue; 734 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 735 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 736 if (Grad) { 737 PetscFVCellGeom *cg; 738 PetscScalar *cx, *cgrad; 739 PetscScalar *xG; 740 PetscReal dx[3]; 741 PetscInt d; 742 743 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 744 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 745 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 746 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 747 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 748 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 749 ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx); 750 if (ierru) { 751 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 752 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 753 goto cleanup; 754 } 755 } else { 756 PetscScalar *xI; 757 PetscScalar *xG; 758 759 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 760 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 761 ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx); 762 if (ierru) { 763 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 764 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 765 goto cleanup; 766 } 767 } 768 } 769 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 770 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 771 } 772 cleanup: 773 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 774 if (Grad) { 775 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 776 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 777 } 778 if (cellGeometry) {ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);} 779 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 780 CHKERRQ(ierru); 781 PetscFunctionReturn(0); 782 } 783 784 #undef __FUNCT__ 785 #define __FUNCT__ "DMPlexInsertBoundaryValues" 786 /*@ 787 DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector 788 789 Input Parameters: 790 + dm - The DM 791 . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions 792 . time - The time 793 . faceGeomFVM - Face geometry data for FV discretizations 794 . cellGeomFVM - Cell geometry data for FV discretizations 795 - gradFVM - Gradient reconstruction data for FV discretizations 796 797 Output Parameters: 798 . locX - Solution updated with boundary values 799 800 Level: developer 801 802 .seealso: DMProjectFunctionLabelLocal() 803 @*/ 804 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 805 { 806 PetscInt numBd, b; 807 PetscErrorCode ierr; 808 809 PetscFunctionBegin; 810 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 811 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 812 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 813 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 814 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 815 ierr = PetscDSGetNumBoundary(dm->prob, &numBd);CHKERRQ(ierr); 816 for (b = 0; b < numBd; ++b) { 817 PetscBool isEssential; 818 const char *labelname; 819 DMLabel label; 820 PetscInt field; 821 PetscObject obj; 822 PetscClassId id; 823 void (*func)(); 824 PetscInt numids; 825 const PetscInt *ids; 826 void *ctx; 827 828 ierr = DMGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 829 if (insertEssential != isEssential) continue; 830 ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr); 831 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 832 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 833 if (id == PETSCFE_CLASSID) { 834 if (!isEssential) continue; /* for FEM, there is no insertion to be done for non-essential boundary conditions */ 835 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 836 ierr = DMPlexInsertBoundaryValues_FEM_Internal(dm, time, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 837 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 838 } else if (id == PETSCFV_CLASSID) { 839 if (!faceGeomFVM) continue; 840 ierr = DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, 841 field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 842 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 843 } 844 PetscFunctionReturn(0); 845 } 846 847 #undef __FUNCT__ 848 #define __FUNCT__ "DMComputeL2Diff_Plex" 849 PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 850 { 851 const PetscInt debug = 0; 852 PetscSection section; 853 PetscQuadrature quad; 854 Vec localX; 855 PetscScalar *funcVal, *interpolant; 856 PetscReal *coords, *v0, *J, *invJ, detJ; 857 PetscReal localDiff = 0.0; 858 const PetscReal *quadPoints, *quadWeights; 859 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 860 PetscErrorCode ierr; 861 862 PetscFunctionBegin; 863 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 864 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 865 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 866 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 867 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 868 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 869 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 870 for (field = 0; field < numFields; ++field) { 871 PetscObject obj; 872 PetscClassId id; 873 PetscInt Nc; 874 875 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 876 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 877 if (id == PETSCFE_CLASSID) { 878 PetscFE fe = (PetscFE) obj; 879 880 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 881 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 882 } else if (id == PETSCFV_CLASSID) { 883 PetscFV fv = (PetscFV) obj; 884 885 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 886 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 887 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 888 numComponents += Nc; 889 } 890 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 891 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 892 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 893 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 894 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 895 for (c = cStart; c < cEnd; ++c) { 896 PetscScalar *x = NULL; 897 PetscReal elemDiff = 0.0; 898 899 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 900 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 901 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 902 903 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 904 PetscObject obj; 905 PetscClassId id; 906 void * const ctx = ctxs ? ctxs[field] : NULL; 907 PetscInt Nb, Nc, q, fc; 908 909 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 910 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 911 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 912 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 913 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 914 if (debug) { 915 char title[1024]; 916 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 917 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 918 } 919 for (q = 0; q < numQuadPoints; ++q) { 920 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 921 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx); 922 if (ierr) { 923 PetscErrorCode ierr2; 924 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 925 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 926 ierr2 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2); 927 CHKERRQ(ierr); 928 } 929 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 930 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 931 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 932 for (fc = 0; fc < Nc; ++fc) { 933 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 934 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 935 } 936 } 937 fieldOffset += Nb*Nc; 938 } 939 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 940 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 941 localDiff += elemDiff; 942 } 943 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 944 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 945 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 946 *diff = PetscSqrtReal(*diff); 947 PetscFunctionReturn(0); 948 } 949 950 #undef __FUNCT__ 951 #define __FUNCT__ "DMComputeL2GradientDiff_Plex" 952 PetscErrorCode DMComputeL2GradientDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 953 { 954 const PetscInt debug = 0; 955 PetscSection section; 956 PetscQuadrature quad; 957 Vec localX; 958 PetscScalar *funcVal, *interpolantVec; 959 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 960 PetscReal localDiff = 0.0; 961 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp; 962 PetscErrorCode ierr; 963 964 PetscFunctionBegin; 965 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 966 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 967 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 968 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 969 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 970 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 971 for (field = 0; field < numFields; ++field) { 972 PetscFE fe; 973 PetscInt Nc; 974 975 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 976 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 977 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 978 numComponents += Nc; 979 } 980 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 981 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 982 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 983 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 984 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 985 for (c = cStart; c < cEnd; ++c) { 986 PetscScalar *x = NULL; 987 PetscReal elemDiff = 0.0; 988 989 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 990 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 991 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 992 993 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 994 PetscFE fe; 995 void * const ctx = ctxs ? ctxs[field] : NULL; 996 const PetscReal *quadPoints, *quadWeights; 997 PetscReal *basisDer; 998 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 999 1000 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 1001 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1002 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1003 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 1004 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 1005 if (debug) { 1006 char title[1024]; 1007 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1008 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 1009 } 1010 for (q = 0; q < numQuadPoints; ++q) { 1011 for (d = 0; d < dim; d++) { 1012 coords[d] = v0[d]; 1013 for (e = 0; e < dim; e++) { 1014 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 1015 } 1016 } 1017 ierr = (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx); 1018 if (ierr) { 1019 PetscErrorCode ierr2; 1020 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 1021 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 1022 ierr2 = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr2); 1023 CHKERRQ(ierr); 1024 } 1025 for (fc = 0; fc < Ncomp; ++fc) { 1026 PetscScalar interpolant = 0.0; 1027 1028 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 1029 for (f = 0; f < Nb; ++f) { 1030 const PetscInt fidx = f*Ncomp+fc; 1031 1032 for (d = 0; d < dim; ++d) { 1033 realSpaceDer[d] = 0.0; 1034 for (g = 0; g < dim; ++g) { 1035 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 1036 } 1037 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 1038 } 1039 } 1040 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 1041 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 1042 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 1043 } 1044 } 1045 comp += Ncomp; 1046 fieldOffset += Nb*Ncomp; 1047 } 1048 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1049 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 1050 localDiff += elemDiff; 1051 } 1052 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 1053 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1054 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1055 *diff = PetscSqrtReal(*diff); 1056 PetscFunctionReturn(0); 1057 } 1058 1059 #undef __FUNCT__ 1060 #define __FUNCT__ "DMComputeL2FieldDiff_Plex" 1061 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 1062 { 1063 const PetscInt debug = 0; 1064 PetscSection section; 1065 PetscQuadrature quad; 1066 Vec localX; 1067 PetscScalar *funcVal, *interpolant; 1068 PetscReal *coords, *v0, *J, *invJ, detJ; 1069 PetscReal *localDiff; 1070 const PetscReal *quadPoints, *quadWeights; 1071 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1072 PetscErrorCode ierr; 1073 1074 PetscFunctionBegin; 1075 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1076 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1077 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1078 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1079 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1080 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1081 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1082 for (field = 0; field < numFields; ++field) { 1083 PetscObject obj; 1084 PetscClassId id; 1085 PetscInt Nc; 1086 1087 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1088 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1089 if (id == PETSCFE_CLASSID) { 1090 PetscFE fe = (PetscFE) obj; 1091 1092 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1093 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1094 } else if (id == PETSCFV_CLASSID) { 1095 PetscFV fv = (PetscFV) obj; 1096 1097 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1098 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1099 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1100 numComponents += Nc; 1101 } 1102 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1103 ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1104 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1105 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1106 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1107 for (c = cStart; c < cEnd; ++c) { 1108 PetscScalar *x = NULL; 1109 1110 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1111 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1112 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1113 1114 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1115 PetscObject obj; 1116 PetscClassId id; 1117 void * const ctx = ctxs ? ctxs[field] : NULL; 1118 PetscInt Nb, Nc, q, fc; 1119 1120 PetscReal elemDiff = 0.0; 1121 1122 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1123 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1124 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1125 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1126 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1127 if (debug) { 1128 char title[1024]; 1129 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1130 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 1131 } 1132 for (q = 0; q < numQuadPoints; ++q) { 1133 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 1134 ierr = (*funcs[field])(dim, time, coords, numFields, funcVal, ctx); 1135 if (ierr) { 1136 PetscErrorCode ierr2; 1137 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 1138 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 1139 ierr2 = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2); 1140 CHKERRQ(ierr); 1141 } 1142 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1143 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1144 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1145 for (fc = 0; fc < Nc; ++fc) { 1146 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d point %g %g %g diff %g\n", c, field, dim > 0 ? coords[0] : 0., dim > 1 ? coords[1] : 0., dim > 2 ? coords[2] : 0., PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 1147 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1148 } 1149 } 1150 fieldOffset += Nb*Nc; 1151 localDiff[field] += elemDiff; 1152 } 1153 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1154 } 1155 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1156 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1157 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 1158 ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1159 PetscFunctionReturn(0); 1160 } 1161 1162 #undef __FUNCT__ 1163 #define __FUNCT__ "DMPlexComputeL2DiffVec" 1164 /*@C 1165 DMPlexComputeL2DiffVec - This function computes the cellwise L_2 difference between a function u and an FEM interpolant solution u_h, and stores it in a Vec. 1166 1167 Input Parameters: 1168 + dm - The DM 1169 . time - The time 1170 . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation 1171 . ctxs - Optional array of contexts to pass to each function, or NULL. 1172 - X - The coefficient vector u_h 1173 1174 Output Parameter: 1175 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 1176 1177 Level: developer 1178 1179 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 1180 @*/ 1181 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 1182 { 1183 PetscSection section; 1184 PetscQuadrature quad; 1185 Vec localX; 1186 PetscScalar *funcVal, *interpolant; 1187 PetscReal *coords, *v0, *J, *invJ, detJ; 1188 const PetscReal *quadPoints, *quadWeights; 1189 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1190 PetscErrorCode ierr; 1191 1192 PetscFunctionBegin; 1193 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 1194 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1195 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1196 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1197 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1198 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1199 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1200 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1201 for (field = 0; field < numFields; ++field) { 1202 PetscObject obj; 1203 PetscClassId id; 1204 PetscInt Nc; 1205 1206 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1207 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1208 if (id == PETSCFE_CLASSID) { 1209 PetscFE fe = (PetscFE) obj; 1210 1211 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1212 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1213 } else if (id == PETSCFV_CLASSID) { 1214 PetscFV fv = (PetscFV) obj; 1215 1216 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1217 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1218 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1219 numComponents += Nc; 1220 } 1221 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1222 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1223 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1224 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1225 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1226 for (c = cStart; c < cEnd; ++c) { 1227 PetscScalar *x = NULL; 1228 PetscScalar elemDiff = 0.0; 1229 1230 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1231 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1232 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1233 1234 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1235 PetscObject obj; 1236 PetscClassId id; 1237 void * const ctx = ctxs ? ctxs[field] : NULL; 1238 PetscInt Nb, Nc, q, fc; 1239 1240 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1241 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1242 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1243 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1244 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1245 if (funcs[field]) { 1246 for (q = 0; q < numQuadPoints; ++q) { 1247 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 1248 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx); 1249 if (ierr) { 1250 PetscErrorCode ierr2; 1251 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 1252 ierr2 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2); 1253 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 1254 CHKERRQ(ierr); 1255 } 1256 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1257 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1258 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1259 for (fc = 0; fc < Nc; ++fc) { 1260 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1261 } 1262 } 1263 } 1264 fieldOffset += Nb*Nc; 1265 } 1266 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1267 ierr = VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);CHKERRQ(ierr); 1268 } 1269 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1270 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1271 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 1272 PetscFunctionReturn(0); 1273 } 1274 1275 #undef __FUNCT__ 1276 #define __FUNCT__ "DMPlexComputeIntegralFEM" 1277 /*@ 1278 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 1279 1280 Input Parameters: 1281 + dm - The mesh 1282 . X - Local input vector 1283 - user - The user context 1284 1285 Output Parameter: 1286 . integral - Local integral for each field 1287 1288 Level: developer 1289 1290 .seealso: DMPlexComputeResidualFEM() 1291 @*/ 1292 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 1293 { 1294 DM_Plex *mesh = (DM_Plex *) dm->data; 1295 DM dmAux, dmGrad; 1296 Vec localX, A, cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL; 1297 PetscDS prob, probAux = NULL; 1298 PetscSection section, sectionAux; 1299 PetscFV fvm = NULL; 1300 PetscFECellGeom *cgeomFEM; 1301 PetscFVFaceGeom *fgeomFVM; 1302 PetscFVCellGeom *cgeomFVM; 1303 PetscScalar *u, *a = NULL; 1304 const PetscScalar *lgrad; 1305 PetscReal *lintegral; 1306 PetscInt *uOff, *uOff_x, *aOff = NULL; 1307 PetscInt dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, cEndInterior, c; 1308 PetscInt totDim, totDimAux; 1309 PetscBool useFVM = PETSC_FALSE; 1310 PetscErrorCode ierr; 1311 1312 PetscFunctionBegin; 1313 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1314 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1315 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1316 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1317 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1318 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1319 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1320 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1321 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1322 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 1323 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 1324 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1325 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1326 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1327 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1328 numCells = cEnd - cStart; 1329 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1330 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1331 if (dmAux) { 1332 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1333 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1334 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1335 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1336 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1337 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, NULL);CHKERRQ(ierr); 1338 } 1339 for (f = 0; f < Nf; ++f) { 1340 PetscObject obj; 1341 PetscClassId id; 1342 1343 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1344 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1345 if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;} 1346 } 1347 if (useFVM) { 1348 Vec grad; 1349 PetscInt fStart, fEnd; 1350 PetscBool compGrad; 1351 1352 ierr = PetscFVGetComputeGradients(fvm, &compGrad);CHKERRQ(ierr); 1353 ierr = PetscFVSetComputeGradients(fvm, PETSC_TRUE);CHKERRQ(ierr); 1354 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 1355 ierr = DMPlexComputeGradientFVM(dm, fvm, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 1356 ierr = PetscFVSetComputeGradients(fvm, compGrad);CHKERRQ(ierr); 1357 ierr = VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1358 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1359 /* Reconstruct and limit cell gradients */ 1360 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1361 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1362 ierr = DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, localX, grad);CHKERRQ(ierr); 1363 /* Communicate gradient values */ 1364 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1365 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1366 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1367 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1368 /* Handle non-essential (e.g. outflow) boundary values */ 1369 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, localX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 1370 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1371 } 1372 ierr = PetscMalloc3(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeomFEM);CHKERRQ(ierr); 1373 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1374 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 1375 for (c = cStart; c < cEnd; ++c) { 1376 PetscScalar *x = NULL; 1377 PetscInt i; 1378 1379 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);CHKERRQ(ierr); 1380 if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c); 1381 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1382 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1383 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1384 if (dmAux) { 1385 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1386 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1387 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1388 } 1389 } 1390 for (f = 0; f < Nf; ++f) { 1391 PetscObject obj; 1392 PetscClassId id; 1393 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1394 1395 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1396 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1397 if (id == PETSCFE_CLASSID) { 1398 PetscFE fe = (PetscFE) obj; 1399 PetscQuadrature q; 1400 PetscInt Nq, Nb; 1401 1402 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1403 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1404 ierr = PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1405 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1406 blockSize = Nb*Nq; 1407 batchSize = numBlocks * blockSize; 1408 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1409 numChunks = numCells / (numBatches*batchSize); 1410 Ne = numChunks*numBatches*batchSize; 1411 Nr = numCells % (numBatches*batchSize); 1412 offset = numCells - Nr; 1413 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, lintegral);CHKERRQ(ierr); 1414 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 1415 } else if (id == PETSCFV_CLASSID) { 1416 /* PetscFV fv = (PetscFV) obj; */ 1417 PetscInt foff; 1418 PetscPointFunc obj_func; 1419 PetscScalar lint; 1420 1421 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1422 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1423 if (obj_func) { 1424 for (c = 0; c < numCells; ++c) { 1425 PetscScalar *u_x; 1426 1427 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1428 obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, a, NULL, NULL, 0.0, cgeomFVM[c].centroid, &lint); 1429 lintegral[f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1430 } 1431 } 1432 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1433 } 1434 if (useFVM) { 1435 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1436 ierr = VecRestoreArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1437 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1438 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1439 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1440 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1441 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1442 } 1443 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1444 if (mesh->printFEM) { 1445 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1446 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1447 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1448 } 1449 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1450 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1451 ierr = PetscFree3(lintegral,u,cgeomFEM);CHKERRQ(ierr); 1452 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1453 PetscFunctionReturn(0); 1454 } 1455 1456 #undef __FUNCT__ 1457 #define __FUNCT__ "DMPlexComputeInterpolatorNested" 1458 /*@ 1459 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1460 1461 Input Parameters: 1462 + dmf - The fine mesh 1463 . dmc - The coarse mesh 1464 - user - The user context 1465 1466 Output Parameter: 1467 . In - The interpolation matrix 1468 1469 Level: developer 1470 1471 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1472 @*/ 1473 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1474 { 1475 DM_Plex *mesh = (DM_Plex *) dmc->data; 1476 const char *name = "Interpolator"; 1477 PetscDS prob; 1478 PetscFE *feRef; 1479 PetscFV *fvRef; 1480 PetscSection fsection, fglobalSection; 1481 PetscSection csection, cglobalSection; 1482 PetscScalar *elemMat; 1483 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1484 PetscInt cTotDim, rTotDim = 0; 1485 PetscErrorCode ierr; 1486 1487 PetscFunctionBegin; 1488 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1489 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1490 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1491 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1492 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1493 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1494 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1495 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1496 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1497 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1498 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1499 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1500 for (f = 0; f < Nf; ++f) { 1501 PetscObject obj; 1502 PetscClassId id; 1503 PetscInt rNb = 0, Nc = 0; 1504 1505 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1506 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1507 if (id == PETSCFE_CLASSID) { 1508 PetscFE fe = (PetscFE) obj; 1509 1510 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1511 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1512 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1513 } else if (id == PETSCFV_CLASSID) { 1514 PetscFV fv = (PetscFV) obj; 1515 PetscDualSpace Q; 1516 1517 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1518 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1519 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1520 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1521 } 1522 rTotDim += rNb*Nc; 1523 } 1524 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1525 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1526 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1527 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1528 PetscDualSpace Qref; 1529 PetscQuadrature f; 1530 const PetscReal *qpoints, *qweights; 1531 PetscReal *points; 1532 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1533 1534 /* Compose points from all dual basis functionals */ 1535 if (feRef[fieldI]) { 1536 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1537 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1538 } else { 1539 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1540 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1541 } 1542 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1543 for (i = 0; i < fpdim; ++i) { 1544 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1545 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1546 npoints += Np; 1547 } 1548 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1549 for (i = 0, k = 0; i < fpdim; ++i) { 1550 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1551 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1552 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1553 } 1554 1555 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1556 PetscObject obj; 1557 PetscClassId id; 1558 PetscReal *B; 1559 PetscInt NcJ = 0, cpdim = 0, j; 1560 1561 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1562 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1563 if (id == PETSCFE_CLASSID) { 1564 PetscFE fe = (PetscFE) obj; 1565 1566 /* Evaluate basis at points */ 1567 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1568 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1569 /* For now, fields only interpolate themselves */ 1570 if (fieldI == fieldJ) { 1571 if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ); 1572 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1573 for (i = 0, k = 0; i < fpdim; ++i) { 1574 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1575 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1576 for (p = 0; p < Np; ++p, ++k) { 1577 for (j = 0; j < cpdim; ++j) { 1578 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p]; 1579 } 1580 } 1581 } 1582 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1583 } 1584 } else if (id == PETSCFV_CLASSID) { 1585 PetscFV fv = (PetscFV) obj; 1586 1587 /* Evaluate constant function at points */ 1588 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1589 cpdim = 1; 1590 /* For now, fields only interpolate themselves */ 1591 if (fieldI == fieldJ) { 1592 if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ); 1593 for (i = 0, k = 0; i < fpdim; ++i) { 1594 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1595 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1596 for (p = 0; p < Np; ++p, ++k) { 1597 for (j = 0; j < cpdim; ++j) { 1598 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p]; 1599 } 1600 } 1601 } 1602 } 1603 } 1604 offsetJ += cpdim*NcJ; 1605 } 1606 offsetI += fpdim*Nc; 1607 ierr = PetscFree(points);CHKERRQ(ierr); 1608 } 1609 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1610 /* Preallocate matrix */ 1611 { 1612 Mat preallocator; 1613 PetscScalar *vals; 1614 PetscInt *cellCIndices, *cellFIndices; 1615 PetscInt locRows, locCols, cell; 1616 1617 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1618 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1619 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1620 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1621 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1622 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1623 for (cell = cStart; cell < cEnd; ++cell) { 1624 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1625 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1626 } 1627 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1628 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1629 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1630 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1631 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1632 } 1633 /* Fill matrix */ 1634 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1635 for (c = cStart; c < cEnd; ++c) { 1636 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1637 } 1638 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1639 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1640 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1641 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1642 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1643 if (mesh->printFEM) { 1644 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1645 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1646 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1647 } 1648 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1649 PetscFunctionReturn(0); 1650 } 1651 1652 #undef __FUNCT__ 1653 #define __FUNCT__ "DMPlexComputeInterpolatorGeneral" 1654 /*@ 1655 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1656 1657 Input Parameters: 1658 + dmf - The fine mesh 1659 . dmc - The coarse mesh 1660 - user - The user context 1661 1662 Output Parameter: 1663 . In - The interpolation matrix 1664 1665 Level: developer 1666 1667 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1668 @*/ 1669 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1670 { 1671 DM_Plex *mesh = (DM_Plex *) dmf->data; 1672 const char *name = "Interpolator"; 1673 PetscDS prob; 1674 PetscSection fsection, csection, globalFSection, globalCSection; 1675 PetscHashJK ht; 1676 PetscLayout rLayout; 1677 PetscInt *dnz, *onz; 1678 PetscInt locRows, rStart, rEnd; 1679 PetscReal *x, *v0, *J, *invJ, detJ; 1680 PetscReal *v0c, *Jc, *invJc, detJc; 1681 PetscScalar *elemMat; 1682 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1683 PetscErrorCode ierr; 1684 1685 PetscFunctionBegin; 1686 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1687 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1688 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1689 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1690 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1691 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1692 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1693 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1694 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1695 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1696 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1697 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1698 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1699 ierr = PetscMalloc1(totDim*totDim, &elemMat);CHKERRQ(ierr); 1700 1701 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1702 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1703 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1704 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1705 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1706 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1707 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1708 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1709 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1710 for (field = 0; field < Nf; ++field) { 1711 PetscObject obj; 1712 PetscClassId id; 1713 PetscDualSpace Q = NULL; 1714 PetscQuadrature f; 1715 const PetscReal *qpoints; 1716 PetscInt Nc, Np, fpdim, i, d; 1717 1718 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1719 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1720 if (id == PETSCFE_CLASSID) { 1721 PetscFE fe = (PetscFE) obj; 1722 1723 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1724 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1725 } else if (id == PETSCFV_CLASSID) { 1726 PetscFV fv = (PetscFV) obj; 1727 1728 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1729 Nc = 1; 1730 } 1731 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1732 /* For each fine grid cell */ 1733 for (cell = cStart; cell < cEnd; ++cell) { 1734 PetscInt *findices, *cindices; 1735 PetscInt numFIndices, numCIndices; 1736 1737 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1738 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1739 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1740 for (i = 0; i < fpdim; ++i) { 1741 Vec pointVec; 1742 PetscScalar *pV; 1743 PetscSF coarseCellSF = NULL; 1744 const PetscSFNode *coarseCells; 1745 PetscInt numCoarseCells, q, r, c; 1746 1747 /* Get points from the dual basis functional quadrature */ 1748 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1749 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1750 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1751 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1752 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1753 for (q = 0; q < Np; ++q) { 1754 /* Transform point to real space */ 1755 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1756 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1757 } 1758 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1759 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1760 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1761 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1762 /* Update preallocation info */ 1763 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1764 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1765 for (r = 0; r < Nc; ++r) { 1766 PetscHashJKKey key; 1767 PetscHashJKIter missing, iter; 1768 1769 key.j = findices[i*Nc+r]; 1770 if (key.j < 0) continue; 1771 /* Get indices for coarse elements */ 1772 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1773 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1774 for (c = 0; c < numCIndices; ++c) { 1775 key.k = cindices[c]; 1776 if (key.k < 0) continue; 1777 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1778 if (missing) { 1779 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1780 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1781 else ++onz[key.j-rStart]; 1782 } 1783 } 1784 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1785 } 1786 } 1787 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1788 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1789 } 1790 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1791 } 1792 } 1793 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1794 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1795 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1796 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1797 for (field = 0; field < Nf; ++field) { 1798 PetscObject obj; 1799 PetscClassId id; 1800 PetscDualSpace Q = NULL; 1801 PetscQuadrature f; 1802 const PetscReal *qpoints, *qweights; 1803 PetscInt Nc, Np, fpdim, i, d; 1804 1805 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1806 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1807 if (id == PETSCFE_CLASSID) { 1808 PetscFE fe = (PetscFE) obj; 1809 1810 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1811 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1812 } else if (id == PETSCFV_CLASSID) { 1813 PetscFV fv = (PetscFV) obj; 1814 1815 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1816 Nc = 1; 1817 } 1818 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1819 /* For each fine grid cell */ 1820 for (cell = cStart; cell < cEnd; ++cell) { 1821 PetscInt *findices, *cindices; 1822 PetscInt numFIndices, numCIndices; 1823 1824 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1825 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1826 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1827 for (i = 0; i < fpdim; ++i) { 1828 Vec pointVec; 1829 PetscScalar *pV; 1830 PetscSF coarseCellSF = NULL; 1831 const PetscSFNode *coarseCells; 1832 PetscInt numCoarseCells, cpdim, q, c, j; 1833 1834 /* Get points from the dual basis functional quadrature */ 1835 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1836 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1837 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1838 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1839 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1840 for (q = 0; q < Np; ++q) { 1841 /* Transform point to real space */ 1842 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1843 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1844 } 1845 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1846 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1847 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1848 /* Update preallocation info */ 1849 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1850 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1851 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1852 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1853 PetscReal pVReal[3]; 1854 1855 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1856 /* Transform points from real space to coarse reference space */ 1857 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1858 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1859 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1860 1861 if (id == PETSCFE_CLASSID) { 1862 PetscFE fe = (PetscFE) obj; 1863 PetscReal *B; 1864 1865 /* Evaluate coarse basis on contained point */ 1866 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1867 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1868 ierr = PetscMemzero(elemMat, cpdim*Nc*Nc * sizeof(PetscScalar));CHKERRQ(ierr); 1869 /* Get elemMat entries by multiplying by weight */ 1870 for (j = 0; j < cpdim; ++j) { 1871 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell]; 1872 } 1873 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1874 } else { 1875 cpdim = 1; 1876 for (j = 0; j < cpdim; ++j) { 1877 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell]; 1878 } 1879 } 1880 /* Update interpolator */ 1881 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, Nc, numCIndices, elemMat);CHKERRQ(ierr);} 1882 if (numCIndices != cpdim*Nc) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D*%D", numCIndices, cpdim, Nc); 1883 ierr = MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1884 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1885 } 1886 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1887 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1888 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1889 } 1890 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1891 } 1892 } 1893 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1894 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1895 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1896 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1897 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1898 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1899 PetscFunctionReturn(0); 1900 } 1901 1902 #undef __FUNCT__ 1903 #define __FUNCT__ "DMPlexComputeInjectorFEM" 1904 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1905 { 1906 PetscDS prob; 1907 PetscFE *feRef; 1908 PetscFV *fvRef; 1909 Vec fv, cv; 1910 IS fis, cis; 1911 PetscSection fsection, fglobalSection, csection, cglobalSection; 1912 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1913 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1914 PetscErrorCode ierr; 1915 1916 PetscFunctionBegin; 1917 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1918 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1919 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1920 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1921 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1922 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1923 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1924 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1925 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1926 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1927 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1928 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1929 for (f = 0; f < Nf; ++f) { 1930 PetscObject obj; 1931 PetscClassId id; 1932 PetscInt fNb = 0, Nc = 0; 1933 1934 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1935 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1936 if (id == PETSCFE_CLASSID) { 1937 PetscFE fe = (PetscFE) obj; 1938 1939 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1940 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1941 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1942 } else if (id == PETSCFV_CLASSID) { 1943 PetscFV fv = (PetscFV) obj; 1944 PetscDualSpace Q; 1945 1946 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1947 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1948 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1949 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1950 } 1951 fTotDim += fNb*Nc; 1952 } 1953 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1954 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1955 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1956 PetscFE feC; 1957 PetscFV fvC; 1958 PetscDualSpace QF, QC; 1959 PetscInt NcF, NcC, fpdim, cpdim; 1960 1961 if (feRef[field]) { 1962 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1963 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1964 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1965 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1966 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1967 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1968 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1969 } else { 1970 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1971 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1972 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1973 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1974 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1975 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1976 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1977 } 1978 if (NcF != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", NcF, NcC); 1979 for (c = 0; c < cpdim; ++c) { 1980 PetscQuadrature cfunc; 1981 const PetscReal *cqpoints; 1982 PetscInt NpC; 1983 PetscBool found = PETSC_FALSE; 1984 1985 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1986 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1987 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1988 for (f = 0; f < fpdim; ++f) { 1989 PetscQuadrature ffunc; 1990 const PetscReal *fqpoints; 1991 PetscReal sum = 0.0; 1992 PetscInt NpF, comp; 1993 1994 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1995 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1996 if (NpC != NpF) continue; 1997 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1998 if (sum > 1.0e-9) continue; 1999 for (comp = 0; comp < NcC; ++comp) { 2000 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 2001 } 2002 found = PETSC_TRUE; 2003 break; 2004 } 2005 if (!found) { 2006 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 2007 if (fvRef[field]) { 2008 PetscInt comp; 2009 for (comp = 0; comp < NcC; ++comp) { 2010 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 2011 } 2012 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 2013 } 2014 } 2015 offsetC += cpdim*NcC; 2016 offsetF += fpdim*NcF; 2017 } 2018 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 2019 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 2020 2021 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 2022 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 2023 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 2024 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 2025 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 2026 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 2027 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 2028 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 2029 for (c = cStart; c < cEnd; ++c) { 2030 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 2031 for (d = 0; d < cTotDim; ++d) { 2032 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 2033 if ((findices[cellCIndices[d]-startC] >= 0) && (findices[cellCIndices[d]-startC] != cellFIndices[cmap[d]])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coarse dof %d maps to both %d and %d", cindices[cellCIndices[d]-startC], findices[cellCIndices[d]-startC], cellFIndices[cmap[d]]); 2034 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 2035 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 2036 } 2037 } 2038 ierr = PetscFree(cmap);CHKERRQ(ierr); 2039 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 2040 2041 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 2042 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 2043 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 2044 ierr = ISDestroy(&cis);CHKERRQ(ierr); 2045 ierr = ISDestroy(&fis);CHKERRQ(ierr); 2046 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 2047 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 2048 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2049 PetscFunctionReturn(0); 2050 } 2051