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