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