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