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