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