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 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 278 } else { 279 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr); 280 if (!sp[f]) continue; 281 } 282 } else if (id == PETSCFV_CLASSID) { 283 PetscFV fv = (PetscFV) obj; 284 285 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 286 ierr = PetscFVGetDualSpace(fv, &sp[f]);CHKERRQ(ierr); 287 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 288 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 289 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 290 totDim += spDim*numComp[f]; 291 } 292 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 293 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim); 294 if (!totDim) continue; 295 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 296 ierr = DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 297 for (f = 0; f < numFields; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 298 for (i = 0; i < numIds; ++i) { 299 IS pointIS; 300 const PetscInt *points; 301 PetscInt n, p; 302 303 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 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 if (h) { 340 for (f = 0; f < numFields; ++f) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 341 } 342 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 343 ierr = DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 344 } 345 for (f = 0; f < numFields; ++f) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 346 ierr = PetscFree2(sp, numComp);CHKERRQ(ierr); 347 if (maxHeight > 0) { 348 ierr = PetscFree(cellsp);CHKERRQ(ierr); 349 } 350 PetscFunctionReturn(0); 351 } 352 353 #undef __FUNCT__ 354 #define __FUNCT__ "DMProjectFunctionLocal_Plex" 355 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 356 { 357 PetscDualSpace *sp, *cellsp; 358 PetscInt *numComp; 359 PetscSection section; 360 PetscScalar *values; 361 PetscInt Nf, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, p, cStart, cEnd, cEndInterior, f, d, v, comp, h, maxHeight; 362 PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE; 363 PetscErrorCode ierr; 364 365 PetscFunctionBegin; 366 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 367 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 368 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 369 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 370 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 371 ierr = PetscMalloc3(Nf, &isFE, Nf, &sp, Nf, &numComp);CHKERRQ(ierr); 372 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 373 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 374 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 375 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);} 376 if (maxHeight > 0) { 377 ierr = PetscMalloc1(Nf,&cellsp);CHKERRQ(ierr); 378 } 379 else { 380 cellsp = sp; 381 } 382 for (f = 0; f < Nf; ++f) { 383 PetscObject obj; 384 PetscClassId id; 385 386 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 387 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 388 if (id == PETSCFE_CLASSID) { 389 PetscFE fe = (PetscFE) obj; 390 391 hasFE = PETSC_TRUE; 392 isFE[f] = PETSC_TRUE; 393 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 394 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 395 } else if (id == PETSCFV_CLASSID) { 396 PetscFV fv = (PetscFV) obj; 397 398 hasFV = PETSC_TRUE; 399 isFE[f] = PETSC_FALSE; 400 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 401 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 402 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 403 } 404 for (h = 0; h <= maxHeight; h++) { 405 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 406 if (!h) {pStart = cStart; pEnd = cEnd;} 407 if (pEnd <= pStart) continue; 408 totDim = 0; 409 for (f = 0; f < Nf; ++f) { 410 if (!h) { 411 sp[f] = cellsp[f]; 412 } 413 else { 414 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr); 415 if (!sp[f]) { 416 continue; 417 } 418 } 419 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 420 totDim += spDim*numComp[f]; 421 } 422 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 423 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim); 424 if (!totDim) continue; 425 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 426 for (p = pStart; p < pEnd; ++p) { 427 PetscFECellGeom fegeom; 428 PetscFVCellGeom fvgeom; 429 430 if (hasFE) { 431 ierr = DMPlexComputeCellGeometryFEM(dm, p, NULL, fegeom.v0, fegeom.J, NULL, &fegeom.detJ);CHKERRQ(ierr); 432 fegeom.dim = dim - h; 433 fegeom.dimEmbed = dimEmbed; 434 } 435 if (hasFV) {ierr = DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);CHKERRQ(ierr);} 436 for (f = 0, v = 0; f < Nf; ++f) { 437 void * const ctx = ctxs ? ctxs[f] : NULL; 438 439 if (!sp[f]) continue; 440 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 441 for (d = 0; d < spDim; ++d) { 442 if (funcs[f]) { 443 if (isFE[f]) {ierr = PetscDualSpaceApply(sp[f], d, time, &fegeom, numComp[f], funcs[f], ctx, &values[v]);} 444 else {ierr = PetscDualSpaceApplyFVM(sp[f], d, time, &fvgeom, numComp[f], funcs[f], ctx, &values[v]);} 445 if (ierr) { 446 PetscErrorCode ierr2; 447 ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2); 448 CHKERRQ(ierr); 449 } 450 } else { 451 for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0; 452 } 453 v += numComp[f]; 454 } 455 } 456 ierr = DMPlexVecSetClosure(dm, section, localX, p, values, mode);CHKERRQ(ierr); 457 } 458 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 459 if (h) { 460 for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 461 } 462 } 463 ierr = PetscFree3(isFE, sp, numComp);CHKERRQ(ierr); 464 if (maxHeight > 0) { 465 ierr = PetscFree(cellsp);CHKERRQ(ierr); 466 } 467 PetscFunctionReturn(0); 468 } 469 470 #undef __FUNCT__ 471 #define __FUNCT__ "DMProjectFieldLocal_Plex" 472 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, Vec localU, 473 void (**funcs)(PetscInt, PetscInt, PetscInt, 474 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 475 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 476 PetscReal, const PetscReal[], PetscScalar[]), 477 InsertMode mode, Vec localX) 478 { 479 DM dmAux; 480 PetscDS prob, probAux = NULL; 481 Vec A; 482 PetscSection section, sectionAux = NULL; 483 PetscDualSpace *sp; 484 PetscInt *Ncf; 485 PetscScalar *values, *u, *u_x, *a, *a_x; 486 PetscReal *x, *v0, *J, *invJ, detJ; 487 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 488 PetscInt Nf, NfAux = 0, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight; 489 PetscErrorCode ierr; 490 491 PetscFunctionBegin; 492 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 493 if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");} 494 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 495 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 496 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 497 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 498 ierr = PetscMalloc2(Nf, &sp, Nf, &Ncf);CHKERRQ(ierr); 499 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 500 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 501 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 502 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 503 ierr = PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);CHKERRQ(ierr); 504 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 505 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 506 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 507 if (dmAux) { 508 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 509 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 510 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 511 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 512 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 513 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 514 } 515 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 516 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 517 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 518 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 519 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 520 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 521 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 522 for (c = cStart; c < cEnd; ++c) { 523 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 524 525 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 526 ierr = DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 527 if (dmAux) {ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 528 for (f = 0, v = 0; f < Nf; ++f) { 529 PetscObject obj; 530 PetscClassId id; 531 532 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 533 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 534 if (id == PETSCFE_CLASSID) { 535 PetscFE fe = (PetscFE) obj; 536 537 ierr = PetscFEGetDualSpace(fe, &sp[f]);CHKERRQ(ierr); 538 ierr = PetscFEGetNumComponents(fe, &Ncf[f]);CHKERRQ(ierr); 539 } else if (id == PETSCFV_CLASSID) { 540 PetscFV fv = (PetscFV) obj; 541 542 ierr = PetscFVGetNumComponents(fv, &Ncf[f]);CHKERRQ(ierr); 543 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 544 } 545 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 546 for (d = 0; d < spDim; ++d) { 547 PetscQuadrature quad; 548 const PetscReal *points, *weights; 549 PetscInt numPoints, q; 550 551 if (funcs[f]) { 552 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 553 ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);CHKERRQ(ierr); 554 for (q = 0; q < numPoints; ++q) { 555 CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x); 556 ierr = EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, coefficients, NULL, u, u_x, NULL);CHKERRQ(ierr); 557 ierr = EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);CHKERRQ(ierr); 558 (*funcs[f])(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &values[v]); 559 } 560 } else { 561 for (comp = 0; comp < Ncf[f]; ++comp) values[v+comp] = 0.0; 562 } 563 v += Ncf[f]; 564 } 565 } 566 ierr = DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 567 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 568 ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 569 } 570 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 571 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 572 ierr = PetscFree2(sp, Ncf);CHKERRQ(ierr); 573 PetscFunctionReturn(0); 574 } 575 576 #undef __FUNCT__ 577 #define __FUNCT__ "DMPlexInsertBoundaryValues_FEM_Internal" 578 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) 579 { 580 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 581 void **ctxs; 582 PetscInt numFields; 583 PetscErrorCode ierr; 584 585 PetscFunctionBegin; 586 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 587 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 588 funcs[field] = func; 589 ctxs[field] = ctx; 590 ierr = DMProjectFunctionLabelLocal(dm, time, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 591 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 592 PetscFunctionReturn(0); 593 } 594 595 #undef __FUNCT__ 596 #define __FUNCT__ "DMPlexInsertBoundaryValues_FVM_Internal" 597 /* This ignores numcomps/comps */ 598 static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, 599 PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 600 { 601 PetscDS prob; 602 PetscSF sf; 603 DM dmFace, dmCell, dmGrad; 604 const PetscScalar *facegeom, *cellgeom = NULL, *grad; 605 const PetscInt *leaves; 606 PetscScalar *x, *fx; 607 PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i; 608 PetscErrorCode ierr, ierru = 0; 609 610 PetscFunctionBegin; 611 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 612 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 613 nleaves = PetscMax(0, nleaves); 614 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 615 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 616 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 617 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 618 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 619 if (cellGeometry) { 620 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 621 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 622 } 623 if (Grad) { 624 PetscFV fv; 625 626 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr); 627 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 628 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 629 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 630 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 631 } 632 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 633 for (i = 0; i < numids; ++i) { 634 IS faceIS; 635 const PetscInt *faces; 636 PetscInt numFaces, f; 637 638 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 639 if (!faceIS) continue; /* No points with that id on this process */ 640 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 641 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 642 for (f = 0; f < numFaces; ++f) { 643 const PetscInt face = faces[f], *cells; 644 PetscFVFaceGeom *fg; 645 646 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 647 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 648 if (loc >= 0) continue; 649 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 650 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 651 if (Grad) { 652 PetscFVCellGeom *cg; 653 PetscScalar *cx, *cgrad; 654 PetscScalar *xG; 655 PetscReal dx[3]; 656 PetscInt d; 657 658 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 659 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 660 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 661 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 662 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 663 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 664 ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx); 665 if (ierru) { 666 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 667 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 668 goto cleanup; 669 } 670 } else { 671 PetscScalar *xI; 672 PetscScalar *xG; 673 674 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 675 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 676 ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx); 677 if (ierru) { 678 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 679 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 680 goto cleanup; 681 } 682 } 683 } 684 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 685 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 686 } 687 cleanup: 688 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 689 if (Grad) { 690 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 691 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 692 } 693 if (cellGeometry) {ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);} 694 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 695 CHKERRQ(ierru); 696 PetscFunctionReturn(0); 697 } 698 699 #undef __FUNCT__ 700 #define __FUNCT__ "DMPlexInsertBoundaryValues" 701 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 702 { 703 PetscInt numBd, b; 704 PetscErrorCode ierr; 705 706 PetscFunctionBegin; 707 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 708 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 709 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 710 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 711 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 712 ierr = DMGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 713 for (b = 0; b < numBd; ++b) { 714 PetscBool isEssential; 715 const char *labelname; 716 DMLabel label; 717 PetscInt field; 718 PetscObject obj; 719 PetscClassId id; 720 void (*func)(); 721 PetscInt numids; 722 const PetscInt *ids; 723 void *ctx; 724 725 ierr = DMGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 726 if (insertEssential != isEssential) continue; 727 ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr); 728 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 729 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 730 if (id == PETSCFE_CLASSID) { 731 if (!isEssential) continue; /* for FEM, there is no insertion to be done for non-essential boundary conditions */ 732 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 733 ierr = DMPlexInsertBoundaryValues_FEM_Internal(dm, time, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 734 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 735 } else if (id == PETSCFV_CLASSID) { 736 if (!faceGeomFVM) continue; 737 ierr = DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, 738 field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 739 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 740 } 741 PetscFunctionReturn(0); 742 } 743 744 #undef __FUNCT__ 745 #define __FUNCT__ "DMComputeL2Diff_Plex" 746 PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 747 { 748 const PetscInt debug = 0; 749 PetscSection section; 750 PetscQuadrature quad; 751 Vec localX; 752 PetscScalar *funcVal, *interpolant; 753 PetscReal *coords, *v0, *J, *invJ, detJ; 754 PetscReal localDiff = 0.0; 755 const PetscReal *quadPoints, *quadWeights; 756 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 757 PetscErrorCode ierr; 758 759 PetscFunctionBegin; 760 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 761 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 762 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 763 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 764 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 765 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 766 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 767 for (field = 0; field < numFields; ++field) { 768 PetscObject obj; 769 PetscClassId id; 770 PetscInt Nc; 771 772 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 773 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 774 if (id == PETSCFE_CLASSID) { 775 PetscFE fe = (PetscFE) obj; 776 777 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 778 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 779 } else if (id == PETSCFV_CLASSID) { 780 PetscFV fv = (PetscFV) obj; 781 782 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 783 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 784 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 785 numComponents += Nc; 786 } 787 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 788 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 789 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 790 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 791 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 792 for (c = cStart; c < cEnd; ++c) { 793 PetscScalar *x = NULL; 794 PetscReal elemDiff = 0.0; 795 796 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 797 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 798 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 799 800 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 801 PetscObject obj; 802 PetscClassId id; 803 void * const ctx = ctxs ? ctxs[field] : NULL; 804 PetscInt Nb, Nc, q, fc; 805 806 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 807 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 808 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 809 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 810 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 811 if (debug) { 812 char title[1024]; 813 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 814 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 815 } 816 for (q = 0; q < numQuadPoints; ++q) { 817 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 818 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx); 819 if (ierr) { 820 PetscErrorCode ierr2; 821 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 822 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 823 ierr2 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2); 824 CHKERRQ(ierr); 825 } 826 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 827 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 828 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 829 for (fc = 0; fc < Nc; ++fc) { 830 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);} 831 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 832 } 833 } 834 fieldOffset += Nb*Nc; 835 } 836 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 837 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 838 localDiff += elemDiff; 839 } 840 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 841 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 842 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 843 *diff = PetscSqrtReal(*diff); 844 PetscFunctionReturn(0); 845 } 846 847 #undef __FUNCT__ 848 #define __FUNCT__ "DMComputeL2GradientDiff_Plex" 849 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) 850 { 851 const PetscInt debug = 0; 852 PetscSection section; 853 PetscQuadrature quad; 854 Vec localX; 855 PetscScalar *funcVal, *interpolantVec; 856 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 857 PetscReal localDiff = 0.0; 858 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp; 859 PetscErrorCode ierr; 860 861 PetscFunctionBegin; 862 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 863 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 864 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 865 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 866 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 867 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 868 for (field = 0; field < numFields; ++field) { 869 PetscFE fe; 870 PetscInt Nc; 871 872 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 873 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 874 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 875 numComponents += Nc; 876 } 877 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 878 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 879 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 880 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 881 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 882 for (c = cStart; c < cEnd; ++c) { 883 PetscScalar *x = NULL; 884 PetscReal elemDiff = 0.0; 885 886 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 887 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 888 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 889 890 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 891 PetscFE fe; 892 void * const ctx = ctxs ? ctxs[field] : NULL; 893 const PetscReal *quadPoints, *quadWeights; 894 PetscReal *basisDer; 895 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 896 897 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 898 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 899 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 900 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 901 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 902 if (debug) { 903 char title[1024]; 904 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 905 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 906 } 907 for (q = 0; q < numQuadPoints; ++q) { 908 for (d = 0; d < dim; d++) { 909 coords[d] = v0[d]; 910 for (e = 0; e < dim; e++) { 911 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 912 } 913 } 914 ierr = (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx); 915 if (ierr) { 916 PetscErrorCode ierr2; 917 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 918 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 919 ierr2 = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr2); 920 CHKERRQ(ierr); 921 } 922 for (fc = 0; fc < Ncomp; ++fc) { 923 PetscScalar interpolant = 0.0; 924 925 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 926 for (f = 0; f < Nb; ++f) { 927 const PetscInt fidx = f*Ncomp+fc; 928 929 for (d = 0; d < dim; ++d) { 930 realSpaceDer[d] = 0.0; 931 for (g = 0; g < dim; ++g) { 932 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 933 } 934 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 935 } 936 } 937 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 938 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);} 939 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 940 } 941 } 942 comp += Ncomp; 943 fieldOffset += Nb*Ncomp; 944 } 945 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 946 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 947 localDiff += elemDiff; 948 } 949 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 950 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 951 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 952 *diff = PetscSqrtReal(*diff); 953 PetscFunctionReturn(0); 954 } 955 956 #undef __FUNCT__ 957 #define __FUNCT__ "DMComputeL2FieldDiff_Plex" 958 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 959 { 960 const PetscInt debug = 0; 961 PetscSection section; 962 PetscQuadrature quad; 963 Vec localX; 964 PetscScalar *funcVal, *interpolant; 965 PetscReal *coords, *v0, *J, *invJ, detJ; 966 PetscReal *localDiff; 967 const PetscReal *quadPoints, *quadWeights; 968 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 969 PetscErrorCode ierr; 970 971 PetscFunctionBegin; 972 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 973 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 974 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 975 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 976 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 977 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 978 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 979 for (field = 0; field < numFields; ++field) { 980 PetscObject obj; 981 PetscClassId id; 982 PetscInt Nc; 983 984 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 985 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 986 if (id == PETSCFE_CLASSID) { 987 PetscFE fe = (PetscFE) obj; 988 989 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 990 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 991 } else if (id == PETSCFV_CLASSID) { 992 PetscFV fv = (PetscFV) obj; 993 994 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 995 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 996 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 997 numComponents += Nc; 998 } 999 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1000 ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1001 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1002 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1003 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1004 for (c = cStart; c < cEnd; ++c) { 1005 PetscScalar *x = NULL; 1006 1007 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1008 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1009 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1010 1011 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1012 PetscObject obj; 1013 PetscClassId id; 1014 void * const ctx = ctxs ? ctxs[field] : NULL; 1015 PetscInt Nb, Nc, q, fc; 1016 1017 PetscReal elemDiff = 0.0; 1018 1019 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1020 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1021 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1022 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1023 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1024 if (debug) { 1025 char title[1024]; 1026 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1027 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 1028 } 1029 for (q = 0; q < numQuadPoints; ++q) { 1030 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 1031 ierr = (*funcs[field])(dim, time, coords, numFields, funcVal, ctx); 1032 if (ierr) { 1033 PetscErrorCode ierr2; 1034 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 1035 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 1036 ierr2 = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2); 1037 CHKERRQ(ierr); 1038 } 1039 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1040 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1041 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1042 for (fc = 0; fc < Nc; ++fc) { 1043 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);} 1044 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1045 } 1046 } 1047 fieldOffset += Nb*Nc; 1048 localDiff[field] += elemDiff; 1049 } 1050 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1051 } 1052 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1053 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1054 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 1055 ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1056 PetscFunctionReturn(0); 1057 } 1058 1059 #undef __FUNCT__ 1060 #define __FUNCT__ "DMPlexComputeL2DiffVec" 1061 /*@C 1062 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. 1063 1064 Input Parameters: 1065 + dm - The DM 1066 . time - The time 1067 . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation 1068 . ctxs - Optional array of contexts to pass to each function, or NULL. 1069 - X - The coefficient vector u_h 1070 1071 Output Parameter: 1072 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 1073 1074 Level: developer 1075 1076 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 1077 @*/ 1078 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 1079 { 1080 PetscSection section; 1081 PetscQuadrature quad; 1082 Vec localX; 1083 PetscScalar *funcVal, *interpolant; 1084 PetscReal *coords, *v0, *J, *invJ, detJ; 1085 const PetscReal *quadPoints, *quadWeights; 1086 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1087 PetscErrorCode ierr; 1088 1089 PetscFunctionBegin; 1090 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 1091 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1092 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1093 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1094 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1095 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1096 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1097 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1098 for (field = 0; field < numFields; ++field) { 1099 PetscObject obj; 1100 PetscClassId id; 1101 PetscInt Nc; 1102 1103 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1104 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1105 if (id == PETSCFE_CLASSID) { 1106 PetscFE fe = (PetscFE) obj; 1107 1108 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1109 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1110 } else if (id == PETSCFV_CLASSID) { 1111 PetscFV fv = (PetscFV) obj; 1112 1113 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1114 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1115 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1116 numComponents += Nc; 1117 } 1118 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1119 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1120 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1121 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1122 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1123 for (c = cStart; c < cEnd; ++c) { 1124 PetscScalar *x = NULL; 1125 PetscScalar elemDiff = 0.0; 1126 1127 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1128 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1129 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1130 1131 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1132 PetscObject obj; 1133 PetscClassId id; 1134 void * const ctx = ctxs ? ctxs[field] : NULL; 1135 PetscInt Nb, Nc, q, fc; 1136 1137 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1138 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1139 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1140 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1141 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1142 if (funcs[field]) { 1143 for (q = 0; q < numQuadPoints; ++q) { 1144 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 1145 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx); 1146 if (ierr) { 1147 PetscErrorCode ierr2; 1148 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 1149 ierr2 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2); 1150 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 1151 CHKERRQ(ierr); 1152 } 1153 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1154 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1155 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1156 for (fc = 0; fc < Nc; ++fc) { 1157 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1158 } 1159 } 1160 } 1161 fieldOffset += Nb*Nc; 1162 } 1163 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1164 ierr = VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);CHKERRQ(ierr); 1165 } 1166 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1167 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1168 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 1169 PetscFunctionReturn(0); 1170 } 1171 1172 #undef __FUNCT__ 1173 #define __FUNCT__ "DMPlexComputeIntegralFEM" 1174 /*@ 1175 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 1176 1177 Input Parameters: 1178 + dm - The mesh 1179 . X - Local input vector 1180 - user - The user context 1181 1182 Output Parameter: 1183 . integral - Local integral for each field 1184 1185 Level: developer 1186 1187 .seealso: DMPlexComputeResidualFEM() 1188 @*/ 1189 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 1190 { 1191 DM_Plex *mesh = (DM_Plex *) dm->data; 1192 DM dmAux; 1193 Vec localX, A; 1194 PetscDS prob, probAux = NULL; 1195 PetscSection section, sectionAux; 1196 PetscFECellGeom *cgeom; 1197 PetscScalar *u, *a = NULL; 1198 PetscReal *lintegral, *vol; 1199 PetscInt dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c; 1200 PetscInt totDim, totDimAux; 1201 PetscErrorCode ierr; 1202 1203 PetscFunctionBegin; 1204 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1205 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1206 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1207 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1208 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1209 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1210 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1211 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1212 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1213 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1214 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1215 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1216 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1217 numCells = cEnd - cStart; 1218 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1219 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1220 if (dmAux) { 1221 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1222 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1223 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1224 } 1225 ierr = PetscMalloc4(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeom,numCells,&vol);CHKERRQ(ierr); 1226 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1227 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 1228 for (c = cStart; c < cEnd; ++c) { 1229 PetscScalar *x = NULL; 1230 PetscInt i; 1231 1232 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);CHKERRQ(ierr); 1233 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol[c], NULL, NULL);CHKERRQ(ierr); 1234 if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c); 1235 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1236 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1237 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1238 if (dmAux) { 1239 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1240 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1241 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1242 } 1243 } 1244 for (f = 0; f < Nf; ++f) { 1245 PetscObject obj; 1246 PetscClassId id; 1247 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1248 1249 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1250 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1251 if (id == PETSCFE_CLASSID) { 1252 PetscFE fe = (PetscFE) obj; 1253 PetscQuadrature q; 1254 PetscInt Nq, Nb; 1255 1256 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1257 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1258 ierr = PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1259 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1260 blockSize = Nb*Nq; 1261 batchSize = numBlocks * blockSize; 1262 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1263 numChunks = numCells / (numBatches*batchSize); 1264 Ne = numChunks*numBatches*batchSize; 1265 Nr = numCells % (numBatches*batchSize); 1266 offset = numCells - Nr; 1267 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);CHKERRQ(ierr); 1268 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 1269 } else if (id == PETSCFV_CLASSID) { 1270 /* PetscFV fv = (PetscFV) obj; */ 1271 PetscInt foff; 1272 PetscPointFunc obj_func; 1273 PetscScalar lint; 1274 1275 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1276 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1277 if (obj_func) { 1278 for (c = 0; c < numCells; ++c) { 1279 /* TODO: Need full pointwise interpolation and get centroid for x argument */ 1280 obj_func(dim, Nf, 0, NULL, NULL, &u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, &lint); 1281 lintegral[f] = PetscRealPart(lint)*vol[c]; 1282 } 1283 } 1284 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1285 } 1286 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1287 if (mesh->printFEM) { 1288 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1289 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1290 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1291 } 1292 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1293 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1294 ierr = PetscFree4(lintegral,u,cgeom,vol);CHKERRQ(ierr); 1295 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1296 PetscFunctionReturn(0); 1297 } 1298 1299 #undef __FUNCT__ 1300 #define __FUNCT__ "DMPlexComputeInterpolatorNested" 1301 /*@ 1302 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1303 1304 Input Parameters: 1305 + dmf - The fine mesh 1306 . dmc - The coarse mesh 1307 - user - The user context 1308 1309 Output Parameter: 1310 . In - The interpolation matrix 1311 1312 Level: developer 1313 1314 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1315 @*/ 1316 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1317 { 1318 DM_Plex *mesh = (DM_Plex *) dmc->data; 1319 const char *name = "Interpolator"; 1320 PetscDS prob; 1321 PetscFE *feRef; 1322 PetscFV *fvRef; 1323 PetscSection fsection, fglobalSection; 1324 PetscSection csection, cglobalSection; 1325 PetscScalar *elemMat; 1326 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1327 PetscInt cTotDim, rTotDim = 0; 1328 PetscErrorCode ierr; 1329 1330 PetscFunctionBegin; 1331 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1332 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1333 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1334 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1335 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1336 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1337 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1338 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1339 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1340 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1341 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1342 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1343 for (f = 0; f < Nf; ++f) { 1344 PetscObject obj; 1345 PetscClassId id; 1346 PetscInt rNb = 0, Nc = 0; 1347 1348 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1349 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1350 if (id == PETSCFE_CLASSID) { 1351 PetscFE fe = (PetscFE) obj; 1352 1353 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1354 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1355 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1356 } else if (id == PETSCFV_CLASSID) { 1357 PetscFV fv = (PetscFV) obj; 1358 PetscDualSpace Q; 1359 1360 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1361 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1362 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1363 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1364 } 1365 rTotDim += rNb*Nc; 1366 } 1367 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1368 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1369 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1370 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1371 PetscDualSpace Qref; 1372 PetscQuadrature f; 1373 const PetscReal *qpoints, *qweights; 1374 PetscReal *points; 1375 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1376 1377 /* Compose points from all dual basis functionals */ 1378 if (feRef[fieldI]) { 1379 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1380 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1381 } else { 1382 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1383 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1384 } 1385 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1386 for (i = 0; i < fpdim; ++i) { 1387 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1388 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1389 npoints += Np; 1390 } 1391 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1392 for (i = 0, k = 0; i < fpdim; ++i) { 1393 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1394 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1395 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1396 } 1397 1398 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1399 PetscObject obj; 1400 PetscClassId id; 1401 PetscReal *B; 1402 PetscInt NcJ = 0, cpdim = 0, j; 1403 1404 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1405 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1406 if (id == PETSCFE_CLASSID) { 1407 PetscFE fe = (PetscFE) obj; 1408 1409 /* Evaluate basis at points */ 1410 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1411 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1412 /* For now, fields only interpolate themselves */ 1413 if (fieldI == fieldJ) { 1414 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); 1415 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1416 for (i = 0, k = 0; i < fpdim; ++i) { 1417 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1418 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1419 for (p = 0; p < Np; ++p, ++k) { 1420 for (j = 0; j < cpdim; ++j) { 1421 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]; 1422 } 1423 } 1424 } 1425 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1426 } 1427 } else if (id == PETSCFV_CLASSID) { 1428 PetscFV fv = (PetscFV) obj; 1429 1430 /* Evaluate constant function at points */ 1431 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1432 cpdim = 1; 1433 /* For now, fields only interpolate themselves */ 1434 if (fieldI == fieldJ) { 1435 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); 1436 for (i = 0, k = 0; i < fpdim; ++i) { 1437 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1438 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1439 for (p = 0; p < Np; ++p, ++k) { 1440 for (j = 0; j < cpdim; ++j) { 1441 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p]; 1442 } 1443 } 1444 } 1445 } 1446 } 1447 offsetJ += cpdim*NcJ; 1448 } 1449 offsetI += fpdim*Nc; 1450 ierr = PetscFree(points);CHKERRQ(ierr); 1451 } 1452 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1453 /* Preallocate matrix */ 1454 { 1455 Mat preallocator; 1456 PetscScalar *vals; 1457 PetscInt *cellCIndices, *cellFIndices; 1458 PetscInt locRows, locCols, cell; 1459 1460 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1461 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1462 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1463 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1464 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1465 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1466 for (cell = cStart; cell < cEnd; ++cell) { 1467 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1468 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1469 } 1470 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1471 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1472 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1473 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1474 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1475 } 1476 /* Fill matrix */ 1477 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1478 for (c = cStart; c < cEnd; ++c) { 1479 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1480 } 1481 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1482 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1483 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1484 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1485 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1486 if (mesh->printFEM) { 1487 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1488 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1489 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1490 } 1491 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1492 PetscFunctionReturn(0); 1493 } 1494 1495 #undef __FUNCT__ 1496 #define __FUNCT__ "DMPlexComputeInterpolatorGeneral" 1497 /*@ 1498 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1499 1500 Input Parameters: 1501 + dmf - The fine mesh 1502 . dmc - The coarse mesh 1503 - user - The user context 1504 1505 Output Parameter: 1506 . In - The interpolation matrix 1507 1508 Level: developer 1509 1510 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1511 @*/ 1512 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1513 { 1514 DM_Plex *mesh = (DM_Plex *) dmf->data; 1515 const char *name = "Interpolator"; 1516 PetscDS prob; 1517 PetscSection fsection, csection, globalFSection, globalCSection; 1518 PetscHashJK ht; 1519 PetscLayout rLayout; 1520 PetscInt *dnz, *onz; 1521 PetscInt locRows, rStart, rEnd; 1522 PetscReal *x, *v0, *J, *invJ, detJ; 1523 PetscReal *v0c, *Jc, *invJc, detJc; 1524 PetscScalar *elemMat; 1525 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1526 PetscErrorCode ierr; 1527 1528 PetscFunctionBegin; 1529 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1530 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1531 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1532 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1533 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1534 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1535 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1536 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1537 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1538 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1539 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1540 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1541 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1542 ierr = PetscMalloc1(totDim*totDim, &elemMat);CHKERRQ(ierr); 1543 1544 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1545 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1546 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1547 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1548 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1549 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1550 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1551 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1552 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1553 for (field = 0; field < Nf; ++field) { 1554 PetscObject obj; 1555 PetscClassId id; 1556 PetscDualSpace Q = NULL; 1557 PetscQuadrature f; 1558 const PetscReal *qpoints; 1559 PetscInt Nc, Np, fpdim, i, d; 1560 1561 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1562 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1563 if (id == PETSCFE_CLASSID) { 1564 PetscFE fe = (PetscFE) obj; 1565 1566 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1567 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1568 } else if (id == PETSCFV_CLASSID) { 1569 PetscFV fv = (PetscFV) obj; 1570 1571 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1572 Nc = 1; 1573 } 1574 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1575 /* For each fine grid cell */ 1576 for (cell = cStart; cell < cEnd; ++cell) { 1577 PetscInt *findices, *cindices; 1578 PetscInt numFIndices, numCIndices; 1579 1580 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1581 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1582 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1583 for (i = 0; i < fpdim; ++i) { 1584 Vec pointVec; 1585 PetscScalar *pV; 1586 IS coarseCellIS; 1587 const PetscInt *coarseCells; 1588 PetscInt numCoarseCells, q, r, c; 1589 1590 /* Get points from the dual basis functional quadrature */ 1591 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1592 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1593 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1594 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1595 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1596 for (q = 0; q < Np; ++q) { 1597 /* Transform point to real space */ 1598 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1599 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1600 } 1601 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1602 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1603 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1604 ierr = ISViewFromOptions(coarseCellIS, NULL, "-interp_is_view");CHKERRQ(ierr); 1605 /* Update preallocation info */ 1606 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1607 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1608 for (r = 0; r < Nc; ++r) { 1609 PetscHashJKKey key; 1610 PetscHashJKIter missing, iter; 1611 1612 key.j = findices[i*Nc+r]; 1613 if (key.j < 0) continue; 1614 /* Get indices for coarse elements */ 1615 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1616 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1617 for (c = 0; c < numCIndices; ++c) { 1618 key.k = cindices[c]; 1619 if (key.k < 0) continue; 1620 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1621 if (missing) { 1622 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1623 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1624 else ++onz[key.j-rStart]; 1625 } 1626 } 1627 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1628 } 1629 } 1630 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1631 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1632 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1633 } 1634 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1635 } 1636 } 1637 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1638 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1639 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1640 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1641 for (field = 0; field < Nf; ++field) { 1642 PetscObject obj; 1643 PetscClassId id; 1644 PetscDualSpace Q = NULL; 1645 PetscQuadrature f; 1646 const PetscReal *qpoints, *qweights; 1647 PetscInt Nc, Np, fpdim, i, d; 1648 1649 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1650 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1651 if (id == PETSCFE_CLASSID) { 1652 PetscFE fe = (PetscFE) obj; 1653 1654 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1655 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1656 } else if (id == PETSCFV_CLASSID) { 1657 PetscFV fv = (PetscFV) obj; 1658 1659 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1660 Nc = 1; 1661 } 1662 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1663 /* For each fine grid cell */ 1664 for (cell = cStart; cell < cEnd; ++cell) { 1665 PetscInt *findices, *cindices; 1666 PetscInt numFIndices, numCIndices; 1667 1668 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1669 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1670 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1671 for (i = 0; i < fpdim; ++i) { 1672 Vec pointVec; 1673 PetscScalar *pV; 1674 IS coarseCellIS; 1675 const PetscInt *coarseCells; 1676 PetscInt numCoarseCells, cpdim, q, c, j; 1677 1678 /* Get points from the dual basis functional quadrature */ 1679 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1680 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1681 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1682 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1683 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1684 for (q = 0; q < Np; ++q) { 1685 /* Transform point to real space */ 1686 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1687 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1688 } 1689 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1690 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1691 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1692 /* Update preallocation info */ 1693 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1694 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1695 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1696 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1697 PetscReal pVReal[3]; 1698 1699 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1700 /* Transform points from real space to coarse reference space */ 1701 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell], NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1702 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1703 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1704 1705 if (id == PETSCFE_CLASSID) { 1706 PetscFE fe = (PetscFE) obj; 1707 PetscReal *B; 1708 1709 /* Evaluate coarse basis on contained point */ 1710 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1711 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1712 /* Get elemMat entries by multiplying by weight */ 1713 for (j = 0; j < cpdim; ++j) { 1714 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell]; 1715 } 1716 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1717 } else { 1718 cpdim = 1; 1719 for (j = 0; j < cpdim; ++j) { 1720 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell]; 1721 } 1722 } 1723 /* Update interpolator */ 1724 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, Nc, numCIndices, elemMat);CHKERRQ(ierr);} 1725 ierr = MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1726 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1727 } 1728 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1729 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1730 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1731 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1732 } 1733 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1734 } 1735 } 1736 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1737 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1738 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1739 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1740 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1741 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1742 PetscFunctionReturn(0); 1743 } 1744 1745 #undef __FUNCT__ 1746 #define __FUNCT__ "DMPlexComputeInjectorFEM" 1747 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1748 { 1749 PetscDS prob; 1750 PetscFE *feRef; 1751 PetscFV *fvRef; 1752 Vec fv, cv; 1753 IS fis, cis; 1754 PetscSection fsection, fglobalSection, csection, cglobalSection; 1755 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1756 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1757 PetscErrorCode ierr; 1758 1759 PetscFunctionBegin; 1760 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1761 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1762 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1763 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1764 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1765 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1766 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1767 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1768 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1769 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1770 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1771 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1772 for (f = 0; f < Nf; ++f) { 1773 PetscObject obj; 1774 PetscClassId id; 1775 PetscInt fNb = 0, Nc = 0; 1776 1777 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1778 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1779 if (id == PETSCFE_CLASSID) { 1780 PetscFE fe = (PetscFE) obj; 1781 1782 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1783 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1784 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1785 } else if (id == PETSCFV_CLASSID) { 1786 PetscFV fv = (PetscFV) obj; 1787 PetscDualSpace Q; 1788 1789 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1790 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1791 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1792 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1793 } 1794 fTotDim += fNb*Nc; 1795 } 1796 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1797 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1798 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1799 PetscFE feC; 1800 PetscFV fvC; 1801 PetscDualSpace QF, QC; 1802 PetscInt NcF, NcC, fpdim, cpdim; 1803 1804 if (feRef[field]) { 1805 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1806 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1807 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1808 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1809 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1810 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1811 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1812 } else { 1813 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1814 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1815 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1816 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1817 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1818 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1819 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1820 } 1821 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); 1822 for (c = 0; c < cpdim; ++c) { 1823 PetscQuadrature cfunc; 1824 const PetscReal *cqpoints; 1825 PetscInt NpC; 1826 PetscBool found = PETSC_FALSE; 1827 1828 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1829 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1830 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1831 for (f = 0; f < fpdim; ++f) { 1832 PetscQuadrature ffunc; 1833 const PetscReal *fqpoints; 1834 PetscReal sum = 0.0; 1835 PetscInt NpF, comp; 1836 1837 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1838 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1839 if (NpC != NpF) continue; 1840 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1841 if (sum > 1.0e-9) continue; 1842 for (comp = 0; comp < NcC; ++comp) { 1843 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1844 } 1845 found = PETSC_TRUE; 1846 break; 1847 } 1848 if (!found) { 1849 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1850 if (fvRef[field]) { 1851 PetscInt comp; 1852 for (comp = 0; comp < NcC; ++comp) { 1853 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1854 } 1855 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1856 } 1857 } 1858 offsetC += cpdim*NcC; 1859 offsetF += fpdim*NcF; 1860 } 1861 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1862 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1863 1864 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1865 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1866 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 1867 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1868 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1869 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1870 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1871 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1872 for (c = cStart; c < cEnd; ++c) { 1873 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1874 for (d = 0; d < cTotDim; ++d) { 1875 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 1876 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]]); 1877 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1878 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1879 } 1880 } 1881 ierr = PetscFree(cmap);CHKERRQ(ierr); 1882 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1883 1884 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1885 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1886 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1887 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1888 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1889 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1890 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1891 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1892 PetscFunctionReturn(0); 1893 } 1894