1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 3 #include <petscds.h> 4 #include <petsc/private/petscimpl.h> 5 #include <petsc/private/petscfeimpl.h> 6 7 /*@C 8 DMInterpolationCreate - Creates a `DMInterpolationInfo` context 9 10 Collective 11 12 Input Parameter: 13 . comm - the communicator 14 15 Output Parameter: 16 . ctx - the context 17 18 Level: beginner 19 20 Developer Note: 21 The naming is incorrect, either the object should be named `DMInterpolation` or all the routines should begin with `DMInterpolationInfo` 22 23 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationDestroy()` 24 @*/ 25 PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx) 26 { 27 PetscFunctionBegin; 28 PetscAssertPointer(ctx, 2); 29 PetscCall(PetscNew(ctx)); 30 31 (*ctx)->comm = comm; 32 (*ctx)->dim = -1; 33 (*ctx)->nInput = 0; 34 (*ctx)->points = NULL; 35 (*ctx)->cells = NULL; 36 (*ctx)->n = -1; 37 (*ctx)->coords = NULL; 38 PetscFunctionReturn(PETSC_SUCCESS); 39 } 40 41 /*@C 42 DMInterpolationSetDim - Sets the spatial dimension for the interpolation context 43 44 Not Collective 45 46 Input Parameters: 47 + ctx - the context 48 - dim - the spatial dimension 49 50 Level: intermediate 51 52 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()` 53 @*/ 54 PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim) 55 { 56 PetscFunctionBegin; 57 PetscCheck(!(dim < 1) && !(dim > 3), ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %" PetscInt_FMT, dim); 58 ctx->dim = dim; 59 PetscFunctionReturn(PETSC_SUCCESS); 60 } 61 62 /*@C 63 DMInterpolationGetDim - Gets the spatial dimension for the interpolation context 64 65 Not Collective 66 67 Input Parameter: 68 . ctx - the context 69 70 Output Parameter: 71 . dim - the spatial dimension 72 73 Level: intermediate 74 75 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()` 76 @*/ 77 PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim) 78 { 79 PetscFunctionBegin; 80 PetscAssertPointer(dim, 2); 81 *dim = ctx->dim; 82 PetscFunctionReturn(PETSC_SUCCESS); 83 } 84 85 /*@C 86 DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context 87 88 Not Collective 89 90 Input Parameters: 91 + ctx - the context 92 - dof - the number of fields 93 94 Level: intermediate 95 96 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()` 97 @*/ 98 PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof) 99 { 100 PetscFunctionBegin; 101 PetscCheck(dof >= 1, ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %" PetscInt_FMT, dof); 102 ctx->dof = dof; 103 PetscFunctionReturn(PETSC_SUCCESS); 104 } 105 106 /*@C 107 DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context 108 109 Not Collective 110 111 Input Parameter: 112 . ctx - the context 113 114 Output Parameter: 115 . dof - the number of fields 116 117 Level: intermediate 118 119 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()` 120 @*/ 121 PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof) 122 { 123 PetscFunctionBegin; 124 PetscAssertPointer(dof, 2); 125 *dof = ctx->dof; 126 PetscFunctionReturn(PETSC_SUCCESS); 127 } 128 129 /*@C 130 DMInterpolationAddPoints - Add points at which we will interpolate the fields 131 132 Not Collective 133 134 Input Parameters: 135 + ctx - the context 136 . n - the number of points 137 - points - the coordinates for each point, an array of size `n` * dim 138 139 Level: intermediate 140 141 Note: 142 The input coordinate information is copied into the object. 143 144 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationCreate()` 145 @*/ 146 PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[]) 147 { 148 PetscFunctionBegin; 149 PetscCheck(ctx->dim >= 0, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set"); 150 PetscCheck(!ctx->points, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet"); 151 ctx->nInput = n; 152 153 PetscCall(PetscMalloc1(n * ctx->dim, &ctx->points)); 154 PetscCall(PetscArraycpy(ctx->points, points, n * ctx->dim)); 155 PetscFunctionReturn(PETSC_SUCCESS); 156 } 157 158 /*@C 159 DMInterpolationSetUp - Compute spatial indices for point location during interpolation 160 161 Collective 162 163 Input Parameters: 164 + ctx - the context 165 . dm - the `DM` for the function space used for interpolation 166 . redundantPoints - If `PETSC_TRUE`, all processes are passing in the same array of points. Otherwise, points need to be communicated among processes. 167 - ignoreOutsideDomain - If `PETSC_TRUE`, ignore points outside the domain, otherwise return an error 168 169 Level: intermediate 170 171 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 172 @*/ 173 PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain) 174 { 175 MPI_Comm comm = ctx->comm; 176 PetscScalar *a; 177 PetscInt p, q, i; 178 PetscMPIInt rank, size; 179 Vec pointVec; 180 PetscSF cellSF; 181 PetscLayout layout; 182 PetscReal *globalPoints; 183 PetscScalar *globalPointsScalar; 184 const PetscInt *ranges; 185 PetscMPIInt *counts, *displs, iN; 186 const PetscSFNode *foundCells; 187 const PetscInt *foundPoints; 188 PetscMPIInt *foundProcs, *globalProcs, in; 189 PetscInt n, N, numFound; 190 191 PetscFunctionBegin; 192 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 193 PetscCallMPI(MPI_Comm_size(comm, &size)); 194 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 195 PetscCheck(ctx->dim >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set"); 196 /* Locate points */ 197 n = ctx->nInput; 198 if (!redundantPoints) { 199 PetscCall(PetscLayoutCreate(comm, &layout)); 200 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 201 PetscCall(PetscLayoutSetLocalSize(layout, n)); 202 PetscCall(PetscLayoutSetUp(layout)); 203 PetscCall(PetscLayoutGetSize(layout, &N)); 204 /* Communicate all points to all processes */ 205 PetscCall(PetscMalloc3(N * ctx->dim, &globalPoints, size, &counts, size, &displs)); 206 PetscCall(PetscLayoutGetRanges(layout, &ranges)); 207 for (p = 0; p < size; ++p) { 208 PetscCall(PetscMPIIntCast((ranges[p + 1] - ranges[p]) * ctx->dim, &counts[p])); 209 PetscCall(PetscMPIIntCast(ranges[p] * ctx->dim, &displs[p])); 210 } 211 PetscCall(PetscMPIIntCast(n * ctx->dim, &in)); 212 PetscCallMPI(MPI_Allgatherv(ctx->points, in, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm)); 213 } else { 214 N = n; 215 globalPoints = ctx->points; 216 counts = displs = NULL; 217 layout = NULL; 218 } 219 #if 0 220 PetscCall(PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs)); 221 /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */ 222 #else 223 #if defined(PETSC_USE_COMPLEX) 224 PetscCall(PetscMalloc1(N * ctx->dim, &globalPointsScalar)); 225 for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i]; 226 #else 227 globalPointsScalar = globalPoints; 228 #endif 229 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec)); 230 PetscCall(PetscMalloc2(N, &foundProcs, N, &globalProcs)); 231 for (p = 0; p < N; ++p) foundProcs[p] = size; 232 cellSF = NULL; 233 PetscCall(DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF)); 234 PetscCall(PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells)); 235 #endif 236 for (p = 0; p < numFound; ++p) { 237 if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank; 238 } 239 /* Let the lowest rank process own each point */ 240 PetscCall(PetscMPIIntCast(N, &iN)); 241 PetscCallMPI(MPIU_Allreduce(foundProcs, globalProcs, iN, MPI_INT, MPI_MIN, comm)); 242 ctx->n = 0; 243 for (p = 0; p < N; ++p) { 244 if (globalProcs[p] == size) { 245 PetscCheck(ignoreOutsideDomain, comm, PETSC_ERR_PLIB, "Point %" PetscInt_FMT ": %g %g %g not located in mesh", p, (double)globalPoints[p * ctx->dim + 0], (double)(ctx->dim > 1 ? globalPoints[p * ctx->dim + 1] : 0.0), 246 (double)(ctx->dim > 2 ? globalPoints[p * ctx->dim + 2] : 0.0)); 247 if (rank == 0) ++ctx->n; 248 } else if (globalProcs[p] == rank) ++ctx->n; 249 } 250 /* Create coordinates vector and array of owned cells */ 251 PetscCall(PetscMalloc1(ctx->n, &ctx->cells)); 252 PetscCall(VecCreate(comm, &ctx->coords)); 253 PetscCall(VecSetSizes(ctx->coords, ctx->n * ctx->dim, PETSC_DECIDE)); 254 PetscCall(VecSetBlockSize(ctx->coords, ctx->dim)); 255 PetscCall(VecSetType(ctx->coords, VECSTANDARD)); 256 PetscCall(VecGetArray(ctx->coords, &a)); 257 for (p = 0, q = 0, i = 0; p < N; ++p) { 258 if (globalProcs[p] == rank) { 259 PetscInt d; 260 261 for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p * ctx->dim + d]; 262 ctx->cells[q] = foundCells[q].index; 263 ++q; 264 } 265 if (globalProcs[p] == size && rank == 0) { 266 PetscInt d; 267 268 for (d = 0; d < ctx->dim; ++d, ++i) a[i] = 0.; 269 ctx->cells[q] = -1; 270 ++q; 271 } 272 } 273 PetscCall(VecRestoreArray(ctx->coords, &a)); 274 #if 0 275 PetscCall(PetscFree3(foundCells,foundProcs,globalProcs)); 276 #else 277 PetscCall(PetscFree2(foundProcs, globalProcs)); 278 PetscCall(PetscSFDestroy(&cellSF)); 279 PetscCall(VecDestroy(&pointVec)); 280 #endif 281 if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar)); 282 if (!redundantPoints) PetscCall(PetscFree3(globalPoints, counts, displs)); 283 PetscCall(PetscLayoutDestroy(&layout)); 284 PetscFunctionReturn(PETSC_SUCCESS); 285 } 286 287 /*@C 288 DMInterpolationGetCoordinates - Gets a `Vec` with the coordinates of each interpolation point 289 290 Collective 291 292 Input Parameter: 293 . ctx - the context 294 295 Output Parameter: 296 . coordinates - the coordinates of interpolation points 297 298 Level: intermediate 299 300 Note: 301 The local vector entries correspond to interpolation points lying on this process, according to the associated `DM`. 302 This is a borrowed vector that the user should not destroy. 303 304 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 305 @*/ 306 PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates) 307 { 308 PetscFunctionBegin; 309 PetscAssertPointer(coordinates, 2); 310 PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup."); 311 *coordinates = ctx->coords; 312 PetscFunctionReturn(PETSC_SUCCESS); 313 } 314 315 /*@C 316 DMInterpolationGetVector - Gets a `Vec` which can hold all the interpolated field values 317 318 Collective 319 320 Input Parameter: 321 . ctx - the context 322 323 Output Parameter: 324 . v - a vector capable of holding the interpolated field values 325 326 Level: intermediate 327 328 Note: 329 This vector should be returned using `DMInterpolationRestoreVector()`. 330 331 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationRestoreVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 332 @*/ 333 PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v) 334 { 335 PetscFunctionBegin; 336 PetscAssertPointer(v, 2); 337 PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup."); 338 PetscCall(VecCreate(ctx->comm, v)); 339 PetscCall(VecSetSizes(*v, ctx->n * ctx->dof, PETSC_DECIDE)); 340 PetscCall(VecSetBlockSize(*v, ctx->dof)); 341 PetscCall(VecSetType(*v, VECSTANDARD)); 342 PetscFunctionReturn(PETSC_SUCCESS); 343 } 344 345 /*@C 346 DMInterpolationRestoreVector - Returns a `Vec` which can hold all the interpolated field values 347 348 Collective 349 350 Input Parameters: 351 + ctx - the context 352 - v - a vector capable of holding the interpolated field values 353 354 Level: intermediate 355 356 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 357 @*/ 358 PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v) 359 { 360 PetscFunctionBegin; 361 PetscAssertPointer(v, 2); 362 PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup."); 363 PetscCall(VecDestroy(v)); 364 PetscFunctionReturn(PETSC_SUCCESS); 365 } 366 367 static inline PetscErrorCode DMInterpolate_Segment_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 368 { 369 const PetscInt c = ctx->cells[p]; 370 const PetscInt dof = ctx->dof; 371 PetscScalar *x = NULL; 372 const PetscScalar *coords; 373 PetscScalar *a; 374 PetscReal v0, J, invJ, detJ, xir[1]; 375 PetscInt xSize, comp; 376 377 PetscFunctionBegin; 378 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 379 PetscCall(VecGetArray(v, &a)); 380 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ)); 381 PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c); 382 xir[0] = invJ * PetscRealPart(coords[p] - v0); 383 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x)); 384 if (2 * dof == xSize) { 385 for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) + x[1 * dof + comp] * xir[0]; 386 } else if (dof == xSize) { 387 for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp]; 388 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Input closure size %" PetscInt_FMT " must be either %" PetscInt_FMT " or %" PetscInt_FMT, xSize, 2 * dof, dof); 389 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x)); 390 PetscCall(VecRestoreArray(v, &a)); 391 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 392 PetscFunctionReturn(PETSC_SUCCESS); 393 } 394 395 static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 396 { 397 const PetscInt c = ctx->cells[p]; 398 PetscScalar *x = NULL; 399 const PetscScalar *coords; 400 PetscScalar *a; 401 PetscReal *v0, *J, *invJ, detJ; 402 PetscReal xi[4]; 403 404 PetscFunctionBegin; 405 PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ)); 406 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 407 PetscCall(VecGetArray(v, &a)); 408 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 409 PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c); 410 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x)); 411 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp]; 412 for (PetscInt d = 0; d < ctx->dim; ++d) { 413 xi[d] = 0.0; 414 for (PetscInt f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]); 415 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] += PetscRealPart(x[(d + 1) * ctx->dof + comp] - x[0 * ctx->dof + comp]) * xi[d]; 416 } 417 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x)); 418 PetscCall(VecRestoreArray(v, &a)); 419 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 420 PetscCall(PetscFree3(v0, J, invJ)); 421 PetscFunctionReturn(PETSC_SUCCESS); 422 } 423 424 static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 425 { 426 const PetscInt c = ctx->cells[p]; 427 const PetscInt order[3] = {2, 1, 3}; 428 PetscScalar *x = NULL; 429 PetscReal *v0, *J, *invJ, detJ; 430 const PetscScalar *coords; 431 PetscScalar *a; 432 PetscReal xi[4]; 433 434 PetscFunctionBegin; 435 PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ)); 436 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 437 PetscCall(VecGetArray(v, &a)); 438 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 439 PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c); 440 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x)); 441 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp]; 442 for (PetscInt d = 0; d < ctx->dim; ++d) { 443 xi[d] = 0.0; 444 for (PetscInt f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]); 445 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] += PetscRealPart(x[order[d] * ctx->dof + comp] - x[0 * ctx->dof + comp]) * xi[d]; 446 } 447 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x)); 448 PetscCall(VecRestoreArray(v, &a)); 449 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 450 PetscCall(PetscFree3(v0, J, invJ)); 451 PetscFunctionReturn(PETSC_SUCCESS); 452 } 453 454 static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx) 455 { 456 const PetscScalar *vertices = (const PetscScalar *)ctx; 457 const PetscScalar x0 = vertices[0]; 458 const PetscScalar y0 = vertices[1]; 459 const PetscScalar x1 = vertices[2]; 460 const PetscScalar y1 = vertices[3]; 461 const PetscScalar x2 = vertices[4]; 462 const PetscScalar y2 = vertices[5]; 463 const PetscScalar x3 = vertices[6]; 464 const PetscScalar y3 = vertices[7]; 465 const PetscScalar f_1 = x1 - x0; 466 const PetscScalar g_1 = y1 - y0; 467 const PetscScalar f_3 = x3 - x0; 468 const PetscScalar g_3 = y3 - y0; 469 const PetscScalar f_01 = x2 - x1 - x3 + x0; 470 const PetscScalar g_01 = y2 - y1 - y3 + y0; 471 const PetscScalar *ref; 472 PetscScalar *real; 473 474 PetscFunctionBegin; 475 PetscCall(VecGetArrayRead(Xref, &ref)); 476 PetscCall(VecGetArray(Xreal, &real)); 477 { 478 const PetscScalar p0 = ref[0]; 479 const PetscScalar p1 = ref[1]; 480 481 real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1; 482 real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1; 483 } 484 PetscCall(PetscLogFlops(28)); 485 PetscCall(VecRestoreArrayRead(Xref, &ref)); 486 PetscCall(VecRestoreArray(Xreal, &real)); 487 PetscFunctionReturn(PETSC_SUCCESS); 488 } 489 490 #include <petsc/private/dmimpl.h> 491 static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx) 492 { 493 const PetscScalar *vertices = (const PetscScalar *)ctx; 494 const PetscScalar x0 = vertices[0]; 495 const PetscScalar y0 = vertices[1]; 496 const PetscScalar x1 = vertices[2]; 497 const PetscScalar y1 = vertices[3]; 498 const PetscScalar x2 = vertices[4]; 499 const PetscScalar y2 = vertices[5]; 500 const PetscScalar x3 = vertices[6]; 501 const PetscScalar y3 = vertices[7]; 502 const PetscScalar f_01 = x2 - x1 - x3 + x0; 503 const PetscScalar g_01 = y2 - y1 - y3 + y0; 504 const PetscScalar *ref; 505 506 PetscFunctionBegin; 507 PetscCall(VecGetArrayRead(Xref, &ref)); 508 { 509 const PetscScalar x = ref[0]; 510 const PetscScalar y = ref[1]; 511 const PetscInt rows[2] = {0, 1}; 512 PetscScalar values[4]; 513 514 values[0] = (x1 - x0 + f_01 * y) * 0.5; 515 values[1] = (x3 - x0 + f_01 * x) * 0.5; 516 values[2] = (y1 - y0 + g_01 * y) * 0.5; 517 values[3] = (y3 - y0 + g_01 * x) * 0.5; 518 PetscCall(MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES)); 519 } 520 PetscCall(PetscLogFlops(30)); 521 PetscCall(VecRestoreArrayRead(Xref, &ref)); 522 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 523 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 524 PetscFunctionReturn(PETSC_SUCCESS); 525 } 526 527 static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 528 { 529 const PetscInt c = ctx->cells[p]; 530 PetscFE fem = NULL; 531 DM dmCoord; 532 SNES snes; 533 KSP ksp; 534 PC pc; 535 Vec coordsLocal, r, ref, real; 536 Mat J; 537 PetscTabulation T = NULL; 538 const PetscScalar *coords; 539 PetscScalar *a; 540 PetscReal xir[2] = {0., 0.}; 541 PetscInt Nf; 542 const PetscInt dof = ctx->dof; 543 PetscScalar *x = NULL, *vertices = NULL; 544 PetscScalar *xi; 545 PetscInt coordSize, xSize; 546 547 PetscFunctionBegin; 548 PetscCall(DMGetNumFields(dm, &Nf)); 549 if (Nf) { 550 PetscObject obj; 551 PetscClassId id; 552 553 PetscCall(DMGetField(dm, 0, NULL, &obj)); 554 PetscCall(PetscObjectGetClassId(obj, &id)); 555 if (id == PETSCFE_CLASSID) { 556 fem = (PetscFE)obj; 557 PetscCall(PetscFECreateTabulation(fem, 1, 1, xir, 0, &T)); 558 } 559 } 560 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 561 PetscCall(DMGetCoordinateDM(dm, &dmCoord)); 562 PetscCall(SNESCreate(PETSC_COMM_SELF, &snes)); 563 PetscCall(SNESSetOptionsPrefix(snes, "quad_interp_")); 564 PetscCall(VecCreate(PETSC_COMM_SELF, &r)); 565 PetscCall(VecSetSizes(r, 2, 2)); 566 PetscCall(VecSetType(r, dm->vectype)); 567 PetscCall(VecDuplicate(r, &ref)); 568 PetscCall(VecDuplicate(r, &real)); 569 PetscCall(MatCreate(PETSC_COMM_SELF, &J)); 570 PetscCall(MatSetSizes(J, 2, 2, 2, 2)); 571 PetscCall(MatSetType(J, MATSEQDENSE)); 572 PetscCall(MatSetUp(J)); 573 PetscCall(SNESSetFunction(snes, r, QuadMap_Private, NULL)); 574 PetscCall(SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL)); 575 PetscCall(SNESGetKSP(snes, &ksp)); 576 PetscCall(KSPGetPC(ksp, &pc)); 577 PetscCall(PCSetType(pc, PCLU)); 578 PetscCall(SNESSetFromOptions(snes)); 579 580 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 581 PetscCall(VecGetArray(v, &a)); 582 PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 583 PetscCheck(4 * 2 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %" PetscInt_FMT " should be %d", coordSize, 4 * 2); 584 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x)); 585 PetscCall(SNESSetFunction(snes, NULL, NULL, vertices)); 586 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices)); 587 PetscCall(VecGetArray(real, &xi)); 588 xi[0] = coords[p * ctx->dim + 0]; 589 xi[1] = coords[p * ctx->dim + 1]; 590 PetscCall(VecRestoreArray(real, &xi)); 591 PetscCall(SNESSolve(snes, real, ref)); 592 PetscCall(VecGetArray(ref, &xi)); 593 xir[0] = PetscRealPart(xi[0]); 594 xir[1] = PetscRealPart(xi[1]); 595 if (4 * dof == xSize) { 596 for (PetscInt comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) * (1 - xir[1]) + x[1 * dof + comp] * xir[0] * (1 - xir[1]) + x[2 * dof + comp] * xir[0] * xir[1] + x[3 * dof + comp] * (1 - xir[0]) * xir[1]; 597 } else if (dof == xSize) { 598 for (PetscInt comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp]; 599 } else { 600 PetscCheck(fem, ctx->comm, PETSC_ERR_ARG_WRONG, "Cannot have a higher order interpolant if the discretization is not PetscFE"); 601 xir[0] = 2.0 * xir[0] - 1.0; 602 xir[1] = 2.0 * xir[1] - 1.0; 603 PetscCall(PetscFEComputeTabulation(fem, 1, xir, 0, T)); 604 for (PetscInt comp = 0; comp < dof; ++comp) { 605 a[p * dof + comp] = 0.0; 606 for (PetscInt d = 0; d < xSize / dof; ++d) a[p * dof + comp] += x[d * dof + comp] * T->T[0][d * dof + comp]; 607 } 608 } 609 PetscCall(VecRestoreArray(ref, &xi)); 610 PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 611 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x)); 612 PetscCall(PetscTabulationDestroy(&T)); 613 PetscCall(VecRestoreArray(v, &a)); 614 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 615 616 PetscCall(SNESDestroy(&snes)); 617 PetscCall(VecDestroy(&r)); 618 PetscCall(VecDestroy(&ref)); 619 PetscCall(VecDestroy(&real)); 620 PetscCall(MatDestroy(&J)); 621 PetscFunctionReturn(PETSC_SUCCESS); 622 } 623 624 static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx) 625 { 626 const PetscScalar *vertices = (const PetscScalar *)ctx; 627 const PetscScalar x0 = vertices[0]; 628 const PetscScalar y0 = vertices[1]; 629 const PetscScalar z0 = vertices[2]; 630 const PetscScalar x1 = vertices[9]; 631 const PetscScalar y1 = vertices[10]; 632 const PetscScalar z1 = vertices[11]; 633 const PetscScalar x2 = vertices[6]; 634 const PetscScalar y2 = vertices[7]; 635 const PetscScalar z2 = vertices[8]; 636 const PetscScalar x3 = vertices[3]; 637 const PetscScalar y3 = vertices[4]; 638 const PetscScalar z3 = vertices[5]; 639 const PetscScalar x4 = vertices[12]; 640 const PetscScalar y4 = vertices[13]; 641 const PetscScalar z4 = vertices[14]; 642 const PetscScalar x5 = vertices[15]; 643 const PetscScalar y5 = vertices[16]; 644 const PetscScalar z5 = vertices[17]; 645 const PetscScalar x6 = vertices[18]; 646 const PetscScalar y6 = vertices[19]; 647 const PetscScalar z6 = vertices[20]; 648 const PetscScalar x7 = vertices[21]; 649 const PetscScalar y7 = vertices[22]; 650 const PetscScalar z7 = vertices[23]; 651 const PetscScalar f_1 = x1 - x0; 652 const PetscScalar g_1 = y1 - y0; 653 const PetscScalar h_1 = z1 - z0; 654 const PetscScalar f_3 = x3 - x0; 655 const PetscScalar g_3 = y3 - y0; 656 const PetscScalar h_3 = z3 - z0; 657 const PetscScalar f_4 = x4 - x0; 658 const PetscScalar g_4 = y4 - y0; 659 const PetscScalar h_4 = z4 - z0; 660 const PetscScalar f_01 = x2 - x1 - x3 + x0; 661 const PetscScalar g_01 = y2 - y1 - y3 + y0; 662 const PetscScalar h_01 = z2 - z1 - z3 + z0; 663 const PetscScalar f_12 = x7 - x3 - x4 + x0; 664 const PetscScalar g_12 = y7 - y3 - y4 + y0; 665 const PetscScalar h_12 = z7 - z3 - z4 + z0; 666 const PetscScalar f_02 = x5 - x1 - x4 + x0; 667 const PetscScalar g_02 = y5 - y1 - y4 + y0; 668 const PetscScalar h_02 = z5 - z1 - z4 + z0; 669 const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7; 670 const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7; 671 const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7; 672 const PetscScalar *ref; 673 PetscScalar *real; 674 675 PetscFunctionBegin; 676 PetscCall(VecGetArrayRead(Xref, &ref)); 677 PetscCall(VecGetArray(Xreal, &real)); 678 { 679 const PetscScalar p0 = ref[0]; 680 const PetscScalar p1 = ref[1]; 681 const PetscScalar p2 = ref[2]; 682 683 real[0] = x0 + f_1 * p0 + f_3 * p1 + f_4 * p2 + f_01 * p0 * p1 + f_12 * p1 * p2 + f_02 * p0 * p2 + f_012 * p0 * p1 * p2; 684 real[1] = y0 + g_1 * p0 + g_3 * p1 + g_4 * p2 + g_01 * p0 * p1 + g_01 * p0 * p1 + g_12 * p1 * p2 + g_02 * p0 * p2 + g_012 * p0 * p1 * p2; 685 real[2] = z0 + h_1 * p0 + h_3 * p1 + h_4 * p2 + h_01 * p0 * p1 + h_01 * p0 * p1 + h_12 * p1 * p2 + h_02 * p0 * p2 + h_012 * p0 * p1 * p2; 686 } 687 PetscCall(PetscLogFlops(114)); 688 PetscCall(VecRestoreArrayRead(Xref, &ref)); 689 PetscCall(VecRestoreArray(Xreal, &real)); 690 PetscFunctionReturn(PETSC_SUCCESS); 691 } 692 693 static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx) 694 { 695 const PetscScalar *vertices = (const PetscScalar *)ctx; 696 const PetscScalar x0 = vertices[0]; 697 const PetscScalar y0 = vertices[1]; 698 const PetscScalar z0 = vertices[2]; 699 const PetscScalar x1 = vertices[9]; 700 const PetscScalar y1 = vertices[10]; 701 const PetscScalar z1 = vertices[11]; 702 const PetscScalar x2 = vertices[6]; 703 const PetscScalar y2 = vertices[7]; 704 const PetscScalar z2 = vertices[8]; 705 const PetscScalar x3 = vertices[3]; 706 const PetscScalar y3 = vertices[4]; 707 const PetscScalar z3 = vertices[5]; 708 const PetscScalar x4 = vertices[12]; 709 const PetscScalar y4 = vertices[13]; 710 const PetscScalar z4 = vertices[14]; 711 const PetscScalar x5 = vertices[15]; 712 const PetscScalar y5 = vertices[16]; 713 const PetscScalar z5 = vertices[17]; 714 const PetscScalar x6 = vertices[18]; 715 const PetscScalar y6 = vertices[19]; 716 const PetscScalar z6 = vertices[20]; 717 const PetscScalar x7 = vertices[21]; 718 const PetscScalar y7 = vertices[22]; 719 const PetscScalar z7 = vertices[23]; 720 const PetscScalar f_xy = x2 - x1 - x3 + x0; 721 const PetscScalar g_xy = y2 - y1 - y3 + y0; 722 const PetscScalar h_xy = z2 - z1 - z3 + z0; 723 const PetscScalar f_yz = x7 - x3 - x4 + x0; 724 const PetscScalar g_yz = y7 - y3 - y4 + y0; 725 const PetscScalar h_yz = z7 - z3 - z4 + z0; 726 const PetscScalar f_xz = x5 - x1 - x4 + x0; 727 const PetscScalar g_xz = y5 - y1 - y4 + y0; 728 const PetscScalar h_xz = z5 - z1 - z4 + z0; 729 const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7; 730 const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7; 731 const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7; 732 const PetscScalar *ref; 733 734 PetscFunctionBegin; 735 PetscCall(VecGetArrayRead(Xref, &ref)); 736 { 737 const PetscScalar x = ref[0]; 738 const PetscScalar y = ref[1]; 739 const PetscScalar z = ref[2]; 740 const PetscInt rows[3] = {0, 1, 2}; 741 PetscScalar values[9]; 742 743 values[0] = (x1 - x0 + f_xy * y + f_xz * z + f_xyz * y * z) / 2.0; 744 values[1] = (x3 - x0 + f_xy * x + f_yz * z + f_xyz * x * z) / 2.0; 745 values[2] = (x4 - x0 + f_yz * y + f_xz * x + f_xyz * x * y) / 2.0; 746 values[3] = (y1 - y0 + g_xy * y + g_xz * z + g_xyz * y * z) / 2.0; 747 values[4] = (y3 - y0 + g_xy * x + g_yz * z + g_xyz * x * z) / 2.0; 748 values[5] = (y4 - y0 + g_yz * y + g_xz * x + g_xyz * x * y) / 2.0; 749 values[6] = (z1 - z0 + h_xy * y + h_xz * z + h_xyz * y * z) / 2.0; 750 values[7] = (z3 - z0 + h_xy * x + h_yz * z + h_xyz * x * z) / 2.0; 751 values[8] = (z4 - z0 + h_yz * y + h_xz * x + h_xyz * x * y) / 2.0; 752 753 PetscCall(MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES)); 754 } 755 PetscCall(PetscLogFlops(152)); 756 PetscCall(VecRestoreArrayRead(Xref, &ref)); 757 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 758 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 759 PetscFunctionReturn(PETSC_SUCCESS); 760 } 761 762 static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 763 { 764 const PetscInt c = ctx->cells[p]; 765 DM dmCoord; 766 SNES snes; 767 KSP ksp; 768 PC pc; 769 Vec coordsLocal, r, ref, real; 770 Mat J; 771 const PetscScalar *coords; 772 PetscScalar *a; 773 PetscScalar *x = NULL, *vertices = NULL; 774 PetscScalar *xi; 775 PetscReal xir[3]; 776 PetscInt coordSize, xSize; 777 778 PetscFunctionBegin; 779 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 780 PetscCall(DMGetCoordinateDM(dm, &dmCoord)); 781 PetscCall(SNESCreate(PETSC_COMM_SELF, &snes)); 782 PetscCall(SNESSetOptionsPrefix(snes, "hex_interp_")); 783 PetscCall(VecCreate(PETSC_COMM_SELF, &r)); 784 PetscCall(VecSetSizes(r, 3, 3)); 785 PetscCall(VecSetType(r, dm->vectype)); 786 PetscCall(VecDuplicate(r, &ref)); 787 PetscCall(VecDuplicate(r, &real)); 788 PetscCall(MatCreate(PETSC_COMM_SELF, &J)); 789 PetscCall(MatSetSizes(J, 3, 3, 3, 3)); 790 PetscCall(MatSetType(J, MATSEQDENSE)); 791 PetscCall(MatSetUp(J)); 792 PetscCall(SNESSetFunction(snes, r, HexMap_Private, NULL)); 793 PetscCall(SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL)); 794 PetscCall(SNESGetKSP(snes, &ksp)); 795 PetscCall(KSPGetPC(ksp, &pc)); 796 PetscCall(PCSetType(pc, PCLU)); 797 PetscCall(SNESSetFromOptions(snes)); 798 799 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 800 PetscCall(VecGetArray(v, &a)); 801 PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 802 PetscCheck(8 * 3 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid coordinate closure size %" PetscInt_FMT " should be %d", coordSize, 8 * 3); 803 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x)); 804 PetscCheck((8 * ctx->dof == xSize) || (ctx->dof == xSize), ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input closure size %" PetscInt_FMT " should be %" PetscInt_FMT " or %" PetscInt_FMT, xSize, 8 * ctx->dof, ctx->dof); 805 PetscCall(SNESSetFunction(snes, NULL, NULL, vertices)); 806 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices)); 807 PetscCall(VecGetArray(real, &xi)); 808 xi[0] = coords[p * ctx->dim + 0]; 809 xi[1] = coords[p * ctx->dim + 1]; 810 xi[2] = coords[p * ctx->dim + 2]; 811 PetscCall(VecRestoreArray(real, &xi)); 812 PetscCall(SNESSolve(snes, real, ref)); 813 PetscCall(VecGetArray(ref, &xi)); 814 xir[0] = PetscRealPart(xi[0]); 815 xir[1] = PetscRealPart(xi[1]); 816 xir[2] = PetscRealPart(xi[2]); 817 if (8 * ctx->dof == xSize) { 818 for (PetscInt comp = 0; comp < ctx->dof; ++comp) { 819 a[p * ctx->dof + comp] = x[0 * ctx->dof + comp] * (1 - xir[0]) * (1 - xir[1]) * (1 - xir[2]) + x[3 * ctx->dof + comp] * xir[0] * (1 - xir[1]) * (1 - xir[2]) + x[2 * ctx->dof + comp] * xir[0] * xir[1] * (1 - xir[2]) + x[1 * ctx->dof + comp] * (1 - xir[0]) * xir[1] * (1 - xir[2]) + 820 x[4 * ctx->dof + comp] * (1 - xir[0]) * (1 - xir[1]) * xir[2] + x[5 * ctx->dof + comp] * xir[0] * (1 - xir[1]) * xir[2] + x[6 * ctx->dof + comp] * xir[0] * xir[1] * xir[2] + x[7 * ctx->dof + comp] * (1 - xir[0]) * xir[1] * xir[2]; 821 } 822 } else { 823 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp]; 824 } 825 PetscCall(VecRestoreArray(ref, &xi)); 826 PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 827 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x)); 828 PetscCall(VecRestoreArray(v, &a)); 829 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 830 831 PetscCall(SNESDestroy(&snes)); 832 PetscCall(VecDestroy(&r)); 833 PetscCall(VecDestroy(&ref)); 834 PetscCall(VecDestroy(&real)); 835 PetscCall(MatDestroy(&J)); 836 PetscFunctionReturn(PETSC_SUCCESS); 837 } 838 839 /*@C 840 DMInterpolationEvaluate - Using the input from `dm` and `x`, calculates interpolated field values at the interpolation points. 841 842 Input Parameters: 843 + ctx - The `DMInterpolationInfo` context obtained with `DMInterpolationCreate()` 844 . dm - The `DM` 845 - x - The local vector containing the field to be interpolated, can be created with `DMCreateGlobalVector()` 846 847 Output Parameter: 848 . v - The vector containing the interpolated values, obtained with `DMInterpolationGetVector()` 849 850 Level: beginner 851 852 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`, `DMInterpolationGetCoordinates()` 853 @*/ 854 PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v) 855 { 856 PetscDS ds; 857 PetscInt n, p, Nf, field; 858 PetscBool useDS = PETSC_FALSE; 859 860 PetscFunctionBegin; 861 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 862 PetscValidHeaderSpecific(x, VEC_CLASSID, 3); 863 PetscValidHeaderSpecific(v, VEC_CLASSID, 4); 864 PetscCall(VecGetLocalSize(v, &n)); 865 PetscCheck(n == ctx->n * ctx->dof, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input vector size %" PetscInt_FMT " should be %" PetscInt_FMT, n, ctx->n * ctx->dof); 866 if (!n) PetscFunctionReturn(PETSC_SUCCESS); 867 PetscCall(DMGetDS(dm, &ds)); 868 if (ds) { 869 useDS = PETSC_TRUE; 870 PetscCall(PetscDSGetNumFields(ds, &Nf)); 871 for (field = 0; field < Nf; ++field) { 872 PetscObject obj; 873 PetscClassId id; 874 875 PetscCall(PetscDSGetDiscretization(ds, field, &obj)); 876 PetscCall(PetscObjectGetClassId(obj, &id)); 877 if (id != PETSCFE_CLASSID && id != PETSCFV_CLASSID) { 878 useDS = PETSC_FALSE; 879 break; 880 } 881 } 882 } 883 if (useDS) { 884 const PetscScalar *coords; 885 PetscScalar *interpolant; 886 PetscInt cdim, d; 887 888 PetscCall(DMGetCoordinateDim(dm, &cdim)); 889 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 890 PetscCall(VecGetArrayWrite(v, &interpolant)); 891 for (p = 0; p < ctx->n; ++p) { 892 PetscReal pcoords[3], xi[3]; 893 PetscScalar *xa = NULL; 894 PetscInt coff = 0, foff = 0, clSize; 895 896 if (ctx->cells[p] < 0) continue; 897 for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]); 898 PetscCall(DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi)); 899 PetscCall(DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa)); 900 for (field = 0; field < Nf; ++field) { 901 PetscTabulation T; 902 PetscObject obj; 903 PetscClassId id; 904 905 PetscCall(PetscDSGetDiscretization(ds, field, &obj)); 906 PetscCall(PetscObjectGetClassId(obj, &id)); 907 if (id == PETSCFE_CLASSID) { 908 PetscFE fe = (PetscFE)obj; 909 910 PetscCall(PetscFECreateTabulation(fe, 1, 1, xi, 0, &T)); 911 { 912 const PetscReal *basis = T->T[0]; 913 const PetscInt Nb = T->Nb; 914 const PetscInt Nc = T->Nc; 915 916 for (PetscInt fc = 0; fc < Nc; ++fc) { 917 interpolant[p * ctx->dof + coff + fc] = 0.0; 918 for (PetscInt f = 0; f < Nb; ++f) interpolant[p * ctx->dof + coff + fc] += xa[foff + f] * basis[(0 * Nb + f) * Nc + fc]; 919 } 920 coff += Nc; 921 foff += Nb; 922 } 923 PetscCall(PetscTabulationDestroy(&T)); 924 } else if (id == PETSCFV_CLASSID) { 925 PetscFV fv = (PetscFV)obj; 926 PetscInt Nc; 927 928 // TODO Could use reconstruction if available 929 PetscCall(PetscFVGetNumComponents(fv, &Nc)); 930 for (PetscInt fc = 0; fc < Nc; ++fc) interpolant[p * ctx->dof + coff + fc] = xa[foff + fc]; 931 coff += Nc; 932 foff += Nc; 933 } 934 } 935 PetscCall(DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa)); 936 PetscCheck(coff == ctx->dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total components %" PetscInt_FMT " != %" PetscInt_FMT " dof specified for interpolation", coff, ctx->dof); 937 PetscCheck(foff == clSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total FE/FV space size %" PetscInt_FMT " != %" PetscInt_FMT " closure size", foff, clSize); 938 } 939 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 940 PetscCall(VecRestoreArrayWrite(v, &interpolant)); 941 } else { 942 for (PetscInt p = 0; p < ctx->n; ++p) { 943 const PetscInt cell = ctx->cells[p]; 944 DMPolytopeType ct; 945 946 PetscCall(DMPlexGetCellType(dm, cell, &ct)); 947 switch (ct) { 948 case DM_POLYTOPE_SEGMENT: 949 PetscCall(DMInterpolate_Segment_Private(ctx, dm, p, x, v)); 950 break; 951 case DM_POLYTOPE_TRIANGLE: 952 PetscCall(DMInterpolate_Triangle_Private(ctx, dm, p, x, v)); 953 break; 954 case DM_POLYTOPE_QUADRILATERAL: 955 PetscCall(DMInterpolate_Quad_Private(ctx, dm, p, x, v)); 956 break; 957 case DM_POLYTOPE_TETRAHEDRON: 958 PetscCall(DMInterpolate_Tetrahedron_Private(ctx, dm, p, x, v)); 959 break; 960 case DM_POLYTOPE_HEXAHEDRON: 961 PetscCall(DMInterpolate_Hex_Private(ctx, dm, cell, x, v)); 962 break; 963 default: 964 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for cell type %s", DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]); 965 } 966 } 967 } 968 PetscFunctionReturn(PETSC_SUCCESS); 969 } 970 971 /*@C 972 DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context 973 974 Collective 975 976 Input Parameter: 977 . ctx - the context 978 979 Level: beginner 980 981 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 982 @*/ 983 PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx) 984 { 985 PetscFunctionBegin; 986 PetscAssertPointer(ctx, 1); 987 PetscCall(VecDestroy(&(*ctx)->coords)); 988 PetscCall(PetscFree((*ctx)->points)); 989 PetscCall(PetscFree((*ctx)->cells)); 990 PetscCall(PetscFree(*ctx)); 991 *ctx = NULL; 992 PetscFunctionReturn(PETSC_SUCCESS); 993 } 994