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, off, 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 = PetscDSGetFieldOffset(prob, field, &off);CHKERRQ(ierr); 630 ierr = DMGetField(dm, field, (PetscObject *) &fv);CHKERRQ(ierr); 631 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 632 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 633 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 634 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 635 if (Grad) { 636 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 637 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 638 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 639 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 640 } 641 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 642 for (i = 0; i < numids; ++i) { 643 IS faceIS; 644 const PetscInt *faces; 645 PetscInt numFaces, f; 646 647 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 648 if (!faceIS) continue; /* No points with that id on this process */ 649 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 650 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 651 for (f = 0; f < numFaces; ++f) { 652 const PetscInt face = faces[f], *cells; 653 const PetscFVFaceGeom *fg; 654 655 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 656 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 657 if (loc >= 0) continue; 658 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 659 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 660 if (Grad) { 661 const PetscFVCellGeom *cg; 662 const PetscScalar *cx, *cgrad; 663 PetscScalar *xG; 664 PetscReal dx[3]; 665 PetscInt d; 666 667 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 668 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 669 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 670 ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr); 671 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 672 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 673 ierr = (*func)(time, fg->centroid, fg->normal, fx, &xG[off], ctx);CHKERRQ(ierr); 674 } else { 675 const PetscScalar *xI; 676 PetscScalar *xG; 677 678 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 679 ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr); 680 ierr = (*func)(time, fg->centroid, fg->normal, xI, &xG[off], ctx);CHKERRQ(ierr); 681 } 682 } 683 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 684 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 685 } 686 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 687 if (Grad) { 688 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 689 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 690 } 691 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 692 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 693 PetscFunctionReturn(0); 694 } 695 696 #undef __FUNCT__ 697 #define __FUNCT__ "DMPlexInsertBoundaryValues" 698 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 699 { 700 PetscInt numBd, b; 701 PetscErrorCode ierr; 702 703 PetscFunctionBegin; 704 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 705 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 706 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 707 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 708 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 709 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 710 for (b = 0; b < numBd; ++b) { 711 PetscBool isEssential; 712 const char *labelname; 713 DMLabel label; 714 PetscInt field; 715 PetscObject obj; 716 PetscClassId id; 717 void (*func)(); 718 PetscInt numids; 719 const PetscInt *ids; 720 void *ctx; 721 722 ierr = DMPlexGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 723 if (!isEssential) continue; 724 ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr); 725 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 726 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 727 if (id == PETSCFE_CLASSID) { 728 ierr = DMPlexInsertBoundaryValues_FEM_Internal(dm, field, label, numids, ids, (void (*)(const PetscReal[], PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 729 } else if (id == PETSCFV_CLASSID) { 730 ierr = DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, 731 field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 732 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 733 } 734 PetscFunctionReturn(0); 735 } 736 737 #undef __FUNCT__ 738 #define __FUNCT__ "DMPlexComputeL2Diff" 739 /*@C 740 DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 741 742 Input Parameters: 743 + dm - The DM 744 . funcs - The functions to evaluate for each field component 745 . ctxs - Optional array of contexts to pass to each function, or NULL. 746 - X - The coefficient vector u_h 747 748 Output Parameter: 749 . diff - The diff ||u - u_h||_2 750 751 Level: developer 752 753 .seealso: DMPlexProjectFunction(), DMPlexComputeL2FieldDiff(), DMPlexComputeL2GradientDiff() 754 @*/ 755 PetscErrorCode DMPlexComputeL2Diff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 756 { 757 const PetscInt debug = 0; 758 PetscSection section; 759 PetscQuadrature quad; 760 Vec localX; 761 PetscScalar *funcVal, *interpolant; 762 PetscReal *coords, *v0, *J, *invJ, detJ; 763 PetscReal localDiff = 0.0; 764 const PetscReal *quadPoints, *quadWeights; 765 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 766 PetscErrorCode ierr; 767 768 PetscFunctionBegin; 769 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 770 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 771 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 772 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 773 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 774 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 775 for (field = 0; field < numFields; ++field) { 776 PetscObject obj; 777 PetscClassId id; 778 PetscInt Nc; 779 780 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 781 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 782 if (id == PETSCFE_CLASSID) { 783 PetscFE fe = (PetscFE) obj; 784 785 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 786 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 787 } else if (id == PETSCFV_CLASSID) { 788 PetscFV fv = (PetscFV) obj; 789 790 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 791 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 792 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 793 numComponents += Nc; 794 } 795 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 796 ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 797 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 798 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 799 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 800 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 801 for (c = cStart; c < cEnd; ++c) { 802 PetscScalar *x = NULL; 803 PetscReal elemDiff = 0.0; 804 805 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 806 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 807 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 808 809 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 810 PetscObject obj; 811 PetscClassId id; 812 void * const ctx = ctxs ? ctxs[field] : NULL; 813 PetscInt Nb, Nc, q, fc; 814 815 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 816 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 817 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 818 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 819 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 820 if (debug) { 821 char title[1024]; 822 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 823 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 824 } 825 for (q = 0; q < numQuadPoints; ++q) { 826 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 827 (*funcs[field])(coords, funcVal, ctx); 828 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 829 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 830 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 831 for (fc = 0; fc < Nc; ++fc) { 832 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);} 833 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 834 } 835 } 836 fieldOffset += Nb*Nc; 837 } 838 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 839 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 840 localDiff += elemDiff; 841 } 842 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 843 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 844 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 845 *diff = PetscSqrtReal(*diff); 846 PetscFunctionReturn(0); 847 } 848 849 #undef __FUNCT__ 850 #define __FUNCT__ "DMPlexComputeL2GradientDiff" 851 /*@C 852 DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h. 853 854 Input Parameters: 855 + dm - The DM 856 . funcs - The gradient functions to evaluate for each field component 857 . ctxs - Optional array of contexts to pass to each function, or NULL. 858 . X - The coefficient vector u_h 859 - n - The vector to project along 860 861 Output Parameter: 862 . diff - The diff ||(grad u - grad u_h) . n||_2 863 864 Level: developer 865 866 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff() 867 @*/ 868 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 869 { 870 const PetscInt debug = 0; 871 PetscSection section; 872 PetscQuadrature quad; 873 Vec localX; 874 PetscScalar *funcVal, *interpolantVec; 875 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 876 PetscReal localDiff = 0.0; 877 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp; 878 PetscErrorCode ierr; 879 880 PetscFunctionBegin; 881 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 882 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 883 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 884 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 885 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 886 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 887 for (field = 0; field < numFields; ++field) { 888 PetscFE fe; 889 PetscInt Nc; 890 891 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 892 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 893 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 894 numComponents += Nc; 895 } 896 /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 897 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 898 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 899 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 900 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 901 for (c = cStart; c < cEnd; ++c) { 902 PetscScalar *x = NULL; 903 PetscReal elemDiff = 0.0; 904 905 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 906 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 907 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 908 909 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 910 PetscFE fe; 911 void * const ctx = ctxs ? ctxs[field] : NULL; 912 const PetscReal *quadPoints, *quadWeights; 913 PetscReal *basisDer; 914 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 915 916 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 917 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 918 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 919 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 920 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 921 if (debug) { 922 char title[1024]; 923 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 924 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 925 } 926 for (q = 0; q < numQuadPoints; ++q) { 927 for (d = 0; d < dim; d++) { 928 coords[d] = v0[d]; 929 for (e = 0; e < dim; e++) { 930 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 931 } 932 } 933 (*funcs[field])(coords, n, funcVal, ctx); 934 for (fc = 0; fc < Ncomp; ++fc) { 935 PetscScalar interpolant = 0.0; 936 937 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 938 for (f = 0; f < Nb; ++f) { 939 const PetscInt fidx = f*Ncomp+fc; 940 941 for (d = 0; d < dim; ++d) { 942 realSpaceDer[d] = 0.0; 943 for (g = 0; g < dim; ++g) { 944 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 945 } 946 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 947 } 948 } 949 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 950 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);} 951 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 952 } 953 } 954 comp += Ncomp; 955 fieldOffset += Nb*Ncomp; 956 } 957 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 958 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 959 localDiff += elemDiff; 960 } 961 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 962 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 963 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 964 *diff = PetscSqrtReal(*diff); 965 PetscFunctionReturn(0); 966 } 967 968 #undef __FUNCT__ 969 #define __FUNCT__ "DMPlexComputeL2FieldDiff" 970 /*@C 971 DMPlexComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components. 972 973 Input Parameters: 974 + dm - The DM 975 . funcs - The functions to evaluate for each field component 976 . ctxs - Optional array of contexts to pass to each function, or NULL. 977 - X - The coefficient vector u_h 978 979 Output Parameter: 980 . diff - The array of differences, ||u^f - u^f_h||_2 981 982 Level: developer 983 984 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff(), DMPlexComputeL2GradientDiff() 985 @*/ 986 PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[]) 987 { 988 const PetscInt debug = 0; 989 PetscSection section; 990 PetscQuadrature quad; 991 Vec localX; 992 PetscScalar *funcVal, *interpolant; 993 PetscReal *coords, *v0, *J, *invJ, detJ; 994 PetscReal *localDiff; 995 const PetscReal *quadPoints, *quadWeights; 996 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 997 PetscErrorCode ierr; 998 999 PetscFunctionBegin; 1000 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1001 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1002 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1003 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1004 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1005 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1006 for (field = 0; field < numFields; ++field) { 1007 PetscObject obj; 1008 PetscClassId id; 1009 PetscInt Nc; 1010 1011 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1012 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1013 if (id == PETSCFE_CLASSID) { 1014 PetscFE fe = (PetscFE) obj; 1015 1016 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1017 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1018 } else if (id == PETSCFV_CLASSID) { 1019 PetscFV fv = (PetscFV) obj; 1020 1021 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1022 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1023 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1024 numComponents += Nc; 1025 } 1026 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1027 ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1028 ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1029 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1030 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1031 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1032 for (c = cStart; c < cEnd; ++c) { 1033 PetscScalar *x = NULL; 1034 1035 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1036 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1037 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1038 1039 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1040 PetscObject obj; 1041 PetscClassId id; 1042 void * const ctx = ctxs ? ctxs[field] : NULL; 1043 PetscInt Nb, Nc, q, fc; 1044 1045 PetscReal elemDiff = 0.0; 1046 1047 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1048 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1049 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1050 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1051 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1052 if (debug) { 1053 char title[1024]; 1054 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1055 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 1056 } 1057 for (q = 0; q < numQuadPoints; ++q) { 1058 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 1059 (*funcs[field])(coords, funcVal, ctx); 1060 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1061 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1062 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1063 for (fc = 0; fc < Nc; ++fc) { 1064 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);} 1065 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1066 } 1067 } 1068 fieldOffset += Nb*Nc; 1069 localDiff[field] += elemDiff; 1070 } 1071 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1072 } 1073 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1074 ierr = MPI_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1075 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 1076 ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1077 PetscFunctionReturn(0); 1078 } 1079 1080 #undef __FUNCT__ 1081 #define __FUNCT__ "DMPlexComputeIntegralFEM" 1082 /*@ 1083 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 1084 1085 Input Parameters: 1086 + dm - The mesh 1087 . X - Local input vector 1088 - user - The user context 1089 1090 Output Parameter: 1091 . integral - Local integral for each field 1092 1093 Level: developer 1094 1095 .seealso: DMPlexComputeResidualFEM() 1096 @*/ 1097 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 1098 { 1099 DM_Plex *mesh = (DM_Plex *) dm->data; 1100 DM dmAux; 1101 Vec localX, A; 1102 PetscDS prob, probAux = NULL; 1103 PetscQuadrature q; 1104 PetscSection section, sectionAux; 1105 PetscFECellGeom *cgeom; 1106 PetscScalar *u, *a = NULL; 1107 PetscInt dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c; 1108 PetscInt totDim, totDimAux; 1109 PetscErrorCode ierr; 1110 1111 PetscFunctionBegin; 1112 /*ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);*/ 1113 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1114 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1115 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1116 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1117 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1118 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1119 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1120 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1121 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1122 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1123 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1124 numCells = cEnd - cStart; 1125 for (f = 0; f < Nf; ++f) {integral[f] = 0.0;} 1126 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1127 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1128 if (dmAux) { 1129 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1130 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1131 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1132 } 1133 ierr = DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1134 ierr = PetscMalloc2(numCells*totDim,&u,numCells,&cgeom);CHKERRQ(ierr); 1135 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1136 for (c = cStart; c < cEnd; ++c) { 1137 PetscScalar *x = NULL; 1138 PetscInt i; 1139 1140 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);CHKERRQ(ierr); 1141 if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c); 1142 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1143 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1144 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1145 if (dmAux) { 1146 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1147 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1148 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1149 } 1150 } 1151 for (f = 0; f < Nf; ++f) { 1152 PetscFE fe; 1153 PetscInt numQuadPoints, Nb; 1154 /* Conforming batches */ 1155 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1156 /* Remainder */ 1157 PetscInt Nr, offset; 1158 1159 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1160 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1161 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1162 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1163 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1164 blockSize = Nb*numQuadPoints; 1165 batchSize = numBlocks * blockSize; 1166 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1167 numChunks = numCells / (numBatches*batchSize); 1168 Ne = numChunks*numBatches*batchSize; 1169 Nr = numCells % (numBatches*batchSize); 1170 offset = numCells - Nr; 1171 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, integral);CHKERRQ(ierr); 1172 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], integral);CHKERRQ(ierr); 1173 } 1174 ierr = PetscFree2(u,cgeom);CHKERRQ(ierr); 1175 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1176 if (mesh->printFEM) { 1177 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1178 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", integral[f]);CHKERRQ(ierr);} 1179 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1180 } 1181 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1182 /* TODO: Allreduce for integral */ 1183 /*ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);*/ 1184 PetscFunctionReturn(0); 1185 } 1186 1187 #undef __FUNCT__ 1188 #define __FUNCT__ "DMPlexComputeInterpolatorFEM" 1189 /*@ 1190 DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1191 1192 Input Parameters: 1193 + dmf - The fine mesh 1194 . dmc - The coarse mesh 1195 - user - The user context 1196 1197 Output Parameter: 1198 . In - The interpolation matrix 1199 1200 Note: 1201 The first member of the user context must be an FEMContext. 1202 1203 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1204 like a GPU, or vectorize on a multicore machine. 1205 1206 Level: developer 1207 1208 .seealso: DMPlexComputeJacobianFEM() 1209 @*/ 1210 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user) 1211 { 1212 DM_Plex *mesh = (DM_Plex *) dmc->data; 1213 const char *name = "Interpolator"; 1214 PetscDS prob; 1215 PetscFE *feRef; 1216 PetscSection fsection, fglobalSection; 1217 PetscSection csection, cglobalSection; 1218 PetscScalar *elemMat; 1219 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1220 PetscInt cTotDim, rTotDim = 0; 1221 PetscErrorCode ierr; 1222 1223 PetscFunctionBegin; 1224 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1225 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1226 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1227 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1228 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1229 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1230 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1231 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1232 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1233 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1234 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1235 ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr); 1236 for (f = 0; f < Nf; ++f) { 1237 PetscFE fe; 1238 PetscInt rNb, Nc; 1239 1240 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1241 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1242 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1243 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1244 rTotDim += rNb*Nc; 1245 } 1246 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1247 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1248 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1249 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1250 PetscDualSpace Qref; 1251 PetscQuadrature f; 1252 const PetscReal *qpoints, *qweights; 1253 PetscReal *points; 1254 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1255 1256 /* Compose points from all dual basis functionals */ 1257 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1258 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1259 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1260 for (i = 0; i < fpdim; ++i) { 1261 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1262 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1263 npoints += Np; 1264 } 1265 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1266 for (i = 0, k = 0; i < fpdim; ++i) { 1267 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1268 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1269 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1270 } 1271 1272 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1273 PetscFE fe; 1274 PetscReal *B; 1275 PetscInt NcJ, cpdim, j; 1276 1277 /* Evaluate basis at points */ 1278 ierr = PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);CHKERRQ(ierr); 1279 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1280 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1281 /* For now, fields only interpolate themselves */ 1282 if (fieldI == fieldJ) { 1283 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); 1284 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1285 for (i = 0, k = 0; i < fpdim; ++i) { 1286 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1287 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1288 for (p = 0; p < Np; ++p, ++k) { 1289 for (j = 0; j < cpdim; ++j) { 1290 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]; 1291 } 1292 } 1293 } 1294 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1295 } 1296 offsetJ += cpdim*NcJ; 1297 } 1298 offsetI += fpdim*Nc; 1299 ierr = PetscFree(points);CHKERRQ(ierr); 1300 } 1301 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1302 /* Preallocate matrix */ 1303 { 1304 PetscHashJK ht; 1305 PetscLayout rLayout; 1306 PetscInt *dnz, *onz, *cellCIndices, *cellFIndices; 1307 PetscInt locRows, rStart, rEnd, cell, r; 1308 1309 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1310 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1311 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1312 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1313 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1314 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1315 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1316 ierr = PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1317 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1318 for (cell = cStart; cell < cEnd; ++cell) { 1319 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1320 for (r = 0; r < rTotDim; ++r) { 1321 PetscHashJKKey key; 1322 PetscHashJKIter missing, iter; 1323 1324 key.j = cellFIndices[r]; 1325 if (key.j < 0) continue; 1326 for (c = 0; c < cTotDim; ++c) { 1327 key.k = cellCIndices[c]; 1328 if (key.k < 0) continue; 1329 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1330 if (missing) { 1331 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1332 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1333 else ++onz[key.j-rStart]; 1334 } 1335 } 1336 } 1337 } 1338 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1339 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1340 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1341 ierr = PetscFree4(dnz,onz,cellCIndices,cellFIndices);CHKERRQ(ierr); 1342 } 1343 /* Fill matrix */ 1344 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1345 for (c = cStart; c < cEnd; ++c) { 1346 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1347 } 1348 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1349 ierr = PetscFree(feRef);CHKERRQ(ierr); 1350 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1351 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1352 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1353 if (mesh->printFEM) { 1354 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1355 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1356 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1357 } 1358 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1359 PetscFunctionReturn(0); 1360 } 1361 1362 #undef __FUNCT__ 1363 #define __FUNCT__ "DMPlexComputeInjectorFEM" 1364 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1365 { 1366 PetscDS prob; 1367 PetscFE *feRef; 1368 Vec fv, cv; 1369 IS fis, cis; 1370 PetscSection fsection, fglobalSection, csection, cglobalSection; 1371 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1372 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m; 1373 PetscErrorCode ierr; 1374 1375 PetscFunctionBegin; 1376 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1377 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1378 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1379 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1380 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1381 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1382 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1383 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1384 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1385 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1386 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1387 ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr); 1388 for (f = 0; f < Nf; ++f) { 1389 PetscFE fe; 1390 PetscInt fNb, Nc; 1391 1392 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1393 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1394 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1395 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1396 fTotDim += fNb*Nc; 1397 } 1398 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1399 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1400 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1401 PetscFE feC; 1402 PetscDualSpace QF, QC; 1403 PetscInt NcF, NcC, fpdim, cpdim; 1404 1405 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1406 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1407 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1408 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); 1409 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1410 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1411 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1412 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1413 for (c = 0; c < cpdim; ++c) { 1414 PetscQuadrature cfunc; 1415 const PetscReal *cqpoints; 1416 PetscInt NpC; 1417 1418 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1419 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1420 if (NpC != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1421 for (f = 0; f < fpdim; ++f) { 1422 PetscQuadrature ffunc; 1423 const PetscReal *fqpoints; 1424 PetscReal sum = 0.0; 1425 PetscInt NpF, comp; 1426 1427 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1428 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1429 if (NpC != NpF) continue; 1430 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1431 if (sum > 1.0e-9) continue; 1432 for (comp = 0; comp < NcC; ++comp) { 1433 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1434 } 1435 break; 1436 } 1437 } 1438 offsetC += cpdim*NcC; 1439 offsetF += fpdim*NcF; 1440 } 1441 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1442 ierr = PetscFree(feRef);CHKERRQ(ierr); 1443 1444 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1445 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1446 ierr = VecGetOwnershipRange(cv, &startC, NULL);CHKERRQ(ierr); 1447 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1448 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1449 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1450 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1451 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1452 for (c = cStart; c < cEnd; ++c) { 1453 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1454 for (d = 0; d < cTotDim; ++d) { 1455 if (cellCIndices[d] < 0) continue; 1456 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]]); 1457 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1458 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1459 } 1460 } 1461 ierr = PetscFree(cmap);CHKERRQ(ierr); 1462 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1463 1464 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1465 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1466 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1467 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1468 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1469 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1470 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1471 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1472 PetscFunctionReturn(0); 1473 } 1474