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