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, PetscInt Nq, const PetscReal points[], PetscReal v[], 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 if (!Nq) { 1048 *detJ = 0.0; 1049 if (numCoords == 12) { 1050 const PetscInt dim = 3; 1051 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 1052 1053 if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);} 1054 ierr = DMPlexComputeProjection3Dto2D(numCoords, coords, R);CHKERRQ(ierr); 1055 if (J) { 1056 const PetscInt pdim = 2; 1057 1058 for (d = 0; d < pdim; d++) { 1059 J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 1060 J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 1061 } 1062 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1063 DMPlex_Det3D_Internal(detJ, J0); 1064 for (d = 0; d < dim; d++) { 1065 for (f = 0; f < dim; f++) { 1066 J[d*dim+f] = 0.0; 1067 for (g = 0; g < dim; g++) { 1068 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 1069 } 1070 } 1071 } 1072 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1073 } 1074 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1075 } else if (numCoords == 8) { 1076 const PetscInt dim = 2; 1077 1078 if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);} 1079 if (J) { 1080 for (d = 0; d < dim; d++) { 1081 J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1082 J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1083 } 1084 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1085 DMPlex_Det2D_Internal(detJ, J); 1086 } 1087 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 1088 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords); 1089 } else { 1090 const PetscInt Nv = 4; 1091 const PetscInt dimR = 2; 1092 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 1093 PetscReal zOrder[12]; 1094 PetscReal zCoeff[12]; 1095 PetscInt i, j, k, l, dim; 1096 1097 if (numCoords == 12) { 1098 dim = 3; 1099 } else if (numCoords == 8) { 1100 dim = 2; 1101 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords); 1102 for (i = 0; i < Nv; i++) { 1103 PetscInt zi = zToPlex[i]; 1104 1105 for (j = 0; j < dim; j++) { 1106 zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]); 1107 } 1108 } 1109 for (j = 0; j < dim; j++) { 1110 zCoeff[dim * 0 + j] = 0.25 * ( zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1111 zCoeff[dim * 1 + j] = 0.25 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1112 zCoeff[dim * 2 + j] = 0.25 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1113 zCoeff[dim * 3 + j] = 0.25 * ( zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1114 } 1115 for (i = 0; i < Nq; i++) { 1116 PetscReal xi = points[dimR * i], eta = points[dimR * i + 1]; 1117 1118 if (v) { 1119 PetscReal extPoint[4]; 1120 1121 extPoint[0] = 1.; 1122 extPoint[1] = xi; 1123 extPoint[2] = eta; 1124 extPoint[3] = xi * eta; 1125 for (j = 0; j < dim; j++) { 1126 PetscReal val = 0.; 1127 1128 for (k = 0; k < Nv; k++) { 1129 val += extPoint[k] * zCoeff[dim * k + j]; 1130 } 1131 v[i * dim + j] = val; 1132 } 1133 } 1134 if (J) { 1135 PetscReal extJ[8]; 1136 1137 extJ[0] = 0.; 1138 extJ[1] = 0.; 1139 extJ[2] = 1.; 1140 extJ[3] = 0.; 1141 extJ[4] = 0.; 1142 extJ[5] = 1.; 1143 extJ[6] = eta; 1144 extJ[7] = xi; 1145 for (j = 0; j < dim; j++) { 1146 for (k = 0; k < dimR; k++) { 1147 PetscReal val = 0.; 1148 1149 for (l = 0; l < Nv; l++) { 1150 val += zCoeff[dim * l + j] * extJ[dimR * l + k]; 1151 } 1152 J[i * dim * dim + dim * j + k] = val; 1153 } 1154 } 1155 if (dim == 3) { /* put the cross product in the third component of the Jacobian */ 1156 PetscReal x, y, z; 1157 PetscReal *iJ = &J[i * dim * dim]; 1158 PetscReal norm; 1159 1160 x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0]; 1161 y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1]; 1162 z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0]; 1163 norm = PetscSqrtReal(x * x + y * y + z * z); 1164 iJ[2] = x / norm; 1165 iJ[5] = y / norm; 1166 iJ[8] = z / norm; 1167 DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]); 1168 if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);} 1169 } else { 1170 DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]); 1171 if (invJ) {DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);} 1172 } 1173 } 1174 } 1175 } 1176 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1177 PetscFunctionReturn(0); 1178 } 1179 1180 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1181 { 1182 PetscSection coordSection; 1183 Vec coordinates; 1184 PetscScalar *coords = NULL; 1185 const PetscInt dim = 3; 1186 PetscInt d; 1187 PetscErrorCode ierr; 1188 1189 PetscFunctionBegin; 1190 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1191 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1192 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1193 *detJ = 0.0; 1194 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1195 if (J) { 1196 for (d = 0; d < dim; d++) { 1197 /* I orient with outward face normals */ 1198 J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d])); 1199 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1200 J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1201 } 1202 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1203 DMPlex_Det3D_Internal(detJ, J); 1204 } 1205 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1206 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1207 PetscFunctionReturn(0); 1208 } 1209 1210 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1211 { 1212 PetscSection coordSection; 1213 Vec coordinates; 1214 PetscScalar *coords = NULL; 1215 const PetscInt dim = 3; 1216 PetscInt d; 1217 PetscErrorCode ierr; 1218 1219 PetscFunctionBegin; 1220 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1221 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1222 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1223 if (!Nq) { 1224 *detJ = 0.0; 1225 if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);} 1226 if (J) { 1227 for (d = 0; d < dim; d++) { 1228 J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1229 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1230 J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d])); 1231 } 1232 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1233 DMPlex_Det3D_Internal(detJ, J); 1234 } 1235 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1236 } else { 1237 const PetscInt Nv = 8; 1238 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 1239 const PetscInt dim = 3; 1240 const PetscInt dimR = 3; 1241 PetscReal zOrder[24]; 1242 PetscReal zCoeff[24]; 1243 PetscInt i, j, k, l; 1244 1245 for (i = 0; i < Nv; i++) { 1246 PetscInt zi = zToPlex[i]; 1247 1248 for (j = 0; j < dim; j++) { 1249 zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]); 1250 } 1251 } 1252 for (j = 0; j < dim; j++) { 1253 zCoeff[dim * 0 + j] = 0.125 * ( zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j] + zOrder[dim * 4 + j] + zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]); 1254 zCoeff[dim * 1 + j] = 0.125 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j] - zOrder[dim * 4 + j] + zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]); 1255 zCoeff[dim * 2 + j] = 0.125 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j] - zOrder[dim * 4 + j] - zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]); 1256 zCoeff[dim * 3 + j] = 0.125 * ( zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j] + zOrder[dim * 4 + j] - zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]); 1257 zCoeff[dim * 4 + j] = 0.125 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] - zOrder[dim * 3 + j] + zOrder[dim * 4 + j] + zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]); 1258 zCoeff[dim * 5 + j] = 0.125 * (+ zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] - zOrder[dim * 3 + j] - zOrder[dim * 4 + j] + zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]); 1259 zCoeff[dim * 6 + j] = 0.125 * (+ zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] - zOrder[dim * 3 + j] - zOrder[dim * 4 + j] - zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]); 1260 zCoeff[dim * 7 + j] = 0.125 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] - zOrder[dim * 3 + j] + zOrder[dim * 4 + j] - zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]); 1261 } 1262 for (i = 0; i < Nq; i++) { 1263 PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2]; 1264 1265 if (v) { 1266 PetscReal extPoint[8]; 1267 1268 extPoint[0] = 1.; 1269 extPoint[1] = xi; 1270 extPoint[2] = eta; 1271 extPoint[3] = xi * eta; 1272 extPoint[4] = theta; 1273 extPoint[5] = theta * xi; 1274 extPoint[6] = theta * eta; 1275 extPoint[7] = theta * eta * xi; 1276 for (j = 0; j < dim; j++) { 1277 PetscReal val = 0.; 1278 1279 for (k = 0; k < Nv; k++) { 1280 val += extPoint[k] * zCoeff[dim * k + j]; 1281 } 1282 v[i * dim + j] = val; 1283 } 1284 } 1285 if (J) { 1286 PetscReal extJ[24]; 1287 1288 extJ[0] = 0. ; extJ[1] = 0. ; extJ[2] = 0. ; 1289 extJ[3] = 1. ; extJ[4] = 0. ; extJ[5] = 0. ; 1290 extJ[6] = 0. ; extJ[7] = 1. ; extJ[8] = 0. ; 1291 extJ[9] = eta ; extJ[10] = xi ; extJ[11] = 0. ; 1292 extJ[12] = 0. ; extJ[13] = 0. ; extJ[14] = 1. ; 1293 extJ[15] = theta ; extJ[16] = 0. ; extJ[17] = xi ; 1294 extJ[18] = 0. ; extJ[19] = theta ; extJ[20] = eta ; 1295 extJ[21] = theta * eta; extJ[22] = theta * xi; extJ[23] = eta * xi; 1296 1297 for (j = 0; j < dim; j++) { 1298 for (k = 0; k < dimR; k++) { 1299 PetscReal val = 0.; 1300 1301 for (l = 0; l < Nv; l++) { 1302 val += zCoeff[dim * l + j] * extJ[dimR * l + k]; 1303 } 1304 J[i * dim * dim + dim * j + k] = val; 1305 } 1306 } 1307 DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]); 1308 if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);} 1309 } 1310 } 1311 } 1312 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1313 PetscFunctionReturn(0); 1314 } 1315 1316 static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1317 { 1318 PetscInt depth, dim, coordDim, coneSize, i; 1319 PetscInt Nq = 0; 1320 const PetscReal *points = NULL; 1321 DMLabel depthLabel; 1322 PetscReal v0[3] = {-1.}, J0[9] = {-1.}, detJ0 = -1.; 1323 PetscBool isAffine = PETSC_TRUE; 1324 PetscErrorCode ierr; 1325 1326 PetscFunctionBegin; 1327 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1328 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 1329 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 1330 ierr = DMLabelGetValue(depthLabel, cell, &dim);CHKERRQ(ierr); 1331 if (depth == 1 && dim == 1) { 1332 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1333 } 1334 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1335 if (coordDim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %D > 3", coordDim); 1336 if (quad) {ierr = PetscQuadratureGetData(quad, NULL, &Nq, &points, NULL);CHKERRQ(ierr);} 1337 switch (dim) { 1338 case 0: 1339 ierr = DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);CHKERRQ(ierr); 1340 isAffine = PETSC_FALSE; 1341 break; 1342 case 1: 1343 if (Nq) { 1344 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);CHKERRQ(ierr); 1345 } else { 1346 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ);CHKERRQ(ierr); 1347 } 1348 break; 1349 case 2: 1350 switch (coneSize) { 1351 case 3: 1352 if (Nq) { 1353 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);CHKERRQ(ierr); 1354 } else { 1355 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ);CHKERRQ(ierr); 1356 } 1357 break; 1358 case 4: 1359 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);CHKERRQ(ierr); 1360 isAffine = PETSC_FALSE; 1361 break; 1362 default: 1363 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 1364 } 1365 break; 1366 case 3: 1367 switch (coneSize) { 1368 case 4: 1369 if (Nq) { 1370 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);CHKERRQ(ierr); 1371 } else { 1372 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ);CHKERRQ(ierr); 1373 } 1374 break; 1375 case 6: /* Faces */ 1376 case 8: /* Vertices */ 1377 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);CHKERRQ(ierr); 1378 isAffine = PETSC_FALSE; 1379 break; 1380 default: 1381 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 1382 } 1383 break; 1384 default: 1385 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1386 } 1387 if (isAffine && Nq) { 1388 if (v) { 1389 for (i = 0; i < Nq; i++) { 1390 CoordinatesRefToReal(coordDim,dim,v0, J0, &points[dim * i], &v[coordDim * i]); 1391 } 1392 } 1393 if (detJ) { 1394 for (i = 0; i < Nq; i++) { 1395 detJ[i] = detJ0; 1396 } 1397 } 1398 if (J) { 1399 PetscInt k; 1400 1401 for (i = 0, k = 0; i < Nq; i++) { 1402 PetscInt j; 1403 1404 for (j = 0; j < coordDim * coordDim; j++, k++) { 1405 J[k] = J0[j]; 1406 } 1407 } 1408 } 1409 if (invJ) { 1410 PetscInt k; 1411 switch (coordDim) { 1412 case 0: 1413 break; 1414 case 1: 1415 invJ[0] = 1./J0[0]; 1416 break; 1417 case 2: 1418 DMPlex_Invert2D_Internal(invJ, J0, detJ0); 1419 break; 1420 case 3: 1421 DMPlex_Invert3D_Internal(invJ, J0, detJ0); 1422 break; 1423 } 1424 for (i = 1, k = coordDim * coordDim; i < Nq; i++) { 1425 PetscInt j; 1426 1427 for (j = 0; j < coordDim * coordDim; j++, k++) { 1428 invJ[k] = invJ[j]; 1429 } 1430 } 1431 } 1432 } 1433 PetscFunctionReturn(0); 1434 } 1435 1436 /*@C 1437 DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell 1438 1439 Collective on DM 1440 1441 Input Arguments: 1442 + dm - the DM 1443 - cell - the cell 1444 1445 Output Arguments: 1446 + v0 - the translation part of this affine transform 1447 . J - the Jacobian of the transform from the reference element 1448 . invJ - the inverse of the Jacobian 1449 - detJ - the Jacobian determinant 1450 1451 Level: advanced 1452 1453 Fortran Notes: 1454 Since it returns arrays, this routine is only available in Fortran 90, and you must 1455 include petsc.h90 in your code. 1456 1457 .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinateVec() 1458 @*/ 1459 PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1460 { 1461 PetscErrorCode ierr; 1462 1463 PetscFunctionBegin; 1464 ierr = DMPlexComputeCellGeometryFEM_Implicit(dm,cell,NULL,v0,J,invJ,detJ);CHKERRQ(ierr); 1465 PetscFunctionReturn(0); 1466 } 1467 1468 #undef __FUNCT__ 1469 #define __FUNCT__ "DMPlexComputeCellGeometryFEM_FE" 1470 static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1471 { 1472 PetscQuadrature feQuad; 1473 PetscSection coordSection; 1474 Vec coordinates; 1475 PetscScalar *coords = NULL; 1476 const PetscReal *quadPoints; 1477 PetscReal *basisDer, *basis, detJt; 1478 PetscInt dim, cdim, pdim, qdim, Nq, numCoords, q; 1479 PetscErrorCode ierr; 1480 1481 PetscFunctionBegin; 1482 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1483 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1484 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1485 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1486 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 1487 if (!quad) { /* use the first point of the first functional of the dual space */ 1488 PetscDualSpace dsp; 1489 1490 ierr = PetscFEGetDualSpace(fe, &dsp);CHKERRQ(ierr); 1491 ierr = PetscDualSpaceGetFunctional(dsp, 0, &quad);CHKERRQ(ierr); 1492 ierr = PetscQuadratureGetData(quad, &qdim, &Nq, &quadPoints, NULL);CHKERRQ(ierr); 1493 Nq = 1; 1494 } else { 1495 ierr = PetscQuadratureGetData(quad, &qdim, &Nq, &quadPoints, NULL);CHKERRQ(ierr); 1496 } 1497 ierr = PetscFEGetDimension(fe, &pdim);CHKERRQ(ierr); 1498 ierr = PetscFEGetQuadrature(fe, &feQuad);CHKERRQ(ierr); 1499 if (feQuad == quad) { 1500 ierr = PetscFEGetDefaultTabulation(fe, &basis, &basisDer, NULL);CHKERRQ(ierr); 1501 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); 1502 } else { 1503 ierr = PetscFEGetTabulation(fe, Nq, quadPoints, &basis, &basisDer, NULL);CHKERRQ(ierr); 1504 } 1505 if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim); 1506 if (v) { 1507 ierr = PetscMemzero(v, Nq*cdim*sizeof(PetscReal));CHKERRQ(ierr); 1508 for (q = 0; q < Nq; ++q) { 1509 PetscInt i, k; 1510 1511 for (k = 0; k < pdim; ++k) 1512 for (i = 0; i < cdim; ++i) 1513 v[q*cdim + i] += basis[q*pdim + k] * PetscRealPart(coords[k*cdim + i]); 1514 ierr = PetscLogFlops(2.0*pdim*cdim);CHKERRQ(ierr); 1515 } 1516 } 1517 if (J) { 1518 ierr = PetscMemzero(J, Nq*cdim*cdim*sizeof(PetscReal));CHKERRQ(ierr); 1519 for (q = 0; q < Nq; ++q) { 1520 PetscInt i, j, k, c, r; 1521 1522 /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */ 1523 for (k = 0; k < pdim; ++k) 1524 for (j = 0; j < dim; ++j) 1525 for (i = 0; i < cdim; ++i) 1526 J[(q*cdim + i)*cdim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]); 1527 ierr = PetscLogFlops(2.0*pdim*dim*cdim);CHKERRQ(ierr); 1528 if (cdim > dim) { 1529 for (c = dim; c < cdim; ++c) 1530 for (r = 0; r < cdim; ++r) 1531 J[r*cdim+c] = r == c ? 1.0 : 0.0; 1532 } 1533 if (!detJ && !invJ) continue; 1534 detJt = 0.; 1535 switch (cdim) { 1536 case 3: 1537 DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]); 1538 if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);} 1539 break; 1540 case 2: 1541 DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]); 1542 if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);} 1543 break; 1544 case 1: 1545 detJt = J[q*cdim*dim]; 1546 if (invJ) invJ[q*cdim*dim] = 1.0/detJt; 1547 } 1548 if (detJ) detJ[q] = detJt; 1549 } 1550 } 1551 else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ"); 1552 if (feQuad != quad) { 1553 ierr = PetscFERestoreTabulation(fe, Nq, quadPoints, &basis, &basisDer, NULL);CHKERRQ(ierr); 1554 } 1555 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1556 PetscFunctionReturn(0); 1557 } 1558 1559 /*@C 1560 DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell 1561 1562 Collective on DM 1563 1564 Input Arguments: 1565 + dm - the DM 1566 . cell - the cell 1567 - quad - the quadrature containing the points in the reference element where the geometry will be evaluated. If quad == NULL, geometry will be 1568 evaluated at the first vertex of the reference element 1569 1570 Output Arguments: 1571 + v - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element 1572 . J - the Jacobian of the transform from the reference element at each quadrature point 1573 . invJ - the inverse of the Jacobian at each quadrature point 1574 - detJ - the Jacobian determinant at each quadrature point 1575 1576 Level: advanced 1577 1578 Fortran Notes: 1579 Since it returns arrays, this routine is only available in Fortran 90, and you must 1580 include petsc.h90 in your code. 1581 1582 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1583 @*/ 1584 PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1585 { 1586 PetscFE fe = NULL; 1587 PetscErrorCode ierr; 1588 1589 PetscFunctionBegin; 1590 if (dm->coordinateDM) { 1591 PetscClassId id; 1592 PetscInt numFields; 1593 PetscDS prob = dm->coordinateDM->prob; 1594 PetscObject disc; 1595 1596 ierr = PetscDSGetNumFields(prob, &numFields);CHKERRQ(ierr); 1597 if (numFields) { 1598 ierr = PetscDSGetDiscretization(prob,0,&disc);CHKERRQ(ierr); 1599 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 1600 if (id == PETSCFE_CLASSID) { 1601 fe = (PetscFE) disc; 1602 } 1603 } 1604 } 1605 if (!fe) {ierr = DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ);CHKERRQ(ierr);} 1606 else {ierr = DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ);CHKERRQ(ierr);} 1607 PetscFunctionReturn(0); 1608 } 1609 1610 static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1611 { 1612 PetscSection coordSection; 1613 Vec coordinates; 1614 PetscScalar *coords = NULL; 1615 PetscScalar tmp[2]; 1616 PetscInt coordSize; 1617 PetscErrorCode ierr; 1618 1619 PetscFunctionBegin; 1620 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1621 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1622 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1623 if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now"); 1624 ierr = DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);CHKERRQ(ierr); 1625 if (centroid) { 1626 centroid[0] = 0.5*PetscRealPart(coords[0] + tmp[0]); 1627 centroid[1] = 0.5*PetscRealPart(coords[1] + tmp[1]); 1628 } 1629 if (normal) { 1630 PetscReal norm; 1631 1632 normal[0] = -PetscRealPart(coords[1] - tmp[1]); 1633 normal[1] = PetscRealPart(coords[0] - tmp[0]); 1634 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]); 1635 normal[0] /= norm; 1636 normal[1] /= norm; 1637 } 1638 if (vol) { 1639 *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - tmp[0])) + PetscSqr(PetscRealPart(coords[1] - tmp[1]))); 1640 } 1641 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1642 PetscFunctionReturn(0); 1643 } 1644 1645 /* Centroid_i = (\sum_n A_n Cn_i ) / A */ 1646 static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1647 { 1648 PetscSection coordSection; 1649 Vec coordinates; 1650 PetscScalar *coords = NULL; 1651 PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9]; 1652 PetscInt tdim = 2, coordSize, numCorners, p, d, e; 1653 PetscErrorCode ierr; 1654 1655 PetscFunctionBegin; 1656 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1657 ierr = DMPlexGetConeSize(dm, cell, &numCorners);CHKERRQ(ierr); 1658 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1659 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1660 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 1661 if (dim > 2 && centroid) { 1662 v0[0] = PetscRealPart(coords[0]); 1663 v0[1] = PetscRealPart(coords[1]); 1664 v0[2] = PetscRealPart(coords[2]); 1665 } 1666 if (normal) { 1667 if (dim > 2) { 1668 const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]); 1669 const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]); 1670 const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]); 1671 PetscReal norm; 1672 1673 normal[0] = y0*z1 - z0*y1; 1674 normal[1] = z0*x1 - x0*z1; 1675 normal[2] = x0*y1 - y0*x1; 1676 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]); 1677 normal[0] /= norm; 1678 normal[1] /= norm; 1679 normal[2] /= norm; 1680 } else { 1681 for (d = 0; d < dim; ++d) normal[d] = 0.0; 1682 } 1683 } 1684 if (dim == 3) {ierr = DMPlexComputeProjection3Dto2D(coordSize, coords, R);CHKERRQ(ierr);} 1685 for (p = 0; p < numCorners; ++p) { 1686 /* Need to do this copy to get types right */ 1687 for (d = 0; d < tdim; ++d) { 1688 ctmp[d] = PetscRealPart(coords[p*tdim+d]); 1689 ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]); 1690 } 1691 Volume_Triangle_Origin_Internal(&vtmp, ctmp); 1692 vsum += vtmp; 1693 for (d = 0; d < tdim; ++d) { 1694 csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp; 1695 } 1696 } 1697 for (d = 0; d < tdim; ++d) { 1698 csum[d] /= (tdim+1)*vsum; 1699 } 1700 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1701 if (vol) *vol = PetscAbsReal(vsum); 1702 if (centroid) { 1703 if (dim > 2) { 1704 for (d = 0; d < dim; ++d) { 1705 centroid[d] = v0[d]; 1706 for (e = 0; e < dim; ++e) { 1707 centroid[d] += R[d*dim+e]*csum[e]; 1708 } 1709 } 1710 } else for (d = 0; d < dim; ++d) centroid[d] = csum[d]; 1711 } 1712 PetscFunctionReturn(0); 1713 } 1714 1715 /* Centroid_i = (\sum_n V_n Cn_i ) / V */ 1716 static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1717 { 1718 PetscSection coordSection; 1719 Vec coordinates; 1720 PetscScalar *coords = NULL; 1721 PetscReal vsum = 0.0, vtmp, coordsTmp[3*3]; 1722 const PetscInt *faces, *facesO; 1723 PetscInt numFaces, f, coordSize, numCorners, p, d; 1724 PetscErrorCode ierr; 1725 1726 PetscFunctionBegin; 1727 if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim); 1728 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1729 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1730 1731 if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0; 1732 ierr = DMPlexGetConeSize(dm, cell, &numFaces);CHKERRQ(ierr); 1733 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 1734 ierr = DMPlexGetConeOrientation(dm, cell, &facesO);CHKERRQ(ierr); 1735 for (f = 0; f < numFaces; ++f) { 1736 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 1737 numCorners = coordSize/dim; 1738 switch (numCorners) { 1739 case 3: 1740 for (d = 0; d < dim; ++d) { 1741 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 1742 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 1743 coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]); 1744 } 1745 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1746 if (facesO[f] < 0) vtmp = -vtmp; 1747 vsum += vtmp; 1748 if (centroid) { /* Centroid of OABC = (a+b+c)/4 */ 1749 for (d = 0; d < dim; ++d) { 1750 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1751 } 1752 } 1753 break; 1754 case 4: 1755 /* DO FOR PYRAMID */ 1756 /* First tet */ 1757 for (d = 0; d < dim; ++d) { 1758 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 1759 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 1760 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 1761 } 1762 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1763 if (facesO[f] < 0) vtmp = -vtmp; 1764 vsum += vtmp; 1765 if (centroid) { 1766 for (d = 0; d < dim; ++d) { 1767 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1768 } 1769 } 1770 /* Second tet */ 1771 for (d = 0; d < dim; ++d) { 1772 coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]); 1773 coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]); 1774 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 1775 } 1776 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1777 if (facesO[f] < 0) vtmp = -vtmp; 1778 vsum += vtmp; 1779 if (centroid) { 1780 for (d = 0; d < dim; ++d) { 1781 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1782 } 1783 } 1784 break; 1785 default: 1786 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners); 1787 } 1788 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 1789 } 1790 if (vol) *vol = PetscAbsReal(vsum); 1791 if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0; 1792 if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4); 1793 PetscFunctionReturn(0); 1794 } 1795 1796 /*@C 1797 DMPlexComputeCellGeometryFVM - Compute the volume for a given cell 1798 1799 Collective on DM 1800 1801 Input Arguments: 1802 + dm - the DM 1803 - cell - the cell 1804 1805 Output Arguments: 1806 + volume - the cell volume 1807 . centroid - the cell centroid 1808 - normal - the cell normal, if appropriate 1809 1810 Level: advanced 1811 1812 Fortran Notes: 1813 Since it returns arrays, this routine is only available in Fortran 90, and you must 1814 include petsc.h90 in your code. 1815 1816 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1817 @*/ 1818 PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1819 { 1820 PetscInt depth, dim; 1821 PetscErrorCode ierr; 1822 1823 PetscFunctionBegin; 1824 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1825 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1826 if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 1827 /* We need to keep a pointer to the depth label */ 1828 ierr = DMGetLabelValue(dm, "depth", cell, &depth);CHKERRQ(ierr); 1829 /* Cone size is now the number of faces */ 1830 switch (depth) { 1831 case 1: 1832 ierr = DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1833 break; 1834 case 2: 1835 ierr = DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1836 break; 1837 case 3: 1838 ierr = DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1839 break; 1840 default: 1841 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1842 } 1843 PetscFunctionReturn(0); 1844 } 1845 1846 /*@ 1847 DMPlexComputeGeometryFEM - Precompute cell geometry for the entire mesh 1848 1849 Collective on dm 1850 1851 Input Parameter: 1852 . dm - The DMPlex 1853 1854 Output Parameter: 1855 . cellgeom - A vector with the cell geometry data for each cell 1856 1857 Level: beginner 1858 1859 .keywords: DMPlexComputeCellGeometryFEM() 1860 @*/ 1861 PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom) 1862 { 1863 DM dmCell; 1864 Vec coordinates; 1865 PetscSection coordSection, sectionCell; 1866 PetscScalar *cgeom; 1867 PetscInt cStart, cEnd, cMax, c; 1868 PetscErrorCode ierr; 1869 1870 PetscFunctionBegin; 1871 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1872 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1873 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1874 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1875 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1876 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1877 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1878 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1879 cEnd = cMax < 0 ? cEnd : cMax; 1880 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1881 /* TODO This needs to be multiplied by Nq for non-affine */ 1882 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFECellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1883 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1884 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1885 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1886 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1887 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1888 for (c = cStart; c < cEnd; ++c) { 1889 PetscFECellGeom *cg; 1890 1891 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1892 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1893 ierr = DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v0, cg->J, cg->invJ, &cg->detJ);CHKERRQ(ierr); 1894 if (cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c); 1895 } 1896 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1897 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1898 PetscFunctionReturn(0); 1899 } 1900 1901 /*@ 1902 DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method 1903 1904 Input Parameter: 1905 . dm - The DM 1906 1907 Output Parameters: 1908 + cellgeom - A Vec of PetscFVCellGeom data 1909 . facegeom - A Vec of PetscFVFaceGeom data 1910 1911 Level: developer 1912 1913 .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM() 1914 @*/ 1915 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 1916 { 1917 DM dmFace, dmCell; 1918 DMLabel ghostLabel; 1919 PetscSection sectionFace, sectionCell; 1920 PetscSection coordSection; 1921 Vec coordinates; 1922 PetscScalar *fgeom, *cgeom; 1923 PetscReal minradius, gminradius; 1924 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 1925 PetscErrorCode ierr; 1926 1927 PetscFunctionBegin; 1928 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1929 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1930 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1931 /* Make cell centroids and volumes */ 1932 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1933 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1934 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1935 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1936 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1937 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1938 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1939 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1940 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1941 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1942 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1943 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1944 if (cEndInterior < 0) { 1945 cEndInterior = cEnd; 1946 } 1947 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1948 for (c = cStart; c < cEndInterior; ++c) { 1949 PetscFVCellGeom *cg; 1950 1951 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1952 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1953 ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr); 1954 } 1955 /* Compute face normals and minimum cell radius */ 1956 ierr = DMClone(dm, &dmFace);CHKERRQ(ierr); 1957 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);CHKERRQ(ierr); 1958 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1959 ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr); 1960 for (f = fStart; f < fEnd; ++f) {ierr = PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1961 ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr); 1962 ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr); 1963 ierr = PetscSectionDestroy(§ionFace);CHKERRQ(ierr); 1964 ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr); 1965 ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr); 1966 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1967 minradius = PETSC_MAX_REAL; 1968 for (f = fStart; f < fEnd; ++f) { 1969 PetscFVFaceGeom *fg; 1970 PetscReal area; 1971 PetscInt ghost = -1, d, numChildren; 1972 1973 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1974 ierr = DMPlexGetTreeChildren(dm,f,&numChildren,NULL);CHKERRQ(ierr); 1975 if (ghost >= 0 || numChildren) continue; 1976 ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr); 1977 ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);CHKERRQ(ierr); 1978 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 1979 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 1980 { 1981 PetscFVCellGeom *cL, *cR; 1982 PetscInt ncells; 1983 const PetscInt *cells; 1984 PetscReal *lcentroid, *rcentroid; 1985 PetscReal l[3], r[3], v[3]; 1986 1987 ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr); 1988 ierr = DMPlexGetSupportSize(dm, f, &ncells);CHKERRQ(ierr); 1989 ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr); 1990 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 1991 if (ncells > 1) { 1992 ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr); 1993 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 1994 } 1995 else { 1996 rcentroid = fg->centroid; 1997 } 1998 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);CHKERRQ(ierr); 1999 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);CHKERRQ(ierr); 2000 DMPlex_WaxpyD_Internal(dim, -1, l, r, v); 2001 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 2002 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 2003 } 2004 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 2005 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]); 2006 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]); 2007 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f); 2008 } 2009 if (cells[0] < cEndInterior) { 2010 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 2011 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2012 } 2013 if (ncells > 1 && cells[1] < cEndInterior) { 2014 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 2015 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2016 } 2017 } 2018 } 2019 ierr = MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 2020 ierr = DMPlexSetMinRadius(dm, gminradius);CHKERRQ(ierr); 2021 /* Compute centroids of ghost cells */ 2022 for (c = cEndInterior; c < cEnd; ++c) { 2023 PetscFVFaceGeom *fg; 2024 const PetscInt *cone, *support; 2025 PetscInt coneSize, supportSize, s; 2026 2027 ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr); 2028 if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize); 2029 ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr); 2030 ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr); 2031 if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize); 2032 ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr); 2033 ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr); 2034 for (s = 0; s < 2; ++s) { 2035 /* Reflect ghost centroid across plane of face */ 2036 if (support[s] == c) { 2037 PetscFVCellGeom *ci; 2038 PetscFVCellGeom *cg; 2039 PetscReal c2f[3], a; 2040 2041 ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr); 2042 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 2043 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 2044 ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr); 2045 DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid); 2046 cg->volume = ci->volume; 2047 } 2048 } 2049 } 2050 ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr); 2051 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 2052 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 2053 ierr = DMDestroy(&dmFace);CHKERRQ(ierr); 2054 PetscFunctionReturn(0); 2055 } 2056 2057 /*@C 2058 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 2059 2060 Not collective 2061 2062 Input Argument: 2063 . dm - the DM 2064 2065 Output Argument: 2066 . minradius - the minium cell radius 2067 2068 Level: developer 2069 2070 .seealso: DMGetCoordinates() 2071 @*/ 2072 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 2073 { 2074 PetscFunctionBegin; 2075 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2076 PetscValidPointer(minradius,2); 2077 *minradius = ((DM_Plex*) dm->data)->minradius; 2078 PetscFunctionReturn(0); 2079 } 2080 2081 /*@C 2082 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 2083 2084 Logically collective 2085 2086 Input Arguments: 2087 + dm - the DM 2088 - minradius - the minium cell radius 2089 2090 Level: developer 2091 2092 .seealso: DMSetCoordinates() 2093 @*/ 2094 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 2095 { 2096 PetscFunctionBegin; 2097 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2098 ((DM_Plex*) dm->data)->minradius = minradius; 2099 PetscFunctionReturn(0); 2100 } 2101 2102 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 2103 { 2104 DMLabel ghostLabel; 2105 PetscScalar *dx, *grad, **gref; 2106 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 2107 PetscErrorCode ierr; 2108 2109 PetscFunctionBegin; 2110 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2111 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2112 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2113 ierr = DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);CHKERRQ(ierr); 2114 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 2115 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 2116 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 2117 for (c = cStart; c < cEndInterior; c++) { 2118 const PetscInt *faces; 2119 PetscInt numFaces, usedFaces, f, d; 2120 PetscFVCellGeom *cg; 2121 PetscBool boundary; 2122 PetscInt ghost; 2123 2124 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2125 ierr = DMPlexGetConeSize(dm, c, &numFaces);CHKERRQ(ierr); 2126 ierr = DMPlexGetCone(dm, c, &faces);CHKERRQ(ierr); 2127 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 2128 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 2129 PetscFVCellGeom *cg1; 2130 PetscFVFaceGeom *fg; 2131 const PetscInt *fcells; 2132 PetscInt ncell, side; 2133 2134 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 2135 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 2136 if ((ghost >= 0) || boundary) continue; 2137 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr); 2138 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 2139 ncell = fcells[!side]; /* the neighbor */ 2140 ierr = DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr); 2141 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 2142 for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d]; 2143 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2144 } 2145 if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 2146 ierr = PetscFVComputeGradient(fvm, usedFaces, dx, grad);CHKERRQ(ierr); 2147 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 2148 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 2149 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 2150 if ((ghost >= 0) || boundary) continue; 2151 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d]; 2152 ++usedFaces; 2153 } 2154 } 2155 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 2156 PetscFunctionReturn(0); 2157 } 2158 2159 static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 2160 { 2161 DMLabel ghostLabel; 2162 PetscScalar *dx, *grad, **gref; 2163 PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0; 2164 PetscSection neighSec; 2165 PetscInt (*neighbors)[2]; 2166 PetscInt *counter; 2167 PetscErrorCode ierr; 2168 2169 PetscFunctionBegin; 2170 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2171 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2172 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2173 if (cEndInterior < 0) { 2174 cEndInterior = cEnd; 2175 } 2176 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);CHKERRQ(ierr); 2177 ierr = PetscSectionSetChart(neighSec,cStart,cEndInterior);CHKERRQ(ierr); 2178 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2179 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 2180 for (f = fStart; f < fEnd; f++) { 2181 const PetscInt *fcells; 2182 PetscBool boundary; 2183 PetscInt ghost = -1; 2184 PetscInt numChildren, numCells, c; 2185 2186 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 2187 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 2188 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 2189 if ((ghost >= 0) || boundary || numChildren) continue; 2190 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 2191 if (numCells == 2) { 2192 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 2193 for (c = 0; c < 2; c++) { 2194 PetscInt cell = fcells[c]; 2195 2196 if (cell >= cStart && cell < cEndInterior) { 2197 ierr = PetscSectionAddDof(neighSec,cell,1);CHKERRQ(ierr); 2198 } 2199 } 2200 } 2201 } 2202 ierr = PetscSectionSetUp(neighSec);CHKERRQ(ierr); 2203 ierr = PetscSectionGetMaxDof(neighSec,&maxNumFaces);CHKERRQ(ierr); 2204 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 2205 nStart = 0; 2206 ierr = PetscSectionGetStorageSize(neighSec,&nEnd);CHKERRQ(ierr); 2207 ierr = PetscMalloc1((nEnd-nStart),&neighbors);CHKERRQ(ierr); 2208 ierr = PetscCalloc1((cEndInterior-cStart),&counter);CHKERRQ(ierr); 2209 for (f = fStart; f < fEnd; f++) { 2210 const PetscInt *fcells; 2211 PetscBool boundary; 2212 PetscInt ghost = -1; 2213 PetscInt numChildren, numCells, c; 2214 2215 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 2216 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 2217 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 2218 if ((ghost >= 0) || boundary || numChildren) continue; 2219 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 2220 if (numCells == 2) { 2221 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 2222 for (c = 0; c < 2; c++) { 2223 PetscInt cell = fcells[c], off; 2224 2225 if (cell >= cStart && cell < cEndInterior) { 2226 ierr = PetscSectionGetOffset(neighSec,cell,&off);CHKERRQ(ierr); 2227 off += counter[cell - cStart]++; 2228 neighbors[off][0] = f; 2229 neighbors[off][1] = fcells[1 - c]; 2230 } 2231 } 2232 } 2233 } 2234 ierr = PetscFree(counter);CHKERRQ(ierr); 2235 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 2236 for (c = cStart; c < cEndInterior; c++) { 2237 PetscInt numFaces, f, d, off, ghost = -1; 2238 PetscFVCellGeom *cg; 2239 2240 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2241 ierr = PetscSectionGetDof(neighSec, c, &numFaces);CHKERRQ(ierr); 2242 ierr = PetscSectionGetOffset(neighSec, c, &off);CHKERRQ(ierr); 2243 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, c, &ghost);CHKERRQ(ierr);} 2244 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); 2245 for (f = 0; f < numFaces; ++f) { 2246 PetscFVCellGeom *cg1; 2247 PetscFVFaceGeom *fg; 2248 const PetscInt *fcells; 2249 PetscInt ncell, side, nface; 2250 2251 nface = neighbors[off + f][0]; 2252 ncell = neighbors[off + f][1]; 2253 ierr = DMPlexGetSupport(dm,nface,&fcells);CHKERRQ(ierr); 2254 side = (c != fcells[0]); 2255 ierr = DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);CHKERRQ(ierr); 2256 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 2257 for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d]; 2258 gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2259 } 2260 ierr = PetscFVComputeGradient(fvm, numFaces, dx, grad);CHKERRQ(ierr); 2261 for (f = 0; f < numFaces; ++f) { 2262 for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d]; 2263 } 2264 } 2265 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 2266 ierr = PetscSectionDestroy(&neighSec);CHKERRQ(ierr); 2267 ierr = PetscFree(neighbors);CHKERRQ(ierr); 2268 PetscFunctionReturn(0); 2269 } 2270 2271 /*@ 2272 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 2273 2274 Collective on DM 2275 2276 Input Arguments: 2277 + dm - The DM 2278 . fvm - The PetscFV 2279 . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM() 2280 - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM() 2281 2282 Output Parameters: 2283 + faceGeometry - The geometric factors for gradient calculation are inserted 2284 - dmGrad - The DM describing the layout of gradient data 2285 2286 Level: developer 2287 2288 .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM() 2289 @*/ 2290 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 2291 { 2292 DM dmFace, dmCell; 2293 PetscScalar *fgeom, *cgeom; 2294 PetscSection sectionGrad, parentSection; 2295 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 2296 PetscErrorCode ierr; 2297 2298 PetscFunctionBegin; 2299 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2300 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 2301 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2302 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2303 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 2304 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 2305 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 2306 ierr = VecGetArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2307 ierr = VecGetArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2308 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 2309 if (!parentSection) { 2310 ierr = BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2311 } else { 2312 ierr = BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2313 } 2314 ierr = VecRestoreArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2315 ierr = VecRestoreArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2316 /* Create storage for gradients */ 2317 ierr = DMClone(dm, dmGrad);CHKERRQ(ierr); 2318 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 2319 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 2320 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 2321 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 2322 ierr = DMSetDefaultSection(*dmGrad, sectionGrad);CHKERRQ(ierr); 2323 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 2324 PetscFunctionReturn(0); 2325 } 2326 2327 /*@ 2328 DMPlexGetDataFVM - Retrieve precomputed cell geometry 2329 2330 Collective on DM 2331 2332 Input Arguments: 2333 + dm - The DM 2334 - fvm - The PetscFV 2335 2336 Output Parameters: 2337 + cellGeometry - The cell geometry 2338 . faceGeometry - The face geometry 2339 - dmGrad - The gradient matrices 2340 2341 Level: developer 2342 2343 .seealso: DMPlexComputeGeometryFVM() 2344 @*/ 2345 PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM) 2346 { 2347 PetscObject cellgeomobj, facegeomobj; 2348 PetscErrorCode ierr; 2349 2350 PetscFunctionBegin; 2351 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2352 if (!cellgeomobj) { 2353 Vec cellgeomInt, facegeomInt; 2354 2355 ierr = DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);CHKERRQ(ierr); 2356 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);CHKERRQ(ierr); 2357 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);CHKERRQ(ierr); 2358 ierr = VecDestroy(&cellgeomInt);CHKERRQ(ierr); 2359 ierr = VecDestroy(&facegeomInt);CHKERRQ(ierr); 2360 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2361 } 2362 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);CHKERRQ(ierr); 2363 if (cellgeom) *cellgeom = (Vec) cellgeomobj; 2364 if (facegeom) *facegeom = (Vec) facegeomobj; 2365 if (gradDM) { 2366 PetscObject gradobj; 2367 PetscBool computeGradients; 2368 2369 ierr = PetscFVGetComputeGradients(fv,&computeGradients);CHKERRQ(ierr); 2370 if (!computeGradients) { 2371 *gradDM = NULL; 2372 PetscFunctionReturn(0); 2373 } 2374 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2375 if (!gradobj) { 2376 DM dmGradInt; 2377 2378 ierr = DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);CHKERRQ(ierr); 2379 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);CHKERRQ(ierr); 2380 ierr = DMDestroy(&dmGradInt);CHKERRQ(ierr); 2381 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2382 } 2383 *gradDM = (DM) gradobj; 2384 } 2385 PetscFunctionReturn(0); 2386 } 2387 2388 static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess) 2389 { 2390 PetscInt l, m; 2391 2392 PetscFunctionBeginHot; 2393 if (dimC == dimR && dimR <= 3) { 2394 /* invert Jacobian, multiply */ 2395 PetscScalar det, idet; 2396 2397 switch (dimR) { 2398 case 1: 2399 invJ[0] = 1./ J[0]; 2400 break; 2401 case 2: 2402 det = J[0] * J[3] - J[1] * J[2]; 2403 idet = 1./det; 2404 invJ[0] = J[3] * idet; 2405 invJ[1] = -J[1] * idet; 2406 invJ[2] = -J[2] * idet; 2407 invJ[3] = J[0] * idet; 2408 break; 2409 case 3: 2410 { 2411 invJ[0] = J[4] * J[8] - J[5] * J[7]; 2412 invJ[1] = J[2] * J[7] - J[1] * J[8]; 2413 invJ[2] = J[1] * J[5] - J[2] * J[4]; 2414 det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6]; 2415 idet = 1./det; 2416 invJ[0] *= idet; 2417 invJ[1] *= idet; 2418 invJ[2] *= idet; 2419 invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]); 2420 invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]); 2421 invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]); 2422 invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]); 2423 invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]); 2424 invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]); 2425 } 2426 break; 2427 } 2428 for (l = 0; l < dimR; l++) { 2429 for (m = 0; m < dimC; m++) { 2430 guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m]; 2431 } 2432 } 2433 } else { 2434 #if defined(PETSC_USE_COMPLEX) 2435 char transpose = 'C'; 2436 #else 2437 char transpose = 'T'; 2438 #endif 2439 PetscBLASInt m = dimR; 2440 PetscBLASInt n = dimC; 2441 PetscBLASInt one = 1; 2442 PetscBLASInt worksize = dimR * dimC, info; 2443 2444 for (l = 0; l < dimC; l++) {invJ[l] = resNeg[l];} 2445 2446 PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&m,&n,&one,J,&m,invJ,&n,work,&worksize, &info)); 2447 if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS"); 2448 2449 for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);} 2450 } 2451 PetscFunctionReturn(0); 2452 } 2453 2454 static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2455 { 2456 PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR); 2457 PetscScalar *coordsScalar = NULL; 2458 PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg; 2459 PetscScalar *J, *invJ, *work; 2460 PetscErrorCode ierr; 2461 2462 PetscFunctionBegin; 2463 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2464 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2465 if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);CHKERRQ(ierr); 2466 ierr = DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, PETSC_REAL, &cellData);CHKERRQ(ierr); 2467 ierr = DMGetWorkArray(dm, 3 * dimR * dimC, PETSC_SCALAR, &J);CHKERRQ(ierr); 2468 cellCoords = &cellData[0]; 2469 cellCoeffs = &cellData[coordSize]; 2470 extJ = &cellData[2 * coordSize]; 2471 resNeg = &cellData[2 * coordSize + dimR]; 2472 invJ = &J[dimR * dimC]; 2473 work = &J[2 * dimR * dimC]; 2474 if (dimR == 2) { 2475 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 2476 2477 for (i = 0; i < 4; i++) { 2478 PetscInt plexI = zToPlex[i]; 2479 2480 for (j = 0; j < dimC; j++) { 2481 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2482 } 2483 } 2484 } else if (dimR == 3) { 2485 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2486 2487 for (i = 0; i < 8; i++) { 2488 PetscInt plexI = zToPlex[i]; 2489 2490 for (j = 0; j < dimC; j++) { 2491 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2492 } 2493 } 2494 } else { 2495 for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);} 2496 } 2497 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 2498 for (i = 0; i < dimR; i++) { 2499 PetscReal *swap; 2500 2501 for (j = 0; j < (numV / 2); j++) { 2502 for (k = 0; k < dimC; k++) { 2503 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 2504 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 2505 } 2506 } 2507 2508 if (i < dimR - 1) { 2509 swap = cellCoeffs; 2510 cellCoeffs = cellCoords; 2511 cellCoords = swap; 2512 } 2513 } 2514 ierr = PetscMemzero(refCoords,numPoints * dimR * sizeof (PetscReal));CHKERRQ(ierr); 2515 for (j = 0; j < numPoints; j++) { 2516 for (i = 0; i < maxIts; i++) { 2517 PetscReal *guess = &refCoords[dimR * j]; 2518 2519 /* compute -residual and Jacobian */ 2520 for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];} 2521 for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;} 2522 for (k = 0; k < numV; k++) { 2523 PetscReal extCoord = 1.; 2524 for (l = 0; l < dimR; l++) { 2525 PetscReal coord = guess[l]; 2526 PetscInt dep = (k & (1 << l)) >> l; 2527 2528 extCoord *= dep * coord + !dep; 2529 extJ[l] = dep; 2530 2531 for (m = 0; m < dimR; m++) { 2532 PetscReal coord = guess[m]; 2533 PetscInt dep = ((k & (1 << m)) >> m) && (m != l); 2534 PetscReal mult = dep * coord + !dep; 2535 2536 extJ[l] *= mult; 2537 } 2538 } 2539 for (l = 0; l < dimC; l++) { 2540 PetscReal coeff = cellCoeffs[dimC * k + l]; 2541 2542 resNeg[l] -= coeff * extCoord; 2543 for (m = 0; m < dimR; m++) { 2544 J[dimR * l + m] += coeff * extJ[m]; 2545 } 2546 } 2547 } 2548 #if 0 && defined(PETSC_USE_DEBUG) 2549 { 2550 PetscReal maxAbs = 0.; 2551 2552 for (l = 0; l < dimC; l++) { 2553 maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l])); 2554 } 2555 ierr = PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,maxAbs);CHKERRQ(ierr); 2556 } 2557 #endif 2558 2559 ierr = DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);CHKERRQ(ierr); 2560 } 2561 } 2562 ierr = DMRestoreWorkArray(dm, 3 * dimR * dimC, PETSC_SCALAR, &J);CHKERRQ(ierr); 2563 ierr = DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, PETSC_REAL, &cellData);CHKERRQ(ierr); 2564 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2565 PetscFunctionReturn(0); 2566 } 2567 2568 static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2569 { 2570 PetscInt coordSize, i, j, k, l, numV = (1 << dimR); 2571 PetscScalar *coordsScalar = NULL; 2572 PetscReal *cellData, *cellCoords, *cellCoeffs; 2573 PetscErrorCode ierr; 2574 2575 PetscFunctionBegin; 2576 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2577 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2578 if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);CHKERRQ(ierr); 2579 ierr = DMGetWorkArray(dm, 2 * coordSize, PETSC_REAL, &cellData);CHKERRQ(ierr); 2580 cellCoords = &cellData[0]; 2581 cellCoeffs = &cellData[coordSize]; 2582 if (dimR == 2) { 2583 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 2584 2585 for (i = 0; i < 4; i++) { 2586 PetscInt plexI = zToPlex[i]; 2587 2588 for (j = 0; j < dimC; j++) { 2589 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2590 } 2591 } 2592 } else if (dimR == 3) { 2593 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2594 2595 for (i = 0; i < 8; i++) { 2596 PetscInt plexI = zToPlex[i]; 2597 2598 for (j = 0; j < dimC; j++) { 2599 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2600 } 2601 } 2602 } else { 2603 for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);} 2604 } 2605 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 2606 for (i = 0; i < dimR; i++) { 2607 PetscReal *swap; 2608 2609 for (j = 0; j < (numV / 2); j++) { 2610 for (k = 0; k < dimC; k++) { 2611 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 2612 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 2613 } 2614 } 2615 2616 if (i < dimR - 1) { 2617 swap = cellCoeffs; 2618 cellCoeffs = cellCoords; 2619 cellCoords = swap; 2620 } 2621 } 2622 ierr = PetscMemzero(realCoords,numPoints * dimC * sizeof (PetscReal));CHKERRQ(ierr); 2623 for (j = 0; j < numPoints; j++) { 2624 const PetscReal *guess = &refCoords[dimR * j]; 2625 PetscReal *mapped = &realCoords[dimC * j]; 2626 2627 for (k = 0; k < numV; k++) { 2628 PetscReal extCoord = 1.; 2629 for (l = 0; l < dimR; l++) { 2630 PetscReal coord = guess[l]; 2631 PetscInt dep = (k & (1 << l)) >> l; 2632 2633 extCoord *= dep * coord + !dep; 2634 } 2635 for (l = 0; l < dimC; l++) { 2636 PetscReal coeff = cellCoeffs[dimC * k + l]; 2637 2638 mapped[l] += coeff * extCoord; 2639 } 2640 } 2641 } 2642 ierr = DMRestoreWorkArray(dm, 2 * coordSize, PETSC_REAL, &cellData);CHKERRQ(ierr); 2643 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2644 PetscFunctionReturn(0); 2645 } 2646 2647 static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2648 { 2649 PetscInt numComp, numDof, i, j, k, l, m, maxIter = 7, coordSize; 2650 PetscScalar *nodes = NULL; 2651 PetscReal *invV, *modes; 2652 PetscReal *B, *D, *resNeg; 2653 PetscScalar *J, *invJ, *work; 2654 PetscErrorCode ierr; 2655 2656 PetscFunctionBegin; 2657 ierr = PetscFEGetDimension(fe, &numDof);CHKERRQ(ierr); 2658 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 2659 if (numComp != dimC) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,dimC); 2660 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2661 /* convert nodes to values in the stable evaluation basis */ 2662 ierr = DMGetWorkArray(dm,dimC * numDof,PETSC_REAL,&modes);CHKERRQ(ierr); 2663 invV = fe->invV; 2664 for (i = 0; i < numDof; i++) { 2665 for (j = 0; j < dimC; j++) { 2666 modes[i * dimC + j] = 0.; 2667 for (k = 0; k < numDof; k++) { 2668 modes[i * dimC + j] += invV[i * numDof + k] * PetscRealPart(nodes[k * dimC + j]); 2669 } 2670 } 2671 } 2672 ierr = DMGetWorkArray(dm,numDof + numDof * dimR + dimC,PETSC_REAL,&B);CHKERRQ(ierr); 2673 D = &B[numDof]; 2674 resNeg = &D[numDof * dimR]; 2675 ierr = DMGetWorkArray(dm,3 * dimC * dimR,PETSC_SCALAR,&J);CHKERRQ(ierr); 2676 invJ = &J[dimC * dimR]; 2677 work = &invJ[dimC * dimR]; 2678 for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;} 2679 for (j = 0; j < numPoints; j++) { 2680 for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */ 2681 PetscReal *guess = &refCoords[j * dimR]; 2682 ierr = PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);CHKERRQ(ierr); 2683 for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[j * dimC + k];} 2684 for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;} 2685 for (k = 0; k < numDof; k++) { 2686 for (l = 0; l < dimC; l++) { 2687 resNeg[l] -= modes[k * dimC + l] * B[k]; 2688 for (m = 0; m < dimR; m++) { 2689 J[l * dimR + m] += modes[k * dimC + l] * D[k * dimR + m]; 2690 } 2691 } 2692 } 2693 #if 0 && defined(PETSC_USE_DEBUG) 2694 { 2695 PetscReal maxAbs = 0.; 2696 2697 for (l = 0; l < dimC; l++) { 2698 maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l])); 2699 } 2700 ierr = PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,maxAbs);CHKERRQ(ierr); 2701 } 2702 #endif 2703 ierr = DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);CHKERRQ(ierr); 2704 } 2705 } 2706 ierr = DMRestoreWorkArray(dm,3 * dimC * dimR,PETSC_SCALAR,&J);CHKERRQ(ierr); 2707 ierr = DMRestoreWorkArray(dm,numDof + numDof * dimR + dimC,PETSC_REAL,&B);CHKERRQ(ierr); 2708 ierr = DMRestoreWorkArray(dm,dimC * numDof,PETSC_REAL,&modes);CHKERRQ(ierr); 2709 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2710 PetscFunctionReturn(0); 2711 } 2712 2713 static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2714 { 2715 PetscInt numComp, numDof, i, j, k, l, coordSize; 2716 PetscScalar *nodes = NULL; 2717 PetscReal *invV, *modes; 2718 PetscReal *B; 2719 PetscErrorCode ierr; 2720 2721 PetscFunctionBegin; 2722 ierr = PetscFEGetDimension(fe, &numDof);CHKERRQ(ierr); 2723 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 2724 if (numComp != dimC) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,dimC); 2725 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2726 /* convert nodes to values in the stable evaluation basis */ 2727 ierr = DMGetWorkArray(dm,dimC * numDof,PETSC_REAL,&modes);CHKERRQ(ierr); 2728 invV = fe->invV; 2729 for (i = 0; i < numDof; i++) { 2730 for (j = 0; j < dimC; j++) { 2731 modes[i * dimC + j] = 0.; 2732 for (k = 0; k < numDof; k++) { 2733 modes[i * dimC + j] += invV[i * numDof + k] * PetscRealPart(nodes[k * dimC + j]); 2734 } 2735 } 2736 } 2737 ierr = DMGetWorkArray(dm,numDof,PETSC_REAL,&B);CHKERRQ(ierr); 2738 for (i = 0; i < numPoints * dimC; i++) {realCoords[i] = 0.;} 2739 for (j = 0; j < numPoints; j++) { 2740 const PetscReal *guess = &refCoords[j * dimR]; 2741 PetscReal *mapped = &realCoords[j * dimC]; 2742 2743 ierr = PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, NULL, NULL);CHKERRQ(ierr); 2744 for (k = 0; k < numDof; k++) { 2745 for (l = 0; l < dimC; l++) { 2746 mapped[l] += modes[k * dimC + l] * B[k]; 2747 } 2748 } 2749 } 2750 ierr = DMRestoreWorkArray(dm,numDof,PETSC_REAL,&B);CHKERRQ(ierr); 2751 ierr = DMRestoreWorkArray(dm,dimC * numDof,PETSC_REAL,&modes);CHKERRQ(ierr); 2752 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2753 PetscFunctionReturn(0); 2754 } 2755 2756 /*@ 2757 DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element 2758 map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not 2759 extend uniquely outside the reference cell (e.g, most non-affine maps) 2760 2761 Not collective 2762 2763 Input Parameters: 2764 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 2765 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 2766 as a multilinear map for tensor-product elements 2767 . cell - the cell whose map is used. 2768 . numPoints - the number of points to locate 2769 - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 2770 2771 Output Parameters: 2772 . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 2773 2774 Level: intermediate 2775 @*/ 2776 PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[]) 2777 { 2778 PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i; 2779 DM coordDM = NULL; 2780 Vec coords; 2781 PetscFE fe = NULL; 2782 PetscErrorCode ierr; 2783 2784 PetscFunctionBegin; 2785 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2786 ierr = DMGetDimension(dm,&dimR);CHKERRQ(ierr); 2787 ierr = DMGetCoordinateDim(dm,&dimC);CHKERRQ(ierr); 2788 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 2789 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 2790 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2791 ierr = DMGetCoordinateDM(dm,&coordDM);CHKERRQ(ierr); 2792 if (coordDM) { 2793 PetscInt coordFields; 2794 2795 ierr = DMGetNumFields(coordDM,&coordFields);CHKERRQ(ierr); 2796 if (coordFields) { 2797 PetscClassId id; 2798 PetscObject disc; 2799 2800 ierr = DMGetField(coordDM,0,&disc);CHKERRQ(ierr); 2801 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 2802 if (id == PETSCFE_CLASSID) { 2803 fe = (PetscFE) disc; 2804 } 2805 } 2806 } 2807 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 2808 ierr = DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 2809 cEnd = cEndInterior > 0 ? cEndInterior : cEnd; 2810 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); 2811 if (!fe) { /* implicit discretization: affine or multilinear */ 2812 PetscInt coneSize; 2813 PetscBool isSimplex, isTensor; 2814 2815 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 2816 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 2817 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 2818 if (isSimplex) { 2819 PetscReal detJ, *v0, *J, *invJ; 2820 2821 ierr = DMGetWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2822 J = &v0[dimC]; 2823 invJ = &J[dimC * dimC]; 2824 ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);CHKERRQ(ierr); 2825 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 2826 CoordinatesRealToRef(dimC, dimR, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]); 2827 } 2828 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2829 } else if (isTensor) { 2830 ierr = DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);CHKERRQ(ierr); 2831 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize); 2832 } else { 2833 ierr = DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);CHKERRQ(ierr); 2834 } 2835 PetscFunctionReturn(0); 2836 } 2837 2838 /*@ 2839 DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map. 2840 2841 Not collective 2842 2843 Input Parameters: 2844 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 2845 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 2846 as a multilinear map for tensor-product elements 2847 . cell - the cell whose map is used. 2848 . numPoints - the number of points to locate 2849 + refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 2850 2851 Output Parameters: 2852 . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 2853 2854 Level: intermediate 2855 @*/ 2856 PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[]) 2857 { 2858 PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i; 2859 DM coordDM = NULL; 2860 Vec coords; 2861 PetscFE fe = NULL; 2862 PetscErrorCode ierr; 2863 2864 PetscFunctionBegin; 2865 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2866 ierr = DMGetDimension(dm,&dimR);CHKERRQ(ierr); 2867 ierr = DMGetCoordinateDim(dm,&dimC);CHKERRQ(ierr); 2868 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 2869 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 2870 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2871 ierr = DMGetCoordinateDM(dm,&coordDM);CHKERRQ(ierr); 2872 if (coordDM) { 2873 PetscInt coordFields; 2874 2875 ierr = DMGetNumFields(coordDM,&coordFields);CHKERRQ(ierr); 2876 if (coordFields) { 2877 PetscClassId id; 2878 PetscObject disc; 2879 2880 ierr = DMGetField(coordDM,0,&disc);CHKERRQ(ierr); 2881 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 2882 if (id == PETSCFE_CLASSID) { 2883 fe = (PetscFE) disc; 2884 } 2885 } 2886 } 2887 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 2888 ierr = DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 2889 cEnd = cEndInterior > 0 ? cEndInterior : cEnd; 2890 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); 2891 if (!fe) { /* implicit discretization: affine or multilinear */ 2892 PetscInt coneSize; 2893 PetscBool isSimplex, isTensor; 2894 2895 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 2896 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 2897 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 2898 if (isSimplex) { 2899 PetscReal detJ, *v0, *J; 2900 2901 ierr = DMGetWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2902 J = &v0[dimC]; 2903 ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ);CHKERRQ(ierr); 2904 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 2905 CoordinatesRefToReal(dimC, dimR, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]); 2906 } 2907 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2908 } else if (isTensor) { 2909 ierr = DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);CHKERRQ(ierr); 2910 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize); 2911 } else { 2912 ierr = DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);CHKERRQ(ierr); 2913 } 2914 PetscFunctionReturn(0); 2915 } 2916