1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 3 #include <petscblaslapack.h> 4 5 static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection) 6 { 7 const PetscReal p0_x = segmentA[0*2+0]; 8 const PetscReal p0_y = segmentA[0*2+1]; 9 const PetscReal p1_x = segmentA[1*2+0]; 10 const PetscReal p1_y = segmentA[1*2+1]; 11 const PetscReal p2_x = segmentB[0*2+0]; 12 const PetscReal p2_y = segmentB[0*2+1]; 13 const PetscReal p3_x = segmentB[1*2+0]; 14 const PetscReal p3_y = segmentB[1*2+1]; 15 const PetscReal s1_x = p1_x - p0_x; 16 const PetscReal s1_y = p1_y - p0_y; 17 const PetscReal s2_x = p3_x - p2_x; 18 const PetscReal s2_y = p3_y - p2_y; 19 const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y); 20 21 PetscFunctionBegin; 22 *hasIntersection = PETSC_FALSE; 23 /* Non-parallel lines */ 24 if (denom != 0.0) { 25 const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom; 26 const PetscReal t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom; 27 28 if (s >= 0 && s <= 1 && t >= 0 && t <= 1) { 29 *hasIntersection = PETSC_TRUE; 30 if (intersection) { 31 intersection[0] = p0_x + (t * s1_x); 32 intersection[1] = p0_y + (t * s1_y); 33 } 34 } 35 } 36 PetscFunctionReturn(0); 37 } 38 39 static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 40 { 41 const PetscInt embedDim = 2; 42 const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON; 43 PetscReal x = PetscRealPart(point[0]); 44 PetscReal y = PetscRealPart(point[1]); 45 PetscReal v0[2], J[4], invJ[4], detJ; 46 PetscReal xi, eta; 47 PetscErrorCode ierr; 48 49 PetscFunctionBegin; 50 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 51 xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]); 52 eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]); 53 54 if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0+eps)) *cell = c; 55 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 56 PetscFunctionReturn(0); 57 } 58 59 static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[]) 60 { 61 const PetscInt embedDim = 2; 62 PetscReal x = PetscRealPart(point[0]); 63 PetscReal y = PetscRealPart(point[1]); 64 PetscReal v0[2], J[4], invJ[4], detJ; 65 PetscReal xi, eta, r; 66 PetscErrorCode ierr; 67 68 PetscFunctionBegin; 69 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 70 xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]); 71 eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]); 72 73 xi = PetscMax(xi, 0.0); 74 eta = PetscMax(eta, 0.0); 75 r = (xi + eta)/2.0; 76 if (xi + eta > 2.0) { 77 r = (xi + eta)/2.0; 78 xi /= r; 79 eta /= r; 80 } 81 cpoint[0] = J[0*embedDim+0]*xi + J[0*embedDim+1]*eta + v0[0]; 82 cpoint[1] = J[1*embedDim+0]*xi + J[1*embedDim+1]*eta + v0[1]; 83 PetscFunctionReturn(0); 84 } 85 86 static PetscErrorCode DMPlexLocatePoint_General_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 87 { 88 PetscSection coordSection; 89 Vec coordsLocal; 90 PetscScalar *coords = NULL; 91 const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0}; 92 PetscReal x = PetscRealPart(point[0]); 93 PetscReal y = PetscRealPart(point[1]); 94 PetscInt crossings = 0, f; 95 PetscErrorCode ierr; 96 97 PetscFunctionBegin; 98 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 99 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 100 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 101 for (f = 0; f < 4; ++f) { 102 PetscReal x_i = PetscRealPart(coords[faces[2*f+0]*2+0]); 103 PetscReal y_i = PetscRealPart(coords[faces[2*f+0]*2+1]); 104 PetscReal x_j = PetscRealPart(coords[faces[2*f+1]*2+0]); 105 PetscReal y_j = PetscRealPart(coords[faces[2*f+1]*2+1]); 106 PetscReal slope = (y_j - y_i) / (x_j - x_i); 107 PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE; 108 PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE; 109 PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE; 110 if ((cond1 || cond2) && above) ++crossings; 111 } 112 if (crossings % 2) *cell = c; 113 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 114 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 115 PetscFunctionReturn(0); 116 } 117 118 static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 119 { 120 const PetscInt embedDim = 3; 121 PetscReal v0[3], J[9], invJ[9], detJ; 122 PetscReal x = PetscRealPart(point[0]); 123 PetscReal y = PetscRealPart(point[1]); 124 PetscReal z = PetscRealPart(point[2]); 125 PetscReal xi, eta, zeta; 126 PetscErrorCode ierr; 127 128 PetscFunctionBegin; 129 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 130 xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]); 131 eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]); 132 zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]); 133 134 if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c; 135 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 136 PetscFunctionReturn(0); 137 } 138 139 static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 140 { 141 PetscSection coordSection; 142 Vec coordsLocal; 143 PetscScalar *coords; 144 const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5, 145 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4}; 146 PetscBool found = PETSC_TRUE; 147 PetscInt f; 148 PetscErrorCode ierr; 149 150 PetscFunctionBegin; 151 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 152 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 153 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 154 for (f = 0; f < 6; ++f) { 155 /* Check the point is under plane */ 156 /* Get face normal */ 157 PetscReal v_i[3]; 158 PetscReal v_j[3]; 159 PetscReal normal[3]; 160 PetscReal pp[3]; 161 PetscReal dot; 162 163 v_i[0] = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]); 164 v_i[1] = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]); 165 v_i[2] = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]); 166 v_j[0] = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]); 167 v_j[1] = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]); 168 v_j[2] = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]); 169 normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1]; 170 normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2]; 171 normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0]; 172 pp[0] = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]); 173 pp[1] = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]); 174 pp[2] = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]); 175 dot = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2]; 176 177 /* Check that projected point is in face (2D location problem) */ 178 if (dot < 0.0) { 179 found = PETSC_FALSE; 180 break; 181 } 182 } 183 if (found) *cell = c; 184 else *cell = DMLOCATEPOINT_POINT_NOT_FOUND; 185 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 186 PetscFunctionReturn(0); 187 } 188 189 static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[]) 190 { 191 PetscInt d; 192 193 PetscFunctionBegin; 194 box->dim = dim; 195 for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]); 196 PetscFunctionReturn(0); 197 } 198 199 PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box) 200 { 201 PetscErrorCode ierr; 202 203 PetscFunctionBegin; 204 ierr = PetscMalloc1(1, box);CHKERRQ(ierr); 205 ierr = PetscGridHashInitialize_Internal(*box, dim, point);CHKERRQ(ierr); 206 PetscFunctionReturn(0); 207 } 208 209 PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[]) 210 { 211 PetscInt d; 212 213 PetscFunctionBegin; 214 for (d = 0; d < box->dim; ++d) { 215 box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d])); 216 box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d])); 217 } 218 PetscFunctionReturn(0); 219 } 220 221 /* 222 PetscGridHashSetGrid - Divide the grid into boxes 223 224 Not collective 225 226 Input Parameters: 227 + box - The grid hash object 228 . n - The number of boxes in each dimension, or PETSC_DETERMINE 229 - h - The box size in each dimension, only used if n[d] == PETSC_DETERMINE 230 231 Level: developer 232 233 .seealso: PetscGridHashCreate() 234 */ 235 PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[]) 236 { 237 PetscInt d; 238 239 PetscFunctionBegin; 240 for (d = 0; d < box->dim; ++d) { 241 box->extent[d] = box->upper[d] - box->lower[d]; 242 if (n[d] == PETSC_DETERMINE) { 243 box->h[d] = h[d]; 244 box->n[d] = PetscCeilReal(box->extent[d]/h[d]); 245 } else { 246 box->n[d] = n[d]; 247 box->h[d] = box->extent[d]/n[d]; 248 } 249 } 250 PetscFunctionReturn(0); 251 } 252 253 /* 254 PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point 255 256 Not collective 257 258 Input Parameters: 259 + box - The grid hash object 260 . numPoints - The number of input points 261 - points - The input point coordinates 262 263 Output Parameters: 264 + dboxes - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim) 265 - boxes - An array of numPoints integers expressing the enclosing box as single number, or NULL 266 267 Level: developer 268 269 .seealso: PetscGridHashCreate() 270 */ 271 PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[]) 272 { 273 const PetscReal *lower = box->lower; 274 const PetscReal *upper = box->upper; 275 const PetscReal *h = box->h; 276 const PetscInt *n = box->n; 277 const PetscInt dim = box->dim; 278 PetscInt d, p; 279 280 PetscFunctionBegin; 281 for (p = 0; p < numPoints; ++p) { 282 for (d = 0; d < dim; ++d) { 283 PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]); 284 285 if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1; 286 if (dbox < 0 || dbox >= n[d]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input point %d (%g, %g, %g) is outside of our bounding box", 287 p, PetscRealPart(points[p*dim+0]), dim > 1 ? PetscRealPart(points[p*dim+1]) : 0.0, dim > 2 ? PetscRealPart(points[p*dim+2]) : 0.0); 288 dboxes[p*dim+d] = dbox; 289 } 290 if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1]; 291 } 292 PetscFunctionReturn(0); 293 } 294 295 PetscErrorCode PetscGridHashDestroy(PetscGridHash *box) 296 { 297 PetscErrorCode ierr; 298 299 PetscFunctionBegin; 300 if (*box) { 301 ierr = PetscSectionDestroy(&(*box)->cellSection);CHKERRQ(ierr); 302 ierr = ISDestroy(&(*box)->cells);CHKERRQ(ierr); 303 ierr = DMLabelDestroy(&(*box)->cellsSparse);CHKERRQ(ierr); 304 } 305 ierr = PetscFree(*box);CHKERRQ(ierr); 306 PetscFunctionReturn(0); 307 } 308 309 PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell) 310 { 311 PetscInt coneSize; 312 PetscErrorCode ierr; 313 314 PetscFunctionBegin; 315 switch (dim) { 316 case 2: 317 ierr = DMPlexGetConeSize(dm, cellStart, &coneSize);CHKERRQ(ierr); 318 switch (coneSize) { 319 case 3: 320 ierr = DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);CHKERRQ(ierr); 321 break; 322 case 4: 323 ierr = DMPlexLocatePoint_General_2D_Internal(dm, point, cellStart, cell);CHKERRQ(ierr); 324 break; 325 default: 326 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize); 327 } 328 break; 329 case 3: 330 ierr = DMPlexGetConeSize(dm, cellStart, &coneSize);CHKERRQ(ierr); 331 switch (coneSize) { 332 case 4: 333 ierr = DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);CHKERRQ(ierr); 334 break; 335 case 6: 336 ierr = DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);CHKERRQ(ierr); 337 break; 338 default: 339 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize); 340 } 341 break; 342 default: 343 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for mesh dimension %D", dim); 344 } 345 PetscFunctionReturn(0); 346 } 347 348 /* 349 DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point 350 */ 351 PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[]) 352 { 353 PetscInt coneSize; 354 PetscErrorCode ierr; 355 356 PetscFunctionBegin; 357 switch (dim) { 358 case 2: 359 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 360 switch (coneSize) { 361 case 3: 362 ierr = DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint);CHKERRQ(ierr); 363 break; 364 #if 0 365 case 4: 366 ierr = DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint);CHKERRQ(ierr); 367 break; 368 #endif 369 default: 370 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell with cone size %D", coneSize); 371 } 372 break; 373 #if 0 374 case 3: 375 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 376 switch (coneSize) { 377 case 4: 378 ierr = DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint);CHKERRQ(ierr); 379 break; 380 case 6: 381 ierr = DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint);CHKERRQ(ierr); 382 break; 383 default: 384 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell with cone size %D", coneSize); 385 } 386 break; 387 #endif 388 default: 389 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for mesh dimension %D", dim); 390 } 391 PetscFunctionReturn(0); 392 } 393 394 /* 395 DMPlexComputeGridHash_Internal - Create a grid hash structure covering the Plex 396 397 Collective on DM 398 399 Input Parameter: 400 . dm - The Plex 401 402 Output Parameter: 403 . localBox - The grid hash object 404 405 Level: developer 406 407 .seealso: PetscGridHashCreate(), PetscGridHashGetEnclosingBox() 408 */ 409 PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox) 410 { 411 MPI_Comm comm; 412 PetscGridHash lbox; 413 Vec coordinates; 414 PetscSection coordSection; 415 Vec coordsLocal; 416 const PetscScalar *coords; 417 PetscInt *dboxes, *boxes; 418 PetscInt n[3] = {10, 10, 10}; 419 PetscInt dim, N, cStart, cEnd, cMax, c, i; 420 PetscErrorCode ierr; 421 422 PetscFunctionBegin; 423 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 424 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 425 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 426 if (dim != 2) SETERRQ(comm, PETSC_ERR_SUP, "I have only coded this for 2D"); 427 ierr = VecGetLocalSize(coordinates, &N);CHKERRQ(ierr); 428 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 429 ierr = PetscGridHashCreate(comm, dim, coords, &lbox);CHKERRQ(ierr); 430 for (i = 0; i < N; i += dim) {ierr = PetscGridHashEnlarge(lbox, &coords[i]);CHKERRQ(ierr);} 431 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 432 ierr = PetscGridHashSetGrid(lbox, n, NULL);CHKERRQ(ierr); 433 #if 0 434 /* Could define a custom reduction to merge these */ 435 ierr = MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);CHKERRQ(ierr); 436 ierr = MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);CHKERRQ(ierr); 437 #endif 438 /* Is there a reason to snap the local bounding box to a division of the global box? */ 439 /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */ 440 /* Create label */ 441 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 442 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 443 if (cMax >= 0) cEnd = PetscMin(cEnd, cMax); 444 ierr = DMLabelCreate("cells", &lbox->cellsSparse);CHKERRQ(ierr); 445 ierr = DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);CHKERRQ(ierr); 446 /* Compute boxes which overlap each cell: http://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */ 447 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 448 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 449 ierr = PetscCalloc2(16 * dim, &dboxes, 16, &boxes);CHKERRQ(ierr); 450 for (c = cStart; c < cEnd; ++c) { 451 const PetscReal *h = lbox->h; 452 PetscScalar *ccoords = NULL; 453 PetscInt csize = 0; 454 PetscScalar point[3]; 455 PetscInt dlim[6], d, e, i, j, k; 456 457 /* Find boxes enclosing each vertex */ 458 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);CHKERRQ(ierr); 459 ierr = PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);CHKERRQ(ierr); 460 /* Mark cells containing the vertices */ 461 for (e = 0; e < csize/dim; ++e) {ierr = DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);CHKERRQ(ierr);} 462 /* Get grid of boxes containing these */ 463 for (d = 0; d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];} 464 for (d = dim; d < 3; ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;} 465 for (e = 1; e < dim+1; ++e) { 466 for (d = 0; d < dim; ++d) { 467 dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]); 468 dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]); 469 } 470 } 471 /* Check for intersection of box with cell */ 472 for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) { 473 for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) { 474 for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) { 475 const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i; 476 PetscScalar cpoint[3]; 477 PetscInt cell, edge, ii, jj, kk; 478 479 /* Check whether cell contains any vertex of these subboxes TODO vectorize this */ 480 for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) { 481 for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) { 482 for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) { 483 484 ierr = DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);CHKERRQ(ierr); 485 if (cell >= 0) {DMLabelSetValue(lbox->cellsSparse, c, box);CHKERRQ(ierr); ii = jj = kk = 2;} 486 } 487 } 488 } 489 /* Check whether cell edge intersects any edge of these subboxes TODO vectorize this */ 490 for (edge = 0; edge < dim+1; ++edge) { 491 PetscReal segA[6], segB[6]; 492 493 for (d = 0; d < dim; ++d) {segA[d] = PetscRealPart(ccoords[edge*dim+d]); segA[dim+d] = PetscRealPart(ccoords[((edge+1)%(dim+1))*dim+d]);} 494 for (kk = 0; kk < (dim > 2 ? 2 : 1); ++kk) { 495 if (dim > 2) {segB[2] = PetscRealPart(point[2]); 496 segB[dim+2] = PetscRealPart(point[2]) + kk*h[2];} 497 for (jj = 0; jj < (dim > 1 ? 2 : 1); ++jj) { 498 if (dim > 1) {segB[1] = PetscRealPart(point[1]); 499 segB[dim+1] = PetscRealPart(point[1]) + jj*h[1];} 500 for (ii = 0; ii < 2; ++ii) { 501 PetscBool intersects; 502 503 segB[0] = PetscRealPart(point[0]); 504 segB[dim+0] = PetscRealPart(point[0]) + ii*h[0]; 505 ierr = DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);CHKERRQ(ierr); 506 if (intersects) {DMLabelSetValue(lbox->cellsSparse, c, box);CHKERRQ(ierr); edge = ii = jj = kk = dim+1;} 507 } 508 } 509 } 510 } 511 } 512 } 513 } 514 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);CHKERRQ(ierr); 515 } 516 ierr = PetscFree2(dboxes, boxes);CHKERRQ(ierr); 517 ierr = DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);CHKERRQ(ierr); 518 ierr = DMLabelDestroy(&lbox->cellsSparse);CHKERRQ(ierr); 519 *localBox = lbox; 520 PetscFunctionReturn(0); 521 } 522 523 PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF) 524 { 525 DM_Plex *mesh = (DM_Plex *) dm->data; 526 PetscBool hash = mesh->useHashLocation; 527 PetscInt bs, numPoints, p, numFound, *found = NULL; 528 PetscInt dim, cStart, cEnd, cMax, numCells, c; 529 const PetscInt *boxCells; 530 PetscSFNode *cells; 531 PetscScalar *a; 532 PetscMPIInt result; 533 PetscErrorCode ierr; 534 535 PetscFunctionBegin; 536 if (ltype == DM_POINTLOCATION_NEAREST && !hash) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Nearest point location only supported with grid hashing. Use -dm_plex_hash_location to enable it."); 537 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 538 ierr = VecGetBlockSize(v, &bs);CHKERRQ(ierr); 539 ierr = MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);CHKERRQ(ierr); 540 if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported"); 541 if (bs != dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Block size for point vector %D must be the mesh coordinate dimension %D", bs, dim); 542 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 543 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 544 if (cMax >= 0) cEnd = PetscMin(cEnd, cMax); 545 ierr = VecGetLocalSize(v, &numPoints);CHKERRQ(ierr); 546 ierr = VecGetArray(v, &a);CHKERRQ(ierr); 547 numPoints /= bs; 548 ierr = PetscMalloc1(numPoints, &cells);CHKERRQ(ierr); 549 if (hash) { 550 if (!mesh->lbox) {ierr = PetscInfo(dm, "Initializing grid hashing");CHKERRQ(ierr);ierr = DMPlexComputeGridHash_Internal(dm, &mesh->lbox);CHKERRQ(ierr);} 551 /* Designate the local box for each point */ 552 /* Send points to correct process */ 553 /* Search cells that lie in each subbox */ 554 /* Should we bin points before doing search? */ 555 ierr = ISGetIndices(mesh->lbox->cells, &boxCells);CHKERRQ(ierr); 556 } 557 for (p = 0, numFound = 0; p < numPoints; ++p) { 558 const PetscScalar *point = &a[p*bs]; 559 PetscInt dbin[3], bin, cell = -1, cellOffset; 560 561 cells[p].rank = 0; 562 cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND; 563 if (hash) { 564 ierr = PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);CHKERRQ(ierr); 565 /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */ 566 ierr = PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);CHKERRQ(ierr); 567 ierr = PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);CHKERRQ(ierr); 568 for (c = cellOffset; c < cellOffset + numCells; ++c) { 569 ierr = DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell);CHKERRQ(ierr); 570 if (cell >= 0) { 571 cells[p].rank = 0; 572 cells[p].index = cell; 573 numFound++; 574 break; 575 } 576 } 577 } else { 578 for (c = cStart; c < cEnd; ++c) { 579 ierr = DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);CHKERRQ(ierr); 580 if (cell >= 0) { 581 cells[p].rank = 0; 582 cells[p].index = cell; 583 numFound++; 584 break; 585 } 586 } 587 } 588 } 589 if (hash) {ierr = ISRestoreIndices(mesh->lbox->cells, &boxCells);CHKERRQ(ierr);} 590 if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) { 591 for (p = 0; p < numPoints; p++) { 592 const PetscScalar *point = &a[p*bs]; 593 PetscReal cpoint[3], diff[3], dist, distMax = PETSC_MAX_REAL; 594 PetscInt dbin[3], bin, cellOffset, d; 595 596 if (cells[p].index < 0) { 597 ++numFound; 598 ierr = PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);CHKERRQ(ierr); 599 ierr = PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);CHKERRQ(ierr); 600 ierr = PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);CHKERRQ(ierr); 601 for (c = cellOffset; c < cellOffset + numCells; ++c) { 602 ierr = DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint);CHKERRQ(ierr); 603 for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]); 604 dist = DMPlex_NormD_Internal(dim, diff); 605 if (dist < distMax) { 606 for (d = 0; d < dim; ++d) a[p*bs+d] = cpoint[d]; 607 cells[p].rank = 0; 608 cells[p].index = boxCells[c]; 609 distMax = dist; 610 break; 611 } 612 } 613 } 614 } 615 } 616 /* This code is only be relevant when interfaced to parallel point location */ 617 /* Check for highest numbered proc that claims a point (do we care?) */ 618 if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) { 619 ierr = PetscMalloc1(numFound,&found);CHKERRQ(ierr); 620 for (p = 0, numFound = 0; p < numPoints; p++) { 621 if (cells[p].rank >= 0 && cells[p].index >= 0) { 622 if (numFound < p) { 623 cells[numFound] = cells[p]; 624 } 625 found[numFound++] = p; 626 } 627 } 628 } 629 ierr = VecRestoreArray(v, &a);CHKERRQ(ierr); 630 ierr = PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);CHKERRQ(ierr); 631 PetscFunctionReturn(0); 632 } 633 634 /*@C 635 DMPlexComputeProjection2Dto1D - Rewrite coordinates to be the 1D projection of the 2D coordinates 636 637 Not collective 638 639 Input Parameter: 640 . coords - The coordinates of a segment 641 642 Output Parameters: 643 + coords - The new y-coordinate, and 0 for x 644 - R - The rotation which accomplishes the projection 645 646 Level: developer 647 648 .seealso: DMPlexComputeProjection3Dto1D(), DMPlexComputeProjection3Dto2D() 649 @*/ 650 PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[]) 651 { 652 const PetscReal x = PetscRealPart(coords[2] - coords[0]); 653 const PetscReal y = PetscRealPart(coords[3] - coords[1]); 654 const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r; 655 656 PetscFunctionBegin; 657 R[0] = c; R[1] = -s; 658 R[2] = s; R[3] = c; 659 coords[0] = 0.0; 660 coords[1] = r; 661 PetscFunctionReturn(0); 662 } 663 664 /*@C 665 DMPlexComputeProjection3Dto1D - Rewrite coordinates to be the 1D projection of the 3D coordinates 666 667 Not collective 668 669 Input Parameter: 670 . coords - The coordinates of a segment 671 672 Output Parameters: 673 + coords - The new y-coordinate, and 0 for x and z 674 - R - The rotation which accomplishes the projection 675 676 Note: This uses the basis completion described by Frisvad in http://www.imm.dtu.dk/~jerf/papers/abstracts/onb.html, DOI:10.1080/2165347X.2012.689606 677 678 Level: developer 679 680 .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto2D() 681 @*/ 682 PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[]) 683 { 684 PetscReal x = PetscRealPart(coords[3] - coords[0]); 685 PetscReal y = PetscRealPart(coords[4] - coords[1]); 686 PetscReal z = PetscRealPart(coords[5] - coords[2]); 687 PetscReal r = PetscSqrtReal(x*x + y*y + z*z); 688 PetscReal rinv = 1. / r; 689 PetscFunctionBegin; 690 691 x *= rinv; y *= rinv; z *= rinv; 692 if (x > 0.) { 693 PetscReal inv1pX = 1./ (1. + x); 694 695 R[0] = x; R[1] = -y; R[2] = -z; 696 R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] = -y*z*inv1pX; 697 R[6] = z; R[7] = -y*z*inv1pX; R[8] = 1. - z*z*inv1pX; 698 } 699 else { 700 PetscReal inv1mX = 1./ (1. - x); 701 702 R[0] = x; R[1] = z; R[2] = y; 703 R[3] = y; R[4] = -y*z*inv1mX; R[5] = 1. - y*y*inv1mX; 704 R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] = -y*z*inv1mX; 705 } 706 coords[0] = 0.0; 707 coords[1] = r; 708 PetscFunctionReturn(0); 709 } 710 711 /*@ 712 DMPlexComputeProjection3Dto2D - Rewrite coordinates to be the 2D projection of the 3D coordinates 713 714 Not collective 715 716 Input Parameter: 717 . coords - The coordinates of a segment 718 719 Output Parameters: 720 + coords - The new y- and z-coordinates, and 0 for x 721 - R - The rotation which accomplishes the projection 722 723 Level: developer 724 725 .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto1D() 726 @*/ 727 PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[]) 728 { 729 PetscReal x1[3], x2[3], n[3], norm; 730 PetscReal x1p[3], x2p[3], xnp[3]; 731 PetscReal sqrtz, alpha; 732 const PetscInt dim = 3; 733 PetscInt d, e, p; 734 735 PetscFunctionBegin; 736 /* 0) Calculate normal vector */ 737 for (d = 0; d < dim; ++d) { 738 x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]); 739 x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]); 740 } 741 n[0] = x1[1]*x2[2] - x1[2]*x2[1]; 742 n[1] = x1[2]*x2[0] - x1[0]*x2[2]; 743 n[2] = x1[0]*x2[1] - x1[1]*x2[0]; 744 norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]); 745 n[0] /= norm; 746 n[1] /= norm; 747 n[2] /= norm; 748 /* 1) Take the normal vector and rotate until it is \hat z 749 750 Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then 751 752 R = / alpha nx nz alpha ny nz -1/alpha \ 753 | -alpha ny alpha nx 0 | 754 \ nx ny nz / 755 756 will rotate the normal vector to \hat z 757 */ 758 sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]); 759 /* Check for n = z */ 760 if (sqrtz < 1.0e-10) { 761 const PetscInt s = PetscSign(n[2]); 762 /* If nz < 0, rotate 180 degrees around x-axis */ 763 for (p = 3; p < coordSize/3; ++p) { 764 coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]); 765 coords[p*2+1] = (PetscRealPart(coords[p*dim+1] - coords[0*dim+1])) * s; 766 } 767 coords[0] = 0.0; 768 coords[1] = 0.0; 769 coords[2] = x1[0]; 770 coords[3] = x1[1] * s; 771 coords[4] = x2[0]; 772 coords[5] = x2[1] * s; 773 R[0] = 1.0; R[1] = 0.0; R[2] = 0.0; 774 R[3] = 0.0; R[4] = 1.0 * s; R[5] = 0.0; 775 R[6] = 0.0; R[7] = 0.0; R[8] = 1.0 * s; 776 PetscFunctionReturn(0); 777 } 778 alpha = 1.0/sqrtz; 779 R[0] = alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz; 780 R[3] = -alpha*n[1]; R[4] = alpha*n[0]; R[5] = 0.0; 781 R[6] = n[0]; R[7] = n[1]; R[8] = n[2]; 782 for (d = 0; d < dim; ++d) { 783 x1p[d] = 0.0; 784 x2p[d] = 0.0; 785 for (e = 0; e < dim; ++e) { 786 x1p[d] += R[d*dim+e]*x1[e]; 787 x2p[d] += R[d*dim+e]*x2[e]; 788 } 789 } 790 if (PetscAbsReal(x1p[2]) > 10. * PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 791 if (PetscAbsReal(x2p[2]) > 10. * PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 792 /* 2) Project to (x, y) */ 793 for (p = 3; p < coordSize/3; ++p) { 794 for (d = 0; d < dim; ++d) { 795 xnp[d] = 0.0; 796 for (e = 0; e < dim; ++e) { 797 xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]); 798 } 799 if (d < dim-1) coords[p*2+d] = xnp[d]; 800 } 801 } 802 coords[0] = 0.0; 803 coords[1] = 0.0; 804 coords[2] = x1p[0]; 805 coords[3] = x1p[1]; 806 coords[4] = x2p[0]; 807 coords[5] = x2p[1]; 808 /* Output R^T which rotates \hat z to the input normal */ 809 for (d = 0; d < dim; ++d) { 810 for (e = d+1; e < dim; ++e) { 811 PetscReal tmp; 812 813 tmp = R[d*dim+e]; 814 R[d*dim+e] = R[e*dim+d]; 815 R[e*dim+d] = tmp; 816 } 817 } 818 PetscFunctionReturn(0); 819 } 820 821 PETSC_UNUSED 822 PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[]) 823 { 824 /* Signed volume is 1/2 the determinant 825 826 | 1 1 1 | 827 | x0 x1 x2 | 828 | y0 y1 y2 | 829 830 but if x0,y0 is the origin, we have 831 832 | x1 x2 | 833 | y1 y2 | 834 */ 835 const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1]; 836 const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1]; 837 PetscReal M[4], detM; 838 M[0] = x1; M[1] = x2; 839 M[2] = y1; M[3] = y2; 840 DMPlex_Det2D_Internal(&detM, M); 841 *vol = 0.5*detM; 842 (void)PetscLogFlops(5.0); 843 } 844 845 PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[]) 846 { 847 DMPlex_Det2D_Internal(vol, coords); 848 *vol *= 0.5; 849 } 850 851 PETSC_UNUSED 852 PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[]) 853 { 854 /* Signed volume is 1/6th of the determinant 855 856 | 1 1 1 1 | 857 | x0 x1 x2 x3 | 858 | y0 y1 y2 y3 | 859 | z0 z1 z2 z3 | 860 861 but if x0,y0,z0 is the origin, we have 862 863 | x1 x2 x3 | 864 | y1 y2 y3 | 865 | z1 z2 z3 | 866 */ 867 const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2]; 868 const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2]; 869 const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2]; 870 PetscReal M[9], detM; 871 M[0] = x1; M[1] = x2; M[2] = x3; 872 M[3] = y1; M[4] = y2; M[5] = y3; 873 M[6] = z1; M[7] = z2; M[8] = z3; 874 DMPlex_Det3D_Internal(&detM, M); 875 *vol = -0.16666666666666666666666*detM; 876 (void)PetscLogFlops(10.0); 877 } 878 879 PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[]) 880 { 881 DMPlex_Det3D_Internal(vol, coords); 882 *vol *= -0.16666666666666666666666; 883 } 884 885 static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 886 { 887 PetscSection coordSection; 888 Vec coordinates; 889 const PetscScalar *coords; 890 PetscInt dim, d, off; 891 PetscErrorCode ierr; 892 893 PetscFunctionBegin; 894 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 895 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 896 ierr = PetscSectionGetDof(coordSection,e,&dim);CHKERRQ(ierr); 897 if (!dim) PetscFunctionReturn(0); 898 ierr = PetscSectionGetOffset(coordSection,e,&off);CHKERRQ(ierr); 899 ierr = VecGetArrayRead(coordinates,&coords);CHKERRQ(ierr); 900 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);} 901 ierr = VecRestoreArrayRead(coordinates,&coords);CHKERRQ(ierr); 902 *detJ = 1.; 903 if (J) { 904 for (d = 0; d < dim * dim; d++) J[d] = 0.; 905 for (d = 0; d < dim; d++) J[d * dim + d] = 1.; 906 if (invJ) { 907 for (d = 0; d < dim * dim; d++) invJ[d] = 0.; 908 for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.; 909 } 910 } 911 PetscFunctionReturn(0); 912 } 913 914 static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 915 { 916 PetscSection coordSection; 917 Vec coordinates; 918 PetscScalar *coords = NULL; 919 PetscInt numCoords, d, pStart, pEnd, numSelfCoords = 0; 920 PetscErrorCode ierr; 921 922 PetscFunctionBegin; 923 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 924 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 925 ierr = PetscSectionGetChart(coordSection,&pStart,&pEnd);CHKERRQ(ierr); 926 if (e >= pStart && e < pEnd) {ierr = PetscSectionGetDof(coordSection,e,&numSelfCoords);CHKERRQ(ierr);} 927 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 928 numCoords = numSelfCoords ? numSelfCoords : numCoords; 929 if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL"); 930 *detJ = 0.0; 931 if (numCoords == 6) { 932 const PetscInt dim = 3; 933 PetscReal R[9], J0; 934 935 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 936 ierr = DMPlexComputeProjection3Dto1D(coords, R);CHKERRQ(ierr); 937 if (J) { 938 J0 = 0.5*PetscRealPart(coords[1]); 939 J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2]; 940 J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5]; 941 J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8]; 942 DMPlex_Det3D_Internal(detJ, J); 943 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 944 } 945 } else if (numCoords == 4) { 946 const PetscInt dim = 2; 947 PetscReal R[4], J0; 948 949 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 950 ierr = DMPlexComputeProjection2Dto1D(coords, R);CHKERRQ(ierr); 951 if (J) { 952 J0 = 0.5*PetscRealPart(coords[1]); 953 J[0] = R[0]*J0; J[1] = R[1]; 954 J[2] = R[2]*J0; J[3] = R[3]; 955 DMPlex_Det2D_Internal(detJ, J); 956 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 957 } 958 } else if (numCoords == 2) { 959 const PetscInt dim = 1; 960 961 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 962 if (J) { 963 J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0])); 964 *detJ = J[0]; 965 ierr = PetscLogFlops(2.0);CHKERRQ(ierr); 966 if (invJ) {invJ[0] = 1.0/J[0]; ierr = PetscLogFlops(1.0);CHKERRQ(ierr);} 967 } 968 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords); 969 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 970 PetscFunctionReturn(0); 971 } 972 973 static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 974 { 975 PetscSection coordSection; 976 Vec coordinates; 977 PetscScalar *coords = NULL; 978 PetscInt numCoords, d, f, g; 979 PetscErrorCode ierr; 980 981 PetscFunctionBegin; 982 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 983 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 984 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 985 *detJ = 0.0; 986 if (numCoords == 9) { 987 const PetscInt dim = 3; 988 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 989 990 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 991 ierr = DMPlexComputeProjection3Dto2D(numCoords, coords, R);CHKERRQ(ierr); 992 if (J) { 993 const PetscInt pdim = 2; 994 995 for (d = 0; d < pdim; d++) { 996 for (f = 0; f < pdim; f++) { 997 J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 998 } 999 } 1000 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1001 DMPlex_Det3D_Internal(detJ, J0); 1002 for (d = 0; d < dim; d++) { 1003 for (f = 0; f < dim; f++) { 1004 J[d*dim+f] = 0.0; 1005 for (g = 0; g < dim; g++) { 1006 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 1007 } 1008 } 1009 } 1010 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1011 } 1012 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1013 } else if (numCoords == 6) { 1014 const PetscInt dim = 2; 1015 1016 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1017 if (J) { 1018 for (d = 0; d < dim; d++) { 1019 for (f = 0; f < dim; f++) { 1020 J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 1021 } 1022 } 1023 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1024 DMPlex_Det2D_Internal(detJ, J); 1025 } 1026 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 1027 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords); 1028 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1029 PetscFunctionReturn(0); 1030 } 1031 1032 static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1033 { 1034 PetscSection coordSection; 1035 Vec coordinates; 1036 PetscScalar *coords = NULL; 1037 PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd; 1038 PetscErrorCode ierr; 1039 1040 PetscFunctionBegin; 1041 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1042 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1043 ierr = PetscSectionGetChart(coordSection,&pStart,&pEnd);CHKERRQ(ierr); 1044 if (e >= pStart && e < pEnd) {ierr = PetscSectionGetDof(coordSection,e,&numSelfCoords);CHKERRQ(ierr);} 1045 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1046 numCoords = numSelfCoords ? numSelfCoords : numCoords; 1047 *detJ = 0.0; 1048 if (numCoords == 12) { 1049 const PetscInt dim = 3; 1050 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 1051 1052 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1053 ierr = DMPlexComputeProjection3Dto2D(numCoords, coords, R);CHKERRQ(ierr); 1054 if (J) { 1055 const PetscInt pdim = 2; 1056 1057 for (d = 0; d < pdim; d++) { 1058 J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 1059 J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 1060 } 1061 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1062 DMPlex_Det3D_Internal(detJ, J0); 1063 for (d = 0; d < dim; d++) { 1064 for (f = 0; f < dim; f++) { 1065 J[d*dim+f] = 0.0; 1066 for (g = 0; g < dim; g++) { 1067 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 1068 } 1069 } 1070 } 1071 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1072 } 1073 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1074 } else if (numCoords == 8) { 1075 const PetscInt dim = 2; 1076 1077 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1078 if (J) { 1079 for (d = 0; d < dim; d++) { 1080 J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1081 J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1082 } 1083 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1084 DMPlex_Det2D_Internal(detJ, J); 1085 } 1086 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 1087 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords); 1088 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1089 PetscFunctionReturn(0); 1090 } 1091 1092 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1093 { 1094 PetscSection coordSection; 1095 Vec coordinates; 1096 PetscScalar *coords = NULL; 1097 const PetscInt dim = 3; 1098 PetscInt d; 1099 PetscErrorCode ierr; 1100 1101 PetscFunctionBegin; 1102 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1103 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1104 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1105 *detJ = 0.0; 1106 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1107 if (J) { 1108 for (d = 0; d < dim; d++) { 1109 /* I orient with outward face normals */ 1110 J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d])); 1111 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1112 J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1113 } 1114 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1115 DMPlex_Det3D_Internal(detJ, J); 1116 } 1117 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1118 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1119 PetscFunctionReturn(0); 1120 } 1121 1122 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1123 { 1124 PetscSection coordSection; 1125 Vec coordinates; 1126 PetscScalar *coords = NULL; 1127 const PetscInt dim = 3; 1128 PetscInt d; 1129 PetscErrorCode ierr; 1130 1131 PetscFunctionBegin; 1132 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1133 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1134 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1135 *detJ = 0.0; 1136 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1137 if (J) { 1138 for (d = 0; d < dim; d++) { 1139 J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1140 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1141 J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d])); 1142 } 1143 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1144 DMPlex_Det3D_Internal(detJ, J); 1145 } 1146 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1147 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1148 PetscFunctionReturn(0); 1149 } 1150 1151 /*@C 1152 DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell 1153 1154 Collective on DM 1155 1156 Input Arguments: 1157 + dm - the DM 1158 - cell - the cell 1159 1160 Output Arguments: 1161 + v0 - the translation part of this affine transform 1162 . J - the Jacobian of the transform from the reference element 1163 . invJ - the inverse of the Jacobian 1164 - detJ - the Jacobian determinant 1165 1166 Level: advanced 1167 1168 Fortran Notes: 1169 Since it returns arrays, this routine is only available in Fortran 90, and you must 1170 include petsc.h90 in your code. 1171 1172 .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinateVec() 1173 @*/ 1174 PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1175 { 1176 PetscInt depth, dim, coneSize; 1177 DMLabel depthLabel; 1178 PetscErrorCode ierr; 1179 1180 PetscFunctionBegin; 1181 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1182 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 1183 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 1184 ierr = DMLabelGetValue(depthLabel, cell, &dim);CHKERRQ(ierr); 1185 if (depth == 1 && dim == 1) { 1186 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1187 } 1188 switch (dim) { 1189 case 0: 1190 ierr = DMPlexComputePointGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1191 break; 1192 case 1: 1193 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1194 break; 1195 case 2: 1196 switch (coneSize) { 1197 case 3: 1198 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1199 break; 1200 case 4: 1201 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1202 break; 1203 default: 1204 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 1205 } 1206 break; 1207 case 3: 1208 switch (coneSize) { 1209 case 4: 1210 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1211 break; 1212 case 6: /* Faces */ 1213 case 8: /* Vertices */ 1214 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1215 break; 1216 default: 1217 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 1218 } 1219 break; 1220 default: 1221 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1222 } 1223 PetscFunctionReturn(0); 1224 } 1225 1226 static PetscErrorCode DMPlexComputeIsoparametricGeometry_Internal(DM dm, PetscFE fe, PetscInt point, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1227 { 1228 PetscQuadrature quad; 1229 PetscSection coordSection; 1230 Vec coordinates; 1231 PetscScalar *coords = NULL; 1232 const PetscReal *quadPoints; 1233 PetscReal *basisDer, *basis, detJt; 1234 PetscInt dim, cdim, pdim, qdim, Nq, numCoords, q; 1235 PetscErrorCode ierr; 1236 1237 PetscFunctionBegin; 1238 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1239 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1240 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1241 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1242 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 1243 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1244 ierr = PetscFEGetDimension(fe, &pdim);CHKERRQ(ierr); 1245 ierr = PetscQuadratureGetData(quad, &qdim, &Nq, &quadPoints, NULL);CHKERRQ(ierr); 1246 ierr = PetscFEGetDefaultTabulation(fe, &basis, &basisDer, NULL);CHKERRQ(ierr); 1247 if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim); 1248 if (numCoords != pdim*cdim) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "There are %d coordinates for point %d != %d*%d", numCoords, point, pdim, cdim); 1249 if (v0) { 1250 ierr = PetscMemzero(v0, Nq*cdim*sizeof(PetscReal));CHKERRQ(ierr); 1251 for (q = 0; q < Nq; ++q) { 1252 PetscInt i, k; 1253 1254 for (k = 0; k < pdim; ++k) 1255 for (i = 0; i < cdim; ++i) 1256 v0[q*cdim + i] += basis[q*pdim + k] * PetscRealPart(coords[k*cdim + i]); 1257 ierr = PetscLogFlops(2.0*pdim*cdim);CHKERRQ(ierr); 1258 } 1259 } 1260 if (J) { 1261 ierr = PetscMemzero(J, Nq*cdim*dim*sizeof(PetscReal));CHKERRQ(ierr); 1262 for (q = 0; q < Nq; ++q) { 1263 PetscInt i, j, k, c, r; 1264 1265 /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */ 1266 for (k = 0; k < pdim; ++k) 1267 for (j = 0; j < dim; ++j) 1268 for (i = 0; i < cdim; ++i) 1269 J[(q*cdim + i)*dim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]); 1270 ierr = PetscLogFlops(2.0*pdim*dim*cdim);CHKERRQ(ierr); 1271 if (cdim > dim) { 1272 for (c = dim; c < cdim; ++c) 1273 for (r = 0; r < cdim; ++r) 1274 J[r*cdim+c] = r == c ? 1.0 : 0.0; 1275 } 1276 if (!detJ && !invJ) continue; 1277 detJt = 0.; 1278 switch (cdim) { 1279 case 3: 1280 DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]); 1281 if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);} 1282 break; 1283 case 2: 1284 DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]); 1285 if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);} 1286 break; 1287 case 1: 1288 detJt = J[q*cdim*dim]; 1289 if (invJ) invJ[q*cdim*dim] = 1.0/detJt; 1290 } 1291 if (detJ) detJ[q] = detJt; 1292 } 1293 } 1294 else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ"); 1295 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1296 PetscFunctionReturn(0); 1297 } 1298 1299 /*@C 1300 DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell 1301 1302 Collective on DM 1303 1304 Input Arguments: 1305 + dm - the DM 1306 . cell - the cell 1307 - fe - the finite element containing the quadrature 1308 1309 Output Arguments: 1310 + v0 - if fe != NULL, the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element 1311 . J - the Jacobian of the transform from the reference element at each quadrature point 1312 . invJ - the inverse of the Jacobian at each quadrature point 1313 - detJ - the Jacobian determinant at each quadrature point 1314 1315 Level: advanced 1316 1317 Fortran Notes: 1318 Since it returns arrays, this routine is only available in Fortran 90, and you must 1319 include petsc.h90 in your code. 1320 1321 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1322 @*/ 1323 PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscFE fe, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1324 { 1325 PetscErrorCode ierr; 1326 1327 PetscFunctionBegin; 1328 if (!fe) {ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);} 1329 else {ierr = DMPlexComputeIsoparametricGeometry_Internal(dm, fe, cell, v0, J, invJ, detJ);CHKERRQ(ierr);} 1330 PetscFunctionReturn(0); 1331 } 1332 1333 static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1334 { 1335 PetscSection coordSection; 1336 Vec coordinates; 1337 PetscScalar *coords = NULL; 1338 PetscScalar tmp[2]; 1339 PetscInt coordSize; 1340 PetscErrorCode ierr; 1341 1342 PetscFunctionBegin; 1343 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1344 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1345 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1346 if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now"); 1347 ierr = DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);CHKERRQ(ierr); 1348 if (centroid) { 1349 centroid[0] = 0.5*PetscRealPart(coords[0] + tmp[0]); 1350 centroid[1] = 0.5*PetscRealPart(coords[1] + tmp[1]); 1351 } 1352 if (normal) { 1353 PetscReal norm; 1354 1355 normal[0] = -PetscRealPart(coords[1] - tmp[1]); 1356 normal[1] = PetscRealPart(coords[0] - tmp[0]); 1357 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]); 1358 normal[0] /= norm; 1359 normal[1] /= norm; 1360 } 1361 if (vol) { 1362 *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - tmp[0])) + PetscSqr(PetscRealPart(coords[1] - tmp[1]))); 1363 } 1364 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1365 PetscFunctionReturn(0); 1366 } 1367 1368 /* Centroid_i = (\sum_n A_n Cn_i ) / A */ 1369 static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1370 { 1371 PetscSection coordSection; 1372 Vec coordinates; 1373 PetscScalar *coords = NULL; 1374 PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9]; 1375 PetscInt tdim = 2, coordSize, numCorners, p, d, e; 1376 PetscErrorCode ierr; 1377 1378 PetscFunctionBegin; 1379 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1380 ierr = DMPlexGetConeSize(dm, cell, &numCorners);CHKERRQ(ierr); 1381 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1382 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1383 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 1384 if (dim > 2 && centroid) { 1385 v0[0] = PetscRealPart(coords[0]); 1386 v0[1] = PetscRealPart(coords[1]); 1387 v0[2] = PetscRealPart(coords[2]); 1388 } 1389 if (normal) { 1390 if (dim > 2) { 1391 const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]); 1392 const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]); 1393 const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]); 1394 PetscReal norm; 1395 1396 normal[0] = y0*z1 - z0*y1; 1397 normal[1] = z0*x1 - x0*z1; 1398 normal[2] = x0*y1 - y0*x1; 1399 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]); 1400 normal[0] /= norm; 1401 normal[1] /= norm; 1402 normal[2] /= norm; 1403 } else { 1404 for (d = 0; d < dim; ++d) normal[d] = 0.0; 1405 } 1406 } 1407 if (dim == 3) {ierr = DMPlexComputeProjection3Dto2D(coordSize, coords, R);CHKERRQ(ierr);} 1408 for (p = 0; p < numCorners; ++p) { 1409 /* Need to do this copy to get types right */ 1410 for (d = 0; d < tdim; ++d) { 1411 ctmp[d] = PetscRealPart(coords[p*tdim+d]); 1412 ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]); 1413 } 1414 Volume_Triangle_Origin_Internal(&vtmp, ctmp); 1415 vsum += vtmp; 1416 for (d = 0; d < tdim; ++d) { 1417 csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp; 1418 } 1419 } 1420 for (d = 0; d < tdim; ++d) { 1421 csum[d] /= (tdim+1)*vsum; 1422 } 1423 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1424 if (vol) *vol = PetscAbsReal(vsum); 1425 if (centroid) { 1426 if (dim > 2) { 1427 for (d = 0; d < dim; ++d) { 1428 centroid[d] = v0[d]; 1429 for (e = 0; e < dim; ++e) { 1430 centroid[d] += R[d*dim+e]*csum[e]; 1431 } 1432 } 1433 } else for (d = 0; d < dim; ++d) centroid[d] = csum[d]; 1434 } 1435 PetscFunctionReturn(0); 1436 } 1437 1438 /* Centroid_i = (\sum_n V_n Cn_i ) / V */ 1439 static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1440 { 1441 PetscSection coordSection; 1442 Vec coordinates; 1443 PetscScalar *coords = NULL; 1444 PetscReal vsum = 0.0, vtmp, coordsTmp[3*3]; 1445 const PetscInt *faces, *facesO; 1446 PetscInt numFaces, f, coordSize, numCorners, p, d; 1447 PetscErrorCode ierr; 1448 1449 PetscFunctionBegin; 1450 if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim); 1451 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1452 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1453 1454 if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0; 1455 ierr = DMPlexGetConeSize(dm, cell, &numFaces);CHKERRQ(ierr); 1456 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 1457 ierr = DMPlexGetConeOrientation(dm, cell, &facesO);CHKERRQ(ierr); 1458 for (f = 0; f < numFaces; ++f) { 1459 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 1460 numCorners = coordSize/dim; 1461 switch (numCorners) { 1462 case 3: 1463 for (d = 0; d < dim; ++d) { 1464 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 1465 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 1466 coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]); 1467 } 1468 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1469 if (facesO[f] < 0) vtmp = -vtmp; 1470 vsum += vtmp; 1471 if (centroid) { /* Centroid of OABC = (a+b+c)/4 */ 1472 for (d = 0; d < dim; ++d) { 1473 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1474 } 1475 } 1476 break; 1477 case 4: 1478 /* DO FOR PYRAMID */ 1479 /* First tet */ 1480 for (d = 0; d < dim; ++d) { 1481 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 1482 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 1483 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 1484 } 1485 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1486 if (facesO[f] < 0) vtmp = -vtmp; 1487 vsum += vtmp; 1488 if (centroid) { 1489 for (d = 0; d < dim; ++d) { 1490 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1491 } 1492 } 1493 /* Second tet */ 1494 for (d = 0; d < dim; ++d) { 1495 coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]); 1496 coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]); 1497 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 1498 } 1499 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1500 if (facesO[f] < 0) vtmp = -vtmp; 1501 vsum += vtmp; 1502 if (centroid) { 1503 for (d = 0; d < dim; ++d) { 1504 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1505 } 1506 } 1507 break; 1508 default: 1509 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners); 1510 } 1511 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 1512 } 1513 if (vol) *vol = PetscAbsReal(vsum); 1514 if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0; 1515 if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4); 1516 PetscFunctionReturn(0); 1517 } 1518 1519 /*@C 1520 DMPlexComputeCellGeometryFVM - Compute the volume for a given cell 1521 1522 Collective on DM 1523 1524 Input Arguments: 1525 + dm - the DM 1526 - cell - the cell 1527 1528 Output Arguments: 1529 + volume - the cell volume 1530 . centroid - the cell centroid 1531 - normal - the cell normal, if appropriate 1532 1533 Level: advanced 1534 1535 Fortran Notes: 1536 Since it returns arrays, this routine is only available in Fortran 90, and you must 1537 include petsc.h90 in your code. 1538 1539 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1540 @*/ 1541 PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1542 { 1543 PetscInt depth, dim; 1544 PetscErrorCode ierr; 1545 1546 PetscFunctionBegin; 1547 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1548 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1549 if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 1550 /* We need to keep a pointer to the depth label */ 1551 ierr = DMGetLabelValue(dm, "depth", cell, &depth);CHKERRQ(ierr); 1552 /* Cone size is now the number of faces */ 1553 switch (depth) { 1554 case 1: 1555 ierr = DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1556 break; 1557 case 2: 1558 ierr = DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1559 break; 1560 case 3: 1561 ierr = DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1562 break; 1563 default: 1564 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1565 } 1566 PetscFunctionReturn(0); 1567 } 1568 1569 /* This should also take a PetscFE argument I think */ 1570 PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom) 1571 { 1572 DM dmCell; 1573 Vec coordinates; 1574 PetscSection coordSection, sectionCell; 1575 PetscScalar *cgeom; 1576 PetscInt cStart, cEnd, cMax, c; 1577 PetscErrorCode ierr; 1578 1579 PetscFunctionBegin; 1580 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1581 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1582 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1583 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1584 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1585 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1586 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1587 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1588 cEnd = cMax < 0 ? cEnd : cMax; 1589 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1590 /* TODO This needs to be multiplied by Nq for non-affine */ 1591 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFECellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1592 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1593 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1594 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1595 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1596 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1597 for (c = cStart; c < cEnd; ++c) { 1598 PetscFECellGeom *cg; 1599 1600 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1601 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1602 ierr = DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v0, cg->J, cg->invJ, &cg->detJ);CHKERRQ(ierr); 1603 if (cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c); 1604 } 1605 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1606 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1607 PetscFunctionReturn(0); 1608 } 1609 1610 /*@ 1611 DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method 1612 1613 Input Parameter: 1614 . dm - The DM 1615 1616 Output Parameters: 1617 + cellgeom - A Vec of PetscFVCellGeom data 1618 . facegeom - A Vec of PetscFVFaceGeom data 1619 1620 Level: developer 1621 1622 .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM() 1623 @*/ 1624 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 1625 { 1626 DM dmFace, dmCell; 1627 DMLabel ghostLabel; 1628 PetscSection sectionFace, sectionCell; 1629 PetscSection coordSection; 1630 Vec coordinates; 1631 PetscScalar *fgeom, *cgeom; 1632 PetscReal minradius, gminradius; 1633 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 1634 PetscErrorCode ierr; 1635 1636 PetscFunctionBegin; 1637 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1638 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1639 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1640 /* Make cell centroids and volumes */ 1641 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1642 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1643 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1644 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1645 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1646 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1647 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1648 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1649 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1650 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1651 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1652 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1653 if (cEndInterior < 0) { 1654 cEndInterior = cEnd; 1655 } 1656 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1657 for (c = cStart; c < cEndInterior; ++c) { 1658 PetscFVCellGeom *cg; 1659 1660 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1661 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1662 ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr); 1663 } 1664 /* Compute face normals and minimum cell radius */ 1665 ierr = DMClone(dm, &dmFace);CHKERRQ(ierr); 1666 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);CHKERRQ(ierr); 1667 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1668 ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr); 1669 for (f = fStart; f < fEnd; ++f) {ierr = PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1670 ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr); 1671 ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr); 1672 ierr = PetscSectionDestroy(§ionFace);CHKERRQ(ierr); 1673 ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr); 1674 ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr); 1675 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1676 minradius = PETSC_MAX_REAL; 1677 for (f = fStart; f < fEnd; ++f) { 1678 PetscFVFaceGeom *fg; 1679 PetscReal area; 1680 PetscInt ghost = -1, d, numChildren; 1681 1682 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1683 ierr = DMPlexGetTreeChildren(dm,f,&numChildren,NULL);CHKERRQ(ierr); 1684 if (ghost >= 0 || numChildren) continue; 1685 ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr); 1686 ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);CHKERRQ(ierr); 1687 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 1688 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 1689 { 1690 PetscFVCellGeom *cL, *cR; 1691 PetscInt ncells; 1692 const PetscInt *cells; 1693 PetscReal *lcentroid, *rcentroid; 1694 PetscReal l[3], r[3], v[3]; 1695 1696 ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr); 1697 ierr = DMPlexGetSupportSize(dm, f, &ncells);CHKERRQ(ierr); 1698 ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr); 1699 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 1700 if (ncells > 1) { 1701 ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr); 1702 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 1703 } 1704 else { 1705 rcentroid = fg->centroid; 1706 } 1707 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);CHKERRQ(ierr); 1708 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);CHKERRQ(ierr); 1709 DMPlex_WaxpyD_Internal(dim, -1, l, r, v); 1710 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 1711 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 1712 } 1713 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 1714 if (dim == 2) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g) v (%g,%g)", f, (double) fg->normal[0], (double) fg->normal[1], (double) v[0], (double) v[1]); 1715 if (dim == 3) SETERRQ7(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", f, (double) fg->normal[0], (double) fg->normal[1], (double) fg->normal[2], (double) v[0], (double) v[1], (double) v[2]); 1716 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f); 1717 } 1718 if (cells[0] < cEndInterior) { 1719 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 1720 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1721 } 1722 if (ncells > 1 && cells[1] < cEndInterior) { 1723 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 1724 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1725 } 1726 } 1727 } 1728 ierr = MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1729 ierr = DMPlexSetMinRadius(dm, gminradius);CHKERRQ(ierr); 1730 /* Compute centroids of ghost cells */ 1731 for (c = cEndInterior; c < cEnd; ++c) { 1732 PetscFVFaceGeom *fg; 1733 const PetscInt *cone, *support; 1734 PetscInt coneSize, supportSize, s; 1735 1736 ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr); 1737 if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize); 1738 ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr); 1739 ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr); 1740 if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize); 1741 ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr); 1742 ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr); 1743 for (s = 0; s < 2; ++s) { 1744 /* Reflect ghost centroid across plane of face */ 1745 if (support[s] == c) { 1746 PetscFVCellGeom *ci; 1747 PetscFVCellGeom *cg; 1748 PetscReal c2f[3], a; 1749 1750 ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr); 1751 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 1752 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 1753 ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr); 1754 DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid); 1755 cg->volume = ci->volume; 1756 } 1757 } 1758 } 1759 ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr); 1760 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1761 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1762 ierr = DMDestroy(&dmFace);CHKERRQ(ierr); 1763 PetscFunctionReturn(0); 1764 } 1765 1766 /*@C 1767 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 1768 1769 Not collective 1770 1771 Input Argument: 1772 . dm - the DM 1773 1774 Output Argument: 1775 . minradius - the minium cell radius 1776 1777 Level: developer 1778 1779 .seealso: DMGetCoordinates() 1780 @*/ 1781 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 1782 { 1783 PetscFunctionBegin; 1784 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1785 PetscValidPointer(minradius,2); 1786 *minradius = ((DM_Plex*) dm->data)->minradius; 1787 PetscFunctionReturn(0); 1788 } 1789 1790 /*@C 1791 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 1792 1793 Logically collective 1794 1795 Input Arguments: 1796 + dm - the DM 1797 - minradius - the minium cell radius 1798 1799 Level: developer 1800 1801 .seealso: DMSetCoordinates() 1802 @*/ 1803 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 1804 { 1805 PetscFunctionBegin; 1806 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1807 ((DM_Plex*) dm->data)->minradius = minradius; 1808 PetscFunctionReturn(0); 1809 } 1810 1811 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 1812 { 1813 DMLabel ghostLabel; 1814 PetscScalar *dx, *grad, **gref; 1815 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 1816 PetscErrorCode ierr; 1817 1818 PetscFunctionBegin; 1819 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1820 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1821 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1822 ierr = DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);CHKERRQ(ierr); 1823 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 1824 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1825 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 1826 for (c = cStart; c < cEndInterior; c++) { 1827 const PetscInt *faces; 1828 PetscInt numFaces, usedFaces, f, d; 1829 PetscFVCellGeom *cg; 1830 PetscBool boundary; 1831 PetscInt ghost; 1832 1833 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1834 ierr = DMPlexGetConeSize(dm, c, &numFaces);CHKERRQ(ierr); 1835 ierr = DMPlexGetCone(dm, c, &faces);CHKERRQ(ierr); 1836 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 1837 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 1838 PetscFVCellGeom *cg1; 1839 PetscFVFaceGeom *fg; 1840 const PetscInt *fcells; 1841 PetscInt ncell, side; 1842 1843 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 1844 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 1845 if ((ghost >= 0) || boundary) continue; 1846 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr); 1847 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 1848 ncell = fcells[!side]; /* the neighbor */ 1849 ierr = DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr); 1850 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 1851 for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d]; 1852 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 1853 } 1854 if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 1855 ierr = PetscFVComputeGradient(fvm, usedFaces, dx, grad);CHKERRQ(ierr); 1856 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 1857 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 1858 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 1859 if ((ghost >= 0) || boundary) continue; 1860 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d]; 1861 ++usedFaces; 1862 } 1863 } 1864 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 1865 PetscFunctionReturn(0); 1866 } 1867 1868 static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 1869 { 1870 DMLabel ghostLabel; 1871 PetscScalar *dx, *grad, **gref; 1872 PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0; 1873 PetscSection neighSec; 1874 PetscInt (*neighbors)[2]; 1875 PetscInt *counter; 1876 PetscErrorCode ierr; 1877 1878 PetscFunctionBegin; 1879 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1880 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1881 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1882 if (cEndInterior < 0) { 1883 cEndInterior = cEnd; 1884 } 1885 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);CHKERRQ(ierr); 1886 ierr = PetscSectionSetChart(neighSec,cStart,cEndInterior);CHKERRQ(ierr); 1887 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1888 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1889 for (f = fStart; f < fEnd; f++) { 1890 const PetscInt *fcells; 1891 PetscBool boundary; 1892 PetscInt ghost = -1; 1893 PetscInt numChildren, numCells, c; 1894 1895 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1896 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 1897 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 1898 if ((ghost >= 0) || boundary || numChildren) continue; 1899 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 1900 if (numCells == 2) { 1901 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 1902 for (c = 0; c < 2; c++) { 1903 PetscInt cell = fcells[c]; 1904 1905 if (cell >= cStart && cell < cEndInterior) { 1906 ierr = PetscSectionAddDof(neighSec,cell,1);CHKERRQ(ierr); 1907 } 1908 } 1909 } 1910 } 1911 ierr = PetscSectionSetUp(neighSec);CHKERRQ(ierr); 1912 ierr = PetscSectionGetMaxDof(neighSec,&maxNumFaces);CHKERRQ(ierr); 1913 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 1914 nStart = 0; 1915 ierr = PetscSectionGetStorageSize(neighSec,&nEnd);CHKERRQ(ierr); 1916 ierr = PetscMalloc1((nEnd-nStart),&neighbors);CHKERRQ(ierr); 1917 ierr = PetscCalloc1((cEndInterior-cStart),&counter);CHKERRQ(ierr); 1918 for (f = fStart; f < fEnd; f++) { 1919 const PetscInt *fcells; 1920 PetscBool boundary; 1921 PetscInt ghost = -1; 1922 PetscInt numChildren, numCells, c; 1923 1924 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1925 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 1926 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 1927 if ((ghost >= 0) || boundary || numChildren) continue; 1928 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 1929 if (numCells == 2) { 1930 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 1931 for (c = 0; c < 2; c++) { 1932 PetscInt cell = fcells[c], off; 1933 1934 if (cell >= cStart && cell < cEndInterior) { 1935 ierr = PetscSectionGetOffset(neighSec,cell,&off);CHKERRQ(ierr); 1936 off += counter[cell - cStart]++; 1937 neighbors[off][0] = f; 1938 neighbors[off][1] = fcells[1 - c]; 1939 } 1940 } 1941 } 1942 } 1943 ierr = PetscFree(counter);CHKERRQ(ierr); 1944 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 1945 for (c = cStart; c < cEndInterior; c++) { 1946 PetscInt numFaces, f, d, off, ghost = -1; 1947 PetscFVCellGeom *cg; 1948 1949 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1950 ierr = PetscSectionGetDof(neighSec, c, &numFaces);CHKERRQ(ierr); 1951 ierr = PetscSectionGetOffset(neighSec, c, &off);CHKERRQ(ierr); 1952 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, c, &ghost);CHKERRQ(ierr);} 1953 if (ghost < 0 && numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 1954 for (f = 0; f < numFaces; ++f) { 1955 PetscFVCellGeom *cg1; 1956 PetscFVFaceGeom *fg; 1957 const PetscInt *fcells; 1958 PetscInt ncell, side, nface; 1959 1960 nface = neighbors[off + f][0]; 1961 ncell = neighbors[off + f][1]; 1962 ierr = DMPlexGetSupport(dm,nface,&fcells);CHKERRQ(ierr); 1963 side = (c != fcells[0]); 1964 ierr = DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);CHKERRQ(ierr); 1965 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 1966 for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d]; 1967 gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */ 1968 } 1969 ierr = PetscFVComputeGradient(fvm, numFaces, dx, grad);CHKERRQ(ierr); 1970 for (f = 0; f < numFaces; ++f) { 1971 for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d]; 1972 } 1973 } 1974 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 1975 ierr = PetscSectionDestroy(&neighSec);CHKERRQ(ierr); 1976 ierr = PetscFree(neighbors);CHKERRQ(ierr); 1977 PetscFunctionReturn(0); 1978 } 1979 1980 /*@ 1981 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 1982 1983 Collective on DM 1984 1985 Input Arguments: 1986 + dm - The DM 1987 . fvm - The PetscFV 1988 . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM() 1989 - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM() 1990 1991 Output Parameters: 1992 + faceGeometry - The geometric factors for gradient calculation are inserted 1993 - dmGrad - The DM describing the layout of gradient data 1994 1995 Level: developer 1996 1997 .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM() 1998 @*/ 1999 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 2000 { 2001 DM dmFace, dmCell; 2002 PetscScalar *fgeom, *cgeom; 2003 PetscSection sectionGrad, parentSection; 2004 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 2005 PetscErrorCode ierr; 2006 2007 PetscFunctionBegin; 2008 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2009 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 2010 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2011 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2012 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 2013 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 2014 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 2015 ierr = VecGetArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2016 ierr = VecGetArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2017 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 2018 if (!parentSection) { 2019 ierr = BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2020 } else { 2021 ierr = BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2022 } 2023 ierr = VecRestoreArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2024 ierr = VecRestoreArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2025 /* Create storage for gradients */ 2026 ierr = DMClone(dm, dmGrad);CHKERRQ(ierr); 2027 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 2028 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 2029 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 2030 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 2031 ierr = DMSetDefaultSection(*dmGrad, sectionGrad);CHKERRQ(ierr); 2032 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 2033 PetscFunctionReturn(0); 2034 } 2035 2036 PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM) 2037 { 2038 PetscObject cellgeomobj, facegeomobj; 2039 PetscErrorCode ierr; 2040 2041 PetscFunctionBegin; 2042 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2043 if (!cellgeomobj) { 2044 Vec cellgeomInt, facegeomInt; 2045 2046 ierr = DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);CHKERRQ(ierr); 2047 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);CHKERRQ(ierr); 2048 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);CHKERRQ(ierr); 2049 ierr = VecDestroy(&cellgeomInt);CHKERRQ(ierr); 2050 ierr = VecDestroy(&facegeomInt);CHKERRQ(ierr); 2051 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2052 } 2053 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);CHKERRQ(ierr); 2054 if (cellgeom) *cellgeom = (Vec) cellgeomobj; 2055 if (facegeom) *facegeom = (Vec) facegeomobj; 2056 if (gradDM) { 2057 PetscObject gradobj; 2058 PetscBool computeGradients; 2059 2060 ierr = PetscFVGetComputeGradients(fv,&computeGradients);CHKERRQ(ierr); 2061 if (!computeGradients) { 2062 *gradDM = NULL; 2063 PetscFunctionReturn(0); 2064 } 2065 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2066 if (!gradobj) { 2067 DM dmGradInt; 2068 2069 ierr = DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);CHKERRQ(ierr); 2070 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);CHKERRQ(ierr); 2071 ierr = DMDestroy(&dmGradInt);CHKERRQ(ierr); 2072 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2073 } 2074 *gradDM = (DM) gradobj; 2075 } 2076 PetscFunctionReturn(0); 2077 } 2078 2079 static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess) 2080 { 2081 PetscInt l, m; 2082 2083 PetscFunctionBeginHot; 2084 if (dimC == dimR && dimR <= 3) { 2085 /* invert Jacobian, multiply */ 2086 PetscScalar det, idet; 2087 2088 switch (dimR) { 2089 case 1: 2090 invJ[0] = 1./ J[0]; 2091 break; 2092 case 2: 2093 det = J[0] * J[3] - J[1] * J[2]; 2094 idet = 1./det; 2095 invJ[0] = J[3] * idet; 2096 invJ[1] = -J[1] * idet; 2097 invJ[2] = -J[2] * idet; 2098 invJ[3] = J[0] * idet; 2099 break; 2100 case 3: 2101 { 2102 invJ[0] = J[4] * J[8] - J[5] * J[7]; 2103 invJ[1] = J[2] * J[7] - J[1] * J[8]; 2104 invJ[2] = J[1] * J[5] - J[2] * J[4]; 2105 det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6]; 2106 idet = 1./det; 2107 invJ[0] *= idet; 2108 invJ[1] *= idet; 2109 invJ[2] *= idet; 2110 invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]); 2111 invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]); 2112 invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]); 2113 invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]); 2114 invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]); 2115 invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]); 2116 } 2117 break; 2118 } 2119 for (l = 0; l < dimR; l++) { 2120 for (m = 0; m < dimC; m++) { 2121 guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m]; 2122 } 2123 } 2124 } else { 2125 #if defined(PETSC_USE_COMPLEX) 2126 char transpose = 'C'; 2127 #else 2128 char transpose = 'T'; 2129 #endif 2130 PetscBLASInt m = dimR; 2131 PetscBLASInt n = dimC; 2132 PetscBLASInt one = 1; 2133 PetscBLASInt worksize = dimR * dimC, info; 2134 2135 for (l = 0; l < dimC; l++) {invJ[l] = resNeg[l];} 2136 2137 PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&m,&n,&one,J,&m,invJ,&n,work,&worksize, &info)); 2138 if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS"); 2139 2140 for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);} 2141 } 2142 PetscFunctionReturn(0); 2143 } 2144 2145 static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2146 { 2147 PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR); 2148 PetscScalar *coordsScalar = NULL; 2149 PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg; 2150 PetscScalar *J, *invJ, *work; 2151 PetscErrorCode ierr; 2152 2153 PetscFunctionBegin; 2154 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2155 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2156 if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);CHKERRQ(ierr); 2157 ierr = DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, PETSC_REAL, &cellData);CHKERRQ(ierr); 2158 ierr = DMGetWorkArray(dm, 3 * dimR * dimC, PETSC_SCALAR, &J);CHKERRQ(ierr); 2159 cellCoords = &cellData[0]; 2160 cellCoeffs = &cellData[coordSize]; 2161 extJ = &cellData[2 * coordSize]; 2162 resNeg = &cellData[2 * coordSize + dimR]; 2163 invJ = &J[dimR * dimC]; 2164 work = &J[2 * dimR * dimC]; 2165 if (dimR == 2) { 2166 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 2167 2168 for (i = 0; i < 4; i++) { 2169 PetscInt plexI = zToPlex[i]; 2170 2171 for (j = 0; j < dimC; j++) { 2172 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2173 } 2174 } 2175 } else if (dimR == 3) { 2176 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2177 2178 for (i = 0; i < 8; i++) { 2179 PetscInt plexI = zToPlex[i]; 2180 2181 for (j = 0; j < dimC; j++) { 2182 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2183 } 2184 } 2185 } else { 2186 for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);} 2187 } 2188 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 2189 for (i = 0; i < dimR; i++) { 2190 PetscReal *swap; 2191 2192 for (j = 0; j < (numV / 2); j++) { 2193 for (k = 0; k < dimC; k++) { 2194 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 2195 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 2196 } 2197 } 2198 2199 if (i < dimR - 1) { 2200 swap = cellCoeffs; 2201 cellCoeffs = cellCoords; 2202 cellCoords = swap; 2203 } 2204 } 2205 ierr = PetscMemzero(refCoords,numPoints * dimR * sizeof (PetscReal));CHKERRQ(ierr); 2206 for (j = 0; j < numPoints; j++) { 2207 for (i = 0; i < maxIts; i++) { 2208 PetscReal *guess = &refCoords[dimR * j]; 2209 2210 /* compute -residual and Jacobian */ 2211 for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];} 2212 for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;} 2213 for (k = 0; k < numV; k++) { 2214 PetscReal extCoord = 1.; 2215 for (l = 0; l < dimR; l++) { 2216 PetscReal coord = guess[l]; 2217 PetscInt dep = (k & (1 << l)) >> l; 2218 2219 extCoord *= dep * coord + !dep; 2220 extJ[l] = dep; 2221 2222 for (m = 0; m < dimR; m++) { 2223 PetscReal coord = guess[m]; 2224 PetscInt dep = ((k & (1 << m)) >> m) && (m != l); 2225 PetscReal mult = dep * coord + !dep; 2226 2227 extJ[l] *= mult; 2228 } 2229 } 2230 for (l = 0; l < dimC; l++) { 2231 PetscReal coeff = cellCoeffs[dimC * k + l]; 2232 2233 resNeg[l] -= coeff * extCoord; 2234 for (m = 0; m < dimR; m++) { 2235 J[dimR * l + m] += coeff * extJ[m]; 2236 } 2237 } 2238 } 2239 #if 0 && defined(PETSC_USE_DEBUG) 2240 { 2241 PetscReal maxAbs = 0.; 2242 2243 for (l = 0; l < dimC; l++) { 2244 maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l])); 2245 } 2246 ierr = PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,maxAbs);CHKERRQ(ierr); 2247 } 2248 #endif 2249 2250 ierr = DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);CHKERRQ(ierr); 2251 } 2252 } 2253 ierr = DMRestoreWorkArray(dm, 3 * dimR * dimC, PETSC_SCALAR, &J);CHKERRQ(ierr); 2254 ierr = DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, PETSC_REAL, &cellData);CHKERRQ(ierr); 2255 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2256 PetscFunctionReturn(0); 2257 } 2258 2259 static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2260 { 2261 PetscInt coordSize, i, j, k, l, numV = (1 << dimR); 2262 PetscScalar *coordsScalar = NULL; 2263 PetscReal *cellData, *cellCoords, *cellCoeffs; 2264 PetscErrorCode ierr; 2265 2266 PetscFunctionBegin; 2267 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2268 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2269 if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);CHKERRQ(ierr); 2270 ierr = DMGetWorkArray(dm, 2 * coordSize, PETSC_REAL, &cellData);CHKERRQ(ierr); 2271 cellCoords = &cellData[0]; 2272 cellCoeffs = &cellData[coordSize]; 2273 if (dimR == 2) { 2274 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 2275 2276 for (i = 0; i < 4; i++) { 2277 PetscInt plexI = zToPlex[i]; 2278 2279 for (j = 0; j < dimC; j++) { 2280 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2281 } 2282 } 2283 } else if (dimR == 3) { 2284 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2285 2286 for (i = 0; i < 8; i++) { 2287 PetscInt plexI = zToPlex[i]; 2288 2289 for (j = 0; j < dimC; j++) { 2290 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2291 } 2292 } 2293 } else { 2294 for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);} 2295 } 2296 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 2297 for (i = 0; i < dimR; i++) { 2298 PetscReal *swap; 2299 2300 for (j = 0; j < (numV / 2); j++) { 2301 for (k = 0; k < dimC; k++) { 2302 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 2303 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 2304 } 2305 } 2306 2307 if (i < dimR - 1) { 2308 swap = cellCoeffs; 2309 cellCoeffs = cellCoords; 2310 cellCoords = swap; 2311 } 2312 } 2313 ierr = PetscMemzero(realCoords,numPoints * dimC * sizeof (PetscReal));CHKERRQ(ierr); 2314 for (j = 0; j < numPoints; j++) { 2315 const PetscReal *guess = &refCoords[dimR * j]; 2316 PetscReal *mapped = &realCoords[dimC * j]; 2317 2318 for (k = 0; k < numV; k++) { 2319 PetscReal extCoord = 1.; 2320 for (l = 0; l < dimR; l++) { 2321 PetscReal coord = guess[l]; 2322 PetscInt dep = (k & (1 << l)) >> l; 2323 2324 extCoord *= dep * coord + !dep; 2325 } 2326 for (l = 0; l < dimC; l++) { 2327 PetscReal coeff = cellCoeffs[dimC * k + l]; 2328 2329 mapped[l] += coeff * extCoord; 2330 } 2331 } 2332 } 2333 ierr = DMRestoreWorkArray(dm, 2 * coordSize, PETSC_REAL, &cellData);CHKERRQ(ierr); 2334 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2335 PetscFunctionReturn(0); 2336 } 2337 2338 static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2339 { 2340 PetscInt numComp, numDof, i, j, k, l, m, maxIter = 7, coordSize; 2341 PetscScalar *nodes = NULL; 2342 PetscReal *invV, *modes; 2343 PetscReal *B, *D, *resNeg; 2344 PetscScalar *J, *invJ, *work; 2345 PetscErrorCode ierr; 2346 2347 PetscFunctionBegin; 2348 ierr = PetscFEGetDimension(fe, &numDof);CHKERRQ(ierr); 2349 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 2350 if (numComp != dimC) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,dimC); 2351 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2352 /* convert nodes to values in the stable evaluation basis */ 2353 ierr = DMGetWorkArray(dm,dimC * numDof,PETSC_REAL,&modes);CHKERRQ(ierr); 2354 invV = fe->invV; 2355 for (i = 0; i < numDof; i++) { 2356 for (j = 0; j < dimC; j++) { 2357 modes[i * dimC + j] = 0.; 2358 for (k = 0; k < numDof; k++) { 2359 modes[i * dimC + j] += invV[i * numDof + k] * PetscRealPart(nodes[k * dimC + j]); 2360 } 2361 } 2362 } 2363 ierr = DMGetWorkArray(dm,numDof + numDof * dimR + dimC,PETSC_REAL,&B);CHKERRQ(ierr); 2364 D = &B[numDof]; 2365 resNeg = &D[numDof * dimR]; 2366 ierr = DMGetWorkArray(dm,3 * dimC * dimR,PETSC_SCALAR,&J);CHKERRQ(ierr); 2367 invJ = &J[dimC * dimR]; 2368 work = &invJ[dimC * dimR]; 2369 for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;} 2370 for (j = 0; j < numPoints; j++) { 2371 for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */ 2372 PetscReal *guess = &refCoords[j * dimR]; 2373 ierr = PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);CHKERRQ(ierr); 2374 for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[j * dimC + k];} 2375 for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;} 2376 for (k = 0; k < numDof; k++) { 2377 for (l = 0; l < dimC; l++) { 2378 resNeg[l] -= modes[k * dimC + l] * B[k]; 2379 for (m = 0; m < dimR; m++) { 2380 J[l * dimR + m] += modes[k * dimC + l] * D[k * dimR + m]; 2381 } 2382 } 2383 } 2384 #if 0 && defined(PETSC_USE_DEBUG) 2385 { 2386 PetscReal maxAbs = 0.; 2387 2388 for (l = 0; l < dimC; l++) { 2389 maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l])); 2390 } 2391 ierr = PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,maxAbs);CHKERRQ(ierr); 2392 } 2393 #endif 2394 ierr = DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);CHKERRQ(ierr); 2395 } 2396 } 2397 ierr = DMRestoreWorkArray(dm,3 * dimC * dimR,PETSC_SCALAR,&J);CHKERRQ(ierr); 2398 ierr = DMRestoreWorkArray(dm,numDof + numDof * dimR + dimC,PETSC_REAL,&B);CHKERRQ(ierr); 2399 ierr = DMRestoreWorkArray(dm,dimC * numDof,PETSC_REAL,&modes);CHKERRQ(ierr); 2400 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2401 PetscFunctionReturn(0); 2402 } 2403 2404 static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2405 { 2406 PetscInt numComp, numDof, i, j, k, l, coordSize; 2407 PetscScalar *nodes = NULL; 2408 PetscReal *invV, *modes; 2409 PetscReal *B; 2410 PetscErrorCode ierr; 2411 2412 PetscFunctionBegin; 2413 ierr = PetscFEGetDimension(fe, &numDof);CHKERRQ(ierr); 2414 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 2415 if (numComp != dimC) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,dimC); 2416 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2417 /* convert nodes to values in the stable evaluation basis */ 2418 ierr = DMGetWorkArray(dm,dimC * numDof,PETSC_REAL,&modes);CHKERRQ(ierr); 2419 invV = fe->invV; 2420 for (i = 0; i < numDof; i++) { 2421 for (j = 0; j < dimC; j++) { 2422 modes[i * dimC + j] = 0.; 2423 for (k = 0; k < numDof; k++) { 2424 modes[i * dimC + j] += invV[i * numDof + k] * PetscRealPart(nodes[k * dimC + j]); 2425 } 2426 } 2427 } 2428 ierr = DMGetWorkArray(dm,numDof,PETSC_REAL,&B);CHKERRQ(ierr); 2429 for (i = 0; i < numPoints * dimC; i++) {realCoords[i] = 0.;} 2430 for (j = 0; j < numPoints; j++) { 2431 const PetscReal *guess = &refCoords[j * dimR]; 2432 PetscReal *mapped = &realCoords[j * dimC]; 2433 2434 ierr = PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, NULL, NULL);CHKERRQ(ierr); 2435 for (k = 0; k < numDof; k++) { 2436 for (l = 0; l < dimC; l++) { 2437 mapped[l] += modes[k * dimC + l] * B[k]; 2438 } 2439 } 2440 } 2441 ierr = DMRestoreWorkArray(dm,numDof,PETSC_REAL,&B);CHKERRQ(ierr); 2442 ierr = DMRestoreWorkArray(dm,dimC * numDof,PETSC_REAL,&modes);CHKERRQ(ierr); 2443 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2444 PetscFunctionReturn(0); 2445 } 2446 2447 /*@ 2448 DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element 2449 map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not 2450 extend uniquely outside the reference cell (e.g, most non-affine maps) 2451 2452 Not collective 2453 2454 Input Parameters: 2455 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 2456 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 2457 as a multilinear map for tensor-product elements 2458 . cell - the cell whose map is used. 2459 . numPoints - the number of points to locate 2460 - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 2461 2462 Output Parameters: 2463 . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 2464 2465 Level: intermediate 2466 @*/ 2467 PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[]) 2468 { 2469 PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i; 2470 DM coordDM = NULL; 2471 Vec coords; 2472 PetscFE fe = NULL; 2473 PetscErrorCode ierr; 2474 2475 PetscFunctionBegin; 2476 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2477 ierr = DMGetDimension(dm,&dimR);CHKERRQ(ierr); 2478 ierr = DMGetCoordinateDim(dm,&dimC);CHKERRQ(ierr); 2479 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 2480 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 2481 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2482 ierr = DMGetCoordinateDM(dm,&coordDM);CHKERRQ(ierr); 2483 if (coordDM) { 2484 PetscInt coordFields; 2485 2486 ierr = DMGetNumFields(coordDM,&coordFields);CHKERRQ(ierr); 2487 if (coordFields) { 2488 PetscClassId id; 2489 PetscObject disc; 2490 2491 ierr = DMGetField(coordDM,0,&disc);CHKERRQ(ierr); 2492 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 2493 if (id == PETSCFE_CLASSID) { 2494 fe = (PetscFE) disc; 2495 } 2496 } 2497 } 2498 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 2499 ierr = DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 2500 cEnd = cEndInterior > 0 ? cEndInterior : cEnd; 2501 if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);CHKERRQ(ierr); 2502 if (!fe) { /* implicit discretization: affine or multilinear */ 2503 PetscInt coneSize; 2504 PetscBool isSimplex, isTensor; 2505 2506 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 2507 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 2508 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 2509 if (isSimplex) { 2510 PetscReal detJ, *v0, *J, *invJ; 2511 2512 ierr = DMGetWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2513 J = &v0[dimC]; 2514 invJ = &J[dimC * dimC]; 2515 ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);CHKERRQ(ierr); 2516 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 2517 CoordinatesRealToRef(dimC, dimR, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]); 2518 } 2519 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2520 } else if (isTensor) { 2521 ierr = DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);CHKERRQ(ierr); 2522 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize); 2523 } else { 2524 ierr = DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);CHKERRQ(ierr); 2525 } 2526 PetscFunctionReturn(0); 2527 } 2528 2529 /*@ 2530 DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map. 2531 2532 Not collective 2533 2534 Input Parameters: 2535 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 2536 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 2537 as a multilinear map for tensor-product elements 2538 . cell - the cell whose map is used. 2539 . numPoints - the number of points to locate 2540 + refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 2541 2542 Output Parameters: 2543 . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 2544 2545 Level: intermediate 2546 @*/ 2547 PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[]) 2548 { 2549 PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i; 2550 DM coordDM = NULL; 2551 Vec coords; 2552 PetscFE fe = NULL; 2553 PetscErrorCode ierr; 2554 2555 PetscFunctionBegin; 2556 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2557 ierr = DMGetDimension(dm,&dimR);CHKERRQ(ierr); 2558 ierr = DMGetCoordinateDim(dm,&dimC);CHKERRQ(ierr); 2559 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 2560 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 2561 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2562 ierr = DMGetCoordinateDM(dm,&coordDM);CHKERRQ(ierr); 2563 if (coordDM) { 2564 PetscInt coordFields; 2565 2566 ierr = DMGetNumFields(coordDM,&coordFields);CHKERRQ(ierr); 2567 if (coordFields) { 2568 PetscClassId id; 2569 PetscObject disc; 2570 2571 ierr = DMGetField(coordDM,0,&disc);CHKERRQ(ierr); 2572 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 2573 if (id == PETSCFE_CLASSID) { 2574 fe = (PetscFE) disc; 2575 } 2576 } 2577 } 2578 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 2579 ierr = DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 2580 cEnd = cEndInterior > 0 ? cEndInterior : cEnd; 2581 if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);CHKERRQ(ierr); 2582 if (!fe) { /* implicit discretization: affine or multilinear */ 2583 PetscInt coneSize; 2584 PetscBool isSimplex, isTensor; 2585 2586 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 2587 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 2588 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 2589 if (isSimplex) { 2590 PetscReal detJ, *v0, *J; 2591 2592 ierr = DMGetWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2593 J = &v0[dimC]; 2594 ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ);CHKERRQ(ierr); 2595 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 2596 CoordinatesRefToReal(dimC, dimR, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]); 2597 } 2598 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2599 } else if (isTensor) { 2600 ierr = DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);CHKERRQ(ierr); 2601 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize); 2602 } else { 2603 ierr = DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);CHKERRQ(ierr); 2604 } 2605 PetscFunctionReturn(0); 2606 } 2607