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