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