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; 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 PetscCallMPI(MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm)); 241 ctx->n = 0; 242 for (p = 0; p < N; ++p) { 243 if (globalProcs[p] == size) { 244 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), 245 (double)(ctx->dim > 2 ? globalPoints[p * ctx->dim + 2] : 0)); 246 if (rank == 0) ++ctx->n; 247 } else if (globalProcs[p] == rank) ++ctx->n; 248 } 249 /* Create coordinates vector and array of owned cells */ 250 PetscCall(PetscMalloc1(ctx->n, &ctx->cells)); 251 PetscCall(VecCreate(comm, &ctx->coords)); 252 PetscCall(VecSetSizes(ctx->coords, ctx->n * ctx->dim, PETSC_DECIDE)); 253 PetscCall(VecSetBlockSize(ctx->coords, ctx->dim)); 254 PetscCall(VecSetType(ctx->coords, VECSTANDARD)); 255 PetscCall(VecGetArray(ctx->coords, &a)); 256 for (p = 0, q = 0, i = 0; p < N; ++p) { 257 if (globalProcs[p] == rank) { 258 PetscInt d; 259 260 for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p * ctx->dim + d]; 261 ctx->cells[q] = foundCells[q].index; 262 ++q; 263 } 264 if (globalProcs[p] == size && rank == 0) { 265 PetscInt d; 266 267 for (d = 0; d < ctx->dim; ++d, ++i) a[i] = 0.; 268 ctx->cells[q] = -1; 269 ++q; 270 } 271 } 272 PetscCall(VecRestoreArray(ctx->coords, &a)); 273 #if 0 274 PetscCall(PetscFree3(foundCells,foundProcs,globalProcs)); 275 #else 276 PetscCall(PetscFree2(foundProcs, globalProcs)); 277 PetscCall(PetscSFDestroy(&cellSF)); 278 PetscCall(VecDestroy(&pointVec)); 279 #endif 280 if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar)); 281 if (!redundantPoints) PetscCall(PetscFree3(globalPoints, counts, displs)); 282 PetscCall(PetscLayoutDestroy(&layout)); 283 PetscFunctionReturn(PETSC_SUCCESS); 284 } 285 286 /*@C 287 DMInterpolationGetCoordinates - Gets a `Vec` with the coordinates of each interpolation point 288 289 Collective 290 291 Input Parameter: 292 . ctx - the context 293 294 Output Parameter: 295 . coordinates - the coordinates of interpolation points 296 297 Level: intermediate 298 299 Note: 300 The local vector entries correspond to interpolation points lying on this process, according to the associated `DM`. 301 This is a borrowed vector that the user should not destroy. 302 303 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 304 @*/ 305 PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates) 306 { 307 PetscFunctionBegin; 308 PetscAssertPointer(coordinates, 2); 309 PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup."); 310 *coordinates = ctx->coords; 311 PetscFunctionReturn(PETSC_SUCCESS); 312 } 313 314 /*@C 315 DMInterpolationGetVector - Gets a `Vec` which can hold all the interpolated field values 316 317 Collective 318 319 Input Parameter: 320 . ctx - the context 321 322 Output Parameter: 323 . v - a vector capable of holding the interpolated field values 324 325 Level: intermediate 326 327 Note: 328 This vector should be returned using `DMInterpolationRestoreVector()`. 329 330 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationRestoreVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 331 @*/ 332 PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v) 333 { 334 PetscFunctionBegin; 335 PetscAssertPointer(v, 2); 336 PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup."); 337 PetscCall(VecCreate(ctx->comm, v)); 338 PetscCall(VecSetSizes(*v, ctx->n * ctx->dof, PETSC_DECIDE)); 339 PetscCall(VecSetBlockSize(*v, ctx->dof)); 340 PetscCall(VecSetType(*v, VECSTANDARD)); 341 PetscFunctionReturn(PETSC_SUCCESS); 342 } 343 344 /*@C 345 DMInterpolationRestoreVector - Returns a `Vec` which can hold all the interpolated field values 346 347 Collective 348 349 Input Parameters: 350 + ctx - the context 351 - v - a vector capable of holding the interpolated field values 352 353 Level: intermediate 354 355 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 356 @*/ 357 PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v) 358 { 359 PetscFunctionBegin; 360 PetscAssertPointer(v, 2); 361 PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup."); 362 PetscCall(VecDestroy(v)); 363 PetscFunctionReturn(PETSC_SUCCESS); 364 } 365 366 static inline PetscErrorCode DMInterpolate_Segment_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 367 { 368 const PetscInt c = ctx->cells[p]; 369 const PetscInt dof = ctx->dof; 370 PetscScalar *x = NULL; 371 const PetscScalar *coords; 372 PetscScalar *a; 373 PetscReal v0, J, invJ, detJ, xir[1]; 374 PetscInt xSize, comp; 375 376 PetscFunctionBegin; 377 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 378 PetscCall(VecGetArray(v, &a)); 379 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ)); 380 PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c); 381 xir[0] = invJ * PetscRealPart(coords[p] - v0); 382 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x)); 383 if (2 * dof == xSize) { 384 for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) + x[1 * dof + comp] * xir[0]; 385 } else if (dof == xSize) { 386 for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp]; 387 } 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); 388 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x)); 389 PetscCall(VecRestoreArray(v, &a)); 390 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 391 PetscFunctionReturn(PETSC_SUCCESS); 392 } 393 394 static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 395 { 396 const PetscInt c = ctx->cells[p]; 397 PetscScalar *x = NULL; 398 const PetscScalar *coords; 399 PetscScalar *a; 400 PetscReal *v0, *J, *invJ, detJ; 401 PetscReal xi[4]; 402 403 PetscFunctionBegin; 404 PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ)); 405 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 406 PetscCall(VecGetArray(v, &a)); 407 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 408 PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c); 409 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x)); 410 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp]; 411 for (PetscInt d = 0; d < ctx->dim; ++d) { 412 xi[d] = 0.0; 413 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]); 414 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]; 415 } 416 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x)); 417 PetscCall(VecRestoreArray(v, &a)); 418 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 419 PetscCall(PetscFree3(v0, J, invJ)); 420 PetscFunctionReturn(PETSC_SUCCESS); 421 } 422 423 static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 424 { 425 const PetscInt c = ctx->cells[p]; 426 const PetscInt order[3] = {2, 1, 3}; 427 PetscScalar *x = NULL; 428 PetscReal *v0, *J, *invJ, detJ; 429 const PetscScalar *coords; 430 PetscScalar *a; 431 PetscReal xi[4]; 432 433 PetscFunctionBegin; 434 PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ)); 435 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 436 PetscCall(VecGetArray(v, &a)); 437 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); 438 PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c); 439 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x)); 440 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp]; 441 for (PetscInt d = 0; d < ctx->dim; ++d) { 442 xi[d] = 0.0; 443 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]); 444 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]; 445 } 446 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x)); 447 PetscCall(VecRestoreArray(v, &a)); 448 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 449 PetscCall(PetscFree3(v0, J, invJ)); 450 PetscFunctionReturn(PETSC_SUCCESS); 451 } 452 453 static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx) 454 { 455 const PetscScalar *vertices = (const PetscScalar *)ctx; 456 const PetscScalar x0 = vertices[0]; 457 const PetscScalar y0 = vertices[1]; 458 const PetscScalar x1 = vertices[2]; 459 const PetscScalar y1 = vertices[3]; 460 const PetscScalar x2 = vertices[4]; 461 const PetscScalar y2 = vertices[5]; 462 const PetscScalar x3 = vertices[6]; 463 const PetscScalar y3 = vertices[7]; 464 const PetscScalar f_1 = x1 - x0; 465 const PetscScalar g_1 = y1 - y0; 466 const PetscScalar f_3 = x3 - x0; 467 const PetscScalar g_3 = y3 - y0; 468 const PetscScalar f_01 = x2 - x1 - x3 + x0; 469 const PetscScalar g_01 = y2 - y1 - y3 + y0; 470 const PetscScalar *ref; 471 PetscScalar *real; 472 473 PetscFunctionBegin; 474 PetscCall(VecGetArrayRead(Xref, &ref)); 475 PetscCall(VecGetArray(Xreal, &real)); 476 { 477 const PetscScalar p0 = ref[0]; 478 const PetscScalar p1 = ref[1]; 479 480 real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1; 481 real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1; 482 } 483 PetscCall(PetscLogFlops(28)); 484 PetscCall(VecRestoreArrayRead(Xref, &ref)); 485 PetscCall(VecRestoreArray(Xreal, &real)); 486 PetscFunctionReturn(PETSC_SUCCESS); 487 } 488 489 #include <petsc/private/dmimpl.h> 490 static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx) 491 { 492 const PetscScalar *vertices = (const PetscScalar *)ctx; 493 const PetscScalar x0 = vertices[0]; 494 const PetscScalar y0 = vertices[1]; 495 const PetscScalar x1 = vertices[2]; 496 const PetscScalar y1 = vertices[3]; 497 const PetscScalar x2 = vertices[4]; 498 const PetscScalar y2 = vertices[5]; 499 const PetscScalar x3 = vertices[6]; 500 const PetscScalar y3 = vertices[7]; 501 const PetscScalar f_01 = x2 - x1 - x3 + x0; 502 const PetscScalar g_01 = y2 - y1 - y3 + y0; 503 const PetscScalar *ref; 504 505 PetscFunctionBegin; 506 PetscCall(VecGetArrayRead(Xref, &ref)); 507 { 508 const PetscScalar x = ref[0]; 509 const PetscScalar y = ref[1]; 510 const PetscInt rows[2] = {0, 1}; 511 PetscScalar values[4]; 512 513 values[0] = (x1 - x0 + f_01 * y) * 0.5; 514 values[1] = (x3 - x0 + f_01 * x) * 0.5; 515 values[2] = (y1 - y0 + g_01 * y) * 0.5; 516 values[3] = (y3 - y0 + g_01 * x) * 0.5; 517 PetscCall(MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES)); 518 } 519 PetscCall(PetscLogFlops(30)); 520 PetscCall(VecRestoreArrayRead(Xref, &ref)); 521 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 522 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 523 PetscFunctionReturn(PETSC_SUCCESS); 524 } 525 526 static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 527 { 528 const PetscInt c = ctx->cells[p]; 529 PetscFE fem = NULL; 530 DM dmCoord; 531 SNES snes; 532 KSP ksp; 533 PC pc; 534 Vec coordsLocal, r, ref, real; 535 Mat J; 536 PetscTabulation T = NULL; 537 const PetscScalar *coords; 538 PetscScalar *a; 539 PetscReal xir[2] = {0., 0.}; 540 PetscInt Nf; 541 const PetscInt dof = ctx->dof; 542 PetscScalar *x = NULL, *vertices = NULL; 543 PetscScalar *xi; 544 PetscInt coordSize, xSize; 545 546 PetscFunctionBegin; 547 PetscCall(DMGetNumFields(dm, &Nf)); 548 if (Nf) { 549 PetscObject obj; 550 PetscClassId id; 551 552 PetscCall(DMGetField(dm, 0, NULL, &obj)); 553 PetscCall(PetscObjectGetClassId(obj, &id)); 554 if (id == PETSCFE_CLASSID) { 555 fem = (PetscFE)obj; 556 PetscCall(PetscFECreateTabulation(fem, 1, 1, xir, 0, &T)); 557 } 558 } 559 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 560 PetscCall(DMGetCoordinateDM(dm, &dmCoord)); 561 PetscCall(SNESCreate(PETSC_COMM_SELF, &snes)); 562 PetscCall(SNESSetOptionsPrefix(snes, "quad_interp_")); 563 PetscCall(VecCreate(PETSC_COMM_SELF, &r)); 564 PetscCall(VecSetSizes(r, 2, 2)); 565 PetscCall(VecSetType(r, dm->vectype)); 566 PetscCall(VecDuplicate(r, &ref)); 567 PetscCall(VecDuplicate(r, &real)); 568 PetscCall(MatCreate(PETSC_COMM_SELF, &J)); 569 PetscCall(MatSetSizes(J, 2, 2, 2, 2)); 570 PetscCall(MatSetType(J, MATSEQDENSE)); 571 PetscCall(MatSetUp(J)); 572 PetscCall(SNESSetFunction(snes, r, QuadMap_Private, NULL)); 573 PetscCall(SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL)); 574 PetscCall(SNESGetKSP(snes, &ksp)); 575 PetscCall(KSPGetPC(ksp, &pc)); 576 PetscCall(PCSetType(pc, PCLU)); 577 PetscCall(SNESSetFromOptions(snes)); 578 579 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 580 PetscCall(VecGetArray(v, &a)); 581 PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 582 PetscCheck(4 * 2 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %" PetscInt_FMT " should be %d", coordSize, 4 * 2); 583 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x)); 584 PetscCall(SNESSetFunction(snes, NULL, NULL, vertices)); 585 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices)); 586 PetscCall(VecGetArray(real, &xi)); 587 xi[0] = coords[p * ctx->dim + 0]; 588 xi[1] = coords[p * ctx->dim + 1]; 589 PetscCall(VecRestoreArray(real, &xi)); 590 PetscCall(SNESSolve(snes, real, ref)); 591 PetscCall(VecGetArray(ref, &xi)); 592 xir[0] = PetscRealPart(xi[0]); 593 xir[1] = PetscRealPart(xi[1]); 594 if (4 * dof == xSize) { 595 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]; 596 } else if (dof == xSize) { 597 for (PetscInt comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp]; 598 } else { 599 PetscCheck(fem, ctx->comm, PETSC_ERR_ARG_WRONG, "Cannot have a higher order interpolant if the discretization is not PetscFE"); 600 xir[0] = 2.0 * xir[0] - 1.0; 601 xir[1] = 2.0 * xir[1] - 1.0; 602 PetscCall(PetscFEComputeTabulation(fem, 1, xir, 0, T)); 603 for (PetscInt comp = 0; comp < dof; ++comp) { 604 a[p * dof + comp] = 0.0; 605 for (PetscInt d = 0; d < xSize / dof; ++d) a[p * dof + comp] += x[d * dof + comp] * T->T[0][d * dof + comp]; 606 } 607 } 608 PetscCall(VecRestoreArray(ref, &xi)); 609 PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 610 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x)); 611 PetscCall(PetscTabulationDestroy(&T)); 612 PetscCall(VecRestoreArray(v, &a)); 613 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 614 615 PetscCall(SNESDestroy(&snes)); 616 PetscCall(VecDestroy(&r)); 617 PetscCall(VecDestroy(&ref)); 618 PetscCall(VecDestroy(&real)); 619 PetscCall(MatDestroy(&J)); 620 PetscFunctionReturn(PETSC_SUCCESS); 621 } 622 623 static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx) 624 { 625 const PetscScalar *vertices = (const PetscScalar *)ctx; 626 const PetscScalar x0 = vertices[0]; 627 const PetscScalar y0 = vertices[1]; 628 const PetscScalar z0 = vertices[2]; 629 const PetscScalar x1 = vertices[9]; 630 const PetscScalar y1 = vertices[10]; 631 const PetscScalar z1 = vertices[11]; 632 const PetscScalar x2 = vertices[6]; 633 const PetscScalar y2 = vertices[7]; 634 const PetscScalar z2 = vertices[8]; 635 const PetscScalar x3 = vertices[3]; 636 const PetscScalar y3 = vertices[4]; 637 const PetscScalar z3 = vertices[5]; 638 const PetscScalar x4 = vertices[12]; 639 const PetscScalar y4 = vertices[13]; 640 const PetscScalar z4 = vertices[14]; 641 const PetscScalar x5 = vertices[15]; 642 const PetscScalar y5 = vertices[16]; 643 const PetscScalar z5 = vertices[17]; 644 const PetscScalar x6 = vertices[18]; 645 const PetscScalar y6 = vertices[19]; 646 const PetscScalar z6 = vertices[20]; 647 const PetscScalar x7 = vertices[21]; 648 const PetscScalar y7 = vertices[22]; 649 const PetscScalar z7 = vertices[23]; 650 const PetscScalar f_1 = x1 - x0; 651 const PetscScalar g_1 = y1 - y0; 652 const PetscScalar h_1 = z1 - z0; 653 const PetscScalar f_3 = x3 - x0; 654 const PetscScalar g_3 = y3 - y0; 655 const PetscScalar h_3 = z3 - z0; 656 const PetscScalar f_4 = x4 - x0; 657 const PetscScalar g_4 = y4 - y0; 658 const PetscScalar h_4 = z4 - z0; 659 const PetscScalar f_01 = x2 - x1 - x3 + x0; 660 const PetscScalar g_01 = y2 - y1 - y3 + y0; 661 const PetscScalar h_01 = z2 - z1 - z3 + z0; 662 const PetscScalar f_12 = x7 - x3 - x4 + x0; 663 const PetscScalar g_12 = y7 - y3 - y4 + y0; 664 const PetscScalar h_12 = z7 - z3 - z4 + z0; 665 const PetscScalar f_02 = x5 - x1 - x4 + x0; 666 const PetscScalar g_02 = y5 - y1 - y4 + y0; 667 const PetscScalar h_02 = z5 - z1 - z4 + z0; 668 const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7; 669 const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7; 670 const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7; 671 const PetscScalar *ref; 672 PetscScalar *real; 673 674 PetscFunctionBegin; 675 PetscCall(VecGetArrayRead(Xref, &ref)); 676 PetscCall(VecGetArray(Xreal, &real)); 677 { 678 const PetscScalar p0 = ref[0]; 679 const PetscScalar p1 = ref[1]; 680 const PetscScalar p2 = ref[2]; 681 682 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; 683 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; 684 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; 685 } 686 PetscCall(PetscLogFlops(114)); 687 PetscCall(VecRestoreArrayRead(Xref, &ref)); 688 PetscCall(VecRestoreArray(Xreal, &real)); 689 PetscFunctionReturn(PETSC_SUCCESS); 690 } 691 692 static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx) 693 { 694 const PetscScalar *vertices = (const PetscScalar *)ctx; 695 const PetscScalar x0 = vertices[0]; 696 const PetscScalar y0 = vertices[1]; 697 const PetscScalar z0 = vertices[2]; 698 const PetscScalar x1 = vertices[9]; 699 const PetscScalar y1 = vertices[10]; 700 const PetscScalar z1 = vertices[11]; 701 const PetscScalar x2 = vertices[6]; 702 const PetscScalar y2 = vertices[7]; 703 const PetscScalar z2 = vertices[8]; 704 const PetscScalar x3 = vertices[3]; 705 const PetscScalar y3 = vertices[4]; 706 const PetscScalar z3 = vertices[5]; 707 const PetscScalar x4 = vertices[12]; 708 const PetscScalar y4 = vertices[13]; 709 const PetscScalar z4 = vertices[14]; 710 const PetscScalar x5 = vertices[15]; 711 const PetscScalar y5 = vertices[16]; 712 const PetscScalar z5 = vertices[17]; 713 const PetscScalar x6 = vertices[18]; 714 const PetscScalar y6 = vertices[19]; 715 const PetscScalar z6 = vertices[20]; 716 const PetscScalar x7 = vertices[21]; 717 const PetscScalar y7 = vertices[22]; 718 const PetscScalar z7 = vertices[23]; 719 const PetscScalar f_xy = x2 - x1 - x3 + x0; 720 const PetscScalar g_xy = y2 - y1 - y3 + y0; 721 const PetscScalar h_xy = z2 - z1 - z3 + z0; 722 const PetscScalar f_yz = x7 - x3 - x4 + x0; 723 const PetscScalar g_yz = y7 - y3 - y4 + y0; 724 const PetscScalar h_yz = z7 - z3 - z4 + z0; 725 const PetscScalar f_xz = x5 - x1 - x4 + x0; 726 const PetscScalar g_xz = y5 - y1 - y4 + y0; 727 const PetscScalar h_xz = z5 - z1 - z4 + z0; 728 const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7; 729 const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7; 730 const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7; 731 const PetscScalar *ref; 732 733 PetscFunctionBegin; 734 PetscCall(VecGetArrayRead(Xref, &ref)); 735 { 736 const PetscScalar x = ref[0]; 737 const PetscScalar y = ref[1]; 738 const PetscScalar z = ref[2]; 739 const PetscInt rows[3] = {0, 1, 2}; 740 PetscScalar values[9]; 741 742 values[0] = (x1 - x0 + f_xy * y + f_xz * z + f_xyz * y * z) / 2.0; 743 values[1] = (x3 - x0 + f_xy * x + f_yz * z + f_xyz * x * z) / 2.0; 744 values[2] = (x4 - x0 + f_yz * y + f_xz * x + f_xyz * x * y) / 2.0; 745 values[3] = (y1 - y0 + g_xy * y + g_xz * z + g_xyz * y * z) / 2.0; 746 values[4] = (y3 - y0 + g_xy * x + g_yz * z + g_xyz * x * z) / 2.0; 747 values[5] = (y4 - y0 + g_yz * y + g_xz * x + g_xyz * x * y) / 2.0; 748 values[6] = (z1 - z0 + h_xy * y + h_xz * z + h_xyz * y * z) / 2.0; 749 values[7] = (z3 - z0 + h_xy * x + h_yz * z + h_xyz * x * z) / 2.0; 750 values[8] = (z4 - z0 + h_yz * y + h_xz * x + h_xyz * x * y) / 2.0; 751 752 PetscCall(MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES)); 753 } 754 PetscCall(PetscLogFlops(152)); 755 PetscCall(VecRestoreArrayRead(Xref, &ref)); 756 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 757 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 758 PetscFunctionReturn(PETSC_SUCCESS); 759 } 760 761 static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v) 762 { 763 const PetscInt c = ctx->cells[p]; 764 DM dmCoord; 765 SNES snes; 766 KSP ksp; 767 PC pc; 768 Vec coordsLocal, r, ref, real; 769 Mat J; 770 const PetscScalar *coords; 771 PetscScalar *a; 772 PetscScalar *x = NULL, *vertices = NULL; 773 PetscScalar *xi; 774 PetscReal xir[3]; 775 PetscInt coordSize, xSize; 776 777 PetscFunctionBegin; 778 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 779 PetscCall(DMGetCoordinateDM(dm, &dmCoord)); 780 PetscCall(SNESCreate(PETSC_COMM_SELF, &snes)); 781 PetscCall(SNESSetOptionsPrefix(snes, "hex_interp_")); 782 PetscCall(VecCreate(PETSC_COMM_SELF, &r)); 783 PetscCall(VecSetSizes(r, 3, 3)); 784 PetscCall(VecSetType(r, dm->vectype)); 785 PetscCall(VecDuplicate(r, &ref)); 786 PetscCall(VecDuplicate(r, &real)); 787 PetscCall(MatCreate(PETSC_COMM_SELF, &J)); 788 PetscCall(MatSetSizes(J, 3, 3, 3, 3)); 789 PetscCall(MatSetType(J, MATSEQDENSE)); 790 PetscCall(MatSetUp(J)); 791 PetscCall(SNESSetFunction(snes, r, HexMap_Private, NULL)); 792 PetscCall(SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL)); 793 PetscCall(SNESGetKSP(snes, &ksp)); 794 PetscCall(KSPGetPC(ksp, &pc)); 795 PetscCall(PCSetType(pc, PCLU)); 796 PetscCall(SNESSetFromOptions(snes)); 797 798 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 799 PetscCall(VecGetArray(v, &a)); 800 PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 801 PetscCheck(8 * 3 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid coordinate closure size %" PetscInt_FMT " should be %d", coordSize, 8 * 3); 802 PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x)); 803 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); 804 PetscCall(SNESSetFunction(snes, NULL, NULL, vertices)); 805 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices)); 806 PetscCall(VecGetArray(real, &xi)); 807 xi[0] = coords[p * ctx->dim + 0]; 808 xi[1] = coords[p * ctx->dim + 1]; 809 xi[2] = coords[p * ctx->dim + 2]; 810 PetscCall(VecRestoreArray(real, &xi)); 811 PetscCall(SNESSolve(snes, real, ref)); 812 PetscCall(VecGetArray(ref, &xi)); 813 xir[0] = PetscRealPart(xi[0]); 814 xir[1] = PetscRealPart(xi[1]); 815 xir[2] = PetscRealPart(xi[2]); 816 if (8 * ctx->dof == xSize) { 817 for (PetscInt comp = 0; comp < ctx->dof; ++comp) { 818 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]) + 819 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]; 820 } 821 } else { 822 for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp]; 823 } 824 PetscCall(VecRestoreArray(ref, &xi)); 825 PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices)); 826 PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x)); 827 PetscCall(VecRestoreArray(v, &a)); 828 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 829 830 PetscCall(SNESDestroy(&snes)); 831 PetscCall(VecDestroy(&r)); 832 PetscCall(VecDestroy(&ref)); 833 PetscCall(VecDestroy(&real)); 834 PetscCall(MatDestroy(&J)); 835 PetscFunctionReturn(PETSC_SUCCESS); 836 } 837 838 /*@C 839 DMInterpolationEvaluate - Using the input from `dm` and `x`, calculates interpolated field values at the interpolation points. 840 841 Input Parameters: 842 + ctx - The `DMInterpolationInfo` context obtained with `DMInterpolationCreate()` 843 . dm - The `DM` 844 - x - The local vector containing the field to be interpolated, can be created with `DMCreateGlobalVector()` 845 846 Output Parameter: 847 . v - The vector containing the interpolated values, obtained with `DMInterpolationGetVector()` 848 849 Level: beginner 850 851 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`, `DMInterpolationGetCoordinates()` 852 @*/ 853 PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v) 854 { 855 PetscDS ds; 856 PetscInt n, p, Nf, field; 857 PetscBool useDS = PETSC_FALSE; 858 859 PetscFunctionBegin; 860 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 861 PetscValidHeaderSpecific(x, VEC_CLASSID, 3); 862 PetscValidHeaderSpecific(v, VEC_CLASSID, 4); 863 PetscCall(VecGetLocalSize(v, &n)); 864 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); 865 if (!n) PetscFunctionReturn(PETSC_SUCCESS); 866 PetscCall(DMGetDS(dm, &ds)); 867 if (ds) { 868 useDS = PETSC_TRUE; 869 PetscCall(PetscDSGetNumFields(ds, &Nf)); 870 for (field = 0; field < Nf; ++field) { 871 PetscObject obj; 872 PetscClassId id; 873 874 PetscCall(PetscDSGetDiscretization(ds, field, &obj)); 875 PetscCall(PetscObjectGetClassId(obj, &id)); 876 if (id != PETSCFE_CLASSID && id != PETSCFV_CLASSID) { 877 useDS = PETSC_FALSE; 878 break; 879 } 880 } 881 } 882 if (useDS) { 883 const PetscScalar *coords; 884 PetscScalar *interpolant; 885 PetscInt cdim, d; 886 887 PetscCall(DMGetCoordinateDim(dm, &cdim)); 888 PetscCall(VecGetArrayRead(ctx->coords, &coords)); 889 PetscCall(VecGetArrayWrite(v, &interpolant)); 890 for (p = 0; p < ctx->n; ++p) { 891 PetscReal pcoords[3], xi[3]; 892 PetscScalar *xa = NULL; 893 PetscInt coff = 0, foff = 0, clSize; 894 895 if (ctx->cells[p] < 0) continue; 896 for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]); 897 PetscCall(DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi)); 898 PetscCall(DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa)); 899 for (field = 0; field < Nf; ++field) { 900 PetscTabulation T; 901 PetscObject obj; 902 PetscClassId id; 903 904 PetscCall(PetscDSGetDiscretization(ds, field, &obj)); 905 PetscCall(PetscObjectGetClassId(obj, &id)); 906 if (id == PETSCFE_CLASSID) { 907 PetscFE fe = (PetscFE)obj; 908 909 PetscCall(PetscFECreateTabulation(fe, 1, 1, xi, 0, &T)); 910 { 911 const PetscReal *basis = T->T[0]; 912 const PetscInt Nb = T->Nb; 913 const PetscInt Nc = T->Nc; 914 915 for (PetscInt fc = 0; fc < Nc; ++fc) { 916 interpolant[p * ctx->dof + coff + fc] = 0.0; 917 for (PetscInt f = 0; f < Nb; ++f) interpolant[p * ctx->dof + coff + fc] += xa[foff + f] * basis[(0 * Nb + f) * Nc + fc]; 918 } 919 coff += Nc; 920 foff += Nb; 921 } 922 PetscCall(PetscTabulationDestroy(&T)); 923 } else if (id == PETSCFV_CLASSID) { 924 PetscFV fv = (PetscFV)obj; 925 PetscInt Nc; 926 927 // TODO Could use reconstruction if available 928 PetscCall(PetscFVGetNumComponents(fv, &Nc)); 929 for (PetscInt fc = 0; fc < Nc; ++fc) interpolant[p * ctx->dof + coff + fc] = xa[foff + fc]; 930 coff += Nc; 931 foff += Nc; 932 } 933 } 934 PetscCall(DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa)); 935 PetscCheck(coff == ctx->dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total components %" PetscInt_FMT " != %" PetscInt_FMT " dof specified for interpolation", coff, ctx->dof); 936 PetscCheck(foff == clSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total FE/FV space size %" PetscInt_FMT " != %" PetscInt_FMT " closure size", foff, clSize); 937 } 938 PetscCall(VecRestoreArrayRead(ctx->coords, &coords)); 939 PetscCall(VecRestoreArrayWrite(v, &interpolant)); 940 } else { 941 for (PetscInt p = 0; p < ctx->n; ++p) { 942 const PetscInt cell = ctx->cells[p]; 943 DMPolytopeType ct; 944 945 PetscCall(DMPlexGetCellType(dm, cell, &ct)); 946 switch (ct) { 947 case DM_POLYTOPE_SEGMENT: 948 PetscCall(DMInterpolate_Segment_Private(ctx, dm, p, x, v)); 949 break; 950 case DM_POLYTOPE_TRIANGLE: 951 PetscCall(DMInterpolate_Triangle_Private(ctx, dm, p, x, v)); 952 break; 953 case DM_POLYTOPE_QUADRILATERAL: 954 PetscCall(DMInterpolate_Quad_Private(ctx, dm, p, x, v)); 955 break; 956 case DM_POLYTOPE_TETRAHEDRON: 957 PetscCall(DMInterpolate_Tetrahedron_Private(ctx, dm, p, x, v)); 958 break; 959 case DM_POLYTOPE_HEXAHEDRON: 960 PetscCall(DMInterpolate_Hex_Private(ctx, dm, cell, x, v)); 961 break; 962 default: 963 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for cell type %s", DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]); 964 } 965 } 966 } 967 PetscFunctionReturn(PETSC_SUCCESS); 968 } 969 970 /*@C 971 DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context 972 973 Collective 974 975 Input Parameter: 976 . ctx - the context 977 978 Level: beginner 979 980 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()` 981 @*/ 982 PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx) 983 { 984 PetscFunctionBegin; 985 PetscAssertPointer(ctx, 1); 986 PetscCall(VecDestroy(&(*ctx)->coords)); 987 PetscCall(PetscFree((*ctx)->points)); 988 PetscCall(PetscFree((*ctx)->cells)); 989 PetscCall(PetscFree(*ctx)); 990 *ctx = NULL; 991 PetscFunctionReturn(PETSC_SUCCESS); 992 } 993