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 PetscSection fsection, fglobalSection; 1240 PetscSection csection, cglobalSection; 1241 PetscScalar *elemMat; 1242 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1243 PetscInt cTotDim, rTotDim = 0; 1244 PetscErrorCode ierr; 1245 1246 PetscFunctionBegin; 1247 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1248 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1249 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1250 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1251 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1252 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1253 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1254 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1255 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1256 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1257 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1258 ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr); 1259 for (f = 0; f < Nf; ++f) { 1260 PetscFE fe; 1261 PetscInt rNb, Nc; 1262 1263 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1264 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1265 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1266 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1267 rTotDim += rNb*Nc; 1268 } 1269 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1270 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1271 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1272 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1273 PetscDualSpace Qref; 1274 PetscQuadrature f; 1275 const PetscReal *qpoints, *qweights; 1276 PetscReal *points; 1277 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1278 1279 /* Compose points from all dual basis functionals */ 1280 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1281 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1282 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1283 for (i = 0; i < fpdim; ++i) { 1284 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1285 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1286 npoints += Np; 1287 } 1288 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1289 for (i = 0, k = 0; i < fpdim; ++i) { 1290 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1291 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1292 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1293 } 1294 1295 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1296 PetscFE fe; 1297 PetscReal *B; 1298 PetscInt NcJ, cpdim, j; 1299 1300 /* Evaluate basis at points */ 1301 ierr = PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);CHKERRQ(ierr); 1302 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1303 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1304 /* For now, fields only interpolate themselves */ 1305 if (fieldI == fieldJ) { 1306 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); 1307 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1308 for (i = 0, k = 0; i < fpdim; ++i) { 1309 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1310 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1311 for (p = 0; p < Np; ++p, ++k) { 1312 for (j = 0; j < cpdim; ++j) { 1313 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]; 1314 } 1315 } 1316 } 1317 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1318 } 1319 offsetJ += cpdim*NcJ; 1320 } 1321 offsetI += fpdim*Nc; 1322 ierr = PetscFree(points);CHKERRQ(ierr); 1323 } 1324 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1325 /* Preallocate matrix */ 1326 { 1327 PetscHashJK ht; 1328 PetscLayout rLayout; 1329 PetscInt *dnz, *onz, *cellCIndices, *cellFIndices; 1330 PetscInt locRows, rStart, rEnd, cell, r; 1331 1332 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1333 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1334 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1335 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1336 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1337 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1338 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1339 ierr = PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1340 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1341 for (cell = cStart; cell < cEnd; ++cell) { 1342 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1343 for (r = 0; r < rTotDim; ++r) { 1344 PetscHashJKKey key; 1345 PetscHashJKIter missing, iter; 1346 1347 key.j = cellFIndices[r]; 1348 if (key.j < 0) continue; 1349 for (c = 0; c < cTotDim; ++c) { 1350 key.k = cellCIndices[c]; 1351 if (key.k < 0) continue; 1352 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1353 if (missing) { 1354 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1355 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1356 else ++onz[key.j-rStart]; 1357 } 1358 } 1359 } 1360 } 1361 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1362 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1363 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1364 ierr = PetscFree4(dnz,onz,cellCIndices,cellFIndices);CHKERRQ(ierr); 1365 } 1366 /* Fill matrix */ 1367 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1368 for (c = cStart; c < cEnd; ++c) { 1369 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1370 } 1371 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1372 ierr = PetscFree(feRef);CHKERRQ(ierr); 1373 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1374 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1375 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1376 if (mesh->printFEM) { 1377 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1378 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1379 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1380 } 1381 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1382 PetscFunctionReturn(0); 1383 } 1384 1385 #undef __FUNCT__ 1386 #define __FUNCT__ "DMPlexComputeInjectorFEM" 1387 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1388 { 1389 PetscDS prob; 1390 PetscFE *feRef; 1391 Vec fv, cv; 1392 IS fis, cis; 1393 PetscSection fsection, fglobalSection, csection, cglobalSection; 1394 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1395 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m; 1396 PetscErrorCode ierr; 1397 1398 PetscFunctionBegin; 1399 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1400 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1401 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1402 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1403 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1404 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1405 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1406 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1407 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1408 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1409 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1410 ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr); 1411 for (f = 0; f < Nf; ++f) { 1412 PetscFE fe; 1413 PetscInt fNb, Nc; 1414 1415 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1416 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1417 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1418 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1419 fTotDim += fNb*Nc; 1420 } 1421 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1422 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1423 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1424 PetscFE feC; 1425 PetscDualSpace QF, QC; 1426 PetscInt NcF, NcC, fpdim, cpdim; 1427 1428 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1429 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1430 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1431 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); 1432 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1433 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1434 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1435 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1436 for (c = 0; c < cpdim; ++c) { 1437 PetscQuadrature cfunc; 1438 const PetscReal *cqpoints; 1439 PetscInt NpC; 1440 1441 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1442 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1443 if (NpC != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1444 for (f = 0; f < fpdim; ++f) { 1445 PetscQuadrature ffunc; 1446 const PetscReal *fqpoints; 1447 PetscReal sum = 0.0; 1448 PetscInt NpF, comp; 1449 1450 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1451 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1452 if (NpC != NpF) continue; 1453 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1454 if (sum > 1.0e-9) continue; 1455 for (comp = 0; comp < NcC; ++comp) { 1456 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1457 } 1458 break; 1459 } 1460 } 1461 offsetC += cpdim*NcC; 1462 offsetF += fpdim*NcF; 1463 } 1464 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1465 ierr = PetscFree(feRef);CHKERRQ(ierr); 1466 1467 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1468 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1469 ierr = VecGetOwnershipRange(cv, &startC, NULL);CHKERRQ(ierr); 1470 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1471 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1472 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1473 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1474 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1475 for (c = cStart; c < cEnd; ++c) { 1476 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1477 for (d = 0; d < cTotDim; ++d) { 1478 if (cellCIndices[d] < 0) continue; 1479 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]]); 1480 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1481 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1482 } 1483 } 1484 ierr = PetscFree(cmap);CHKERRQ(ierr); 1485 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1486 1487 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1488 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1489 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1490 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1491 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1492 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1493 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1494 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1495 PetscFunctionReturn(0); 1496 } 1497