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