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 /* This should also take a PetscFE argument I think */ 1847 PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom) 1848 { 1849 DM dmCell; 1850 Vec coordinates; 1851 PetscSection coordSection, sectionCell; 1852 PetscScalar *cgeom; 1853 PetscInt cStart, cEnd, cMax, c; 1854 PetscErrorCode ierr; 1855 1856 PetscFunctionBegin; 1857 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1858 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1859 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1860 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1861 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1862 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1863 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1864 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1865 cEnd = cMax < 0 ? cEnd : cMax; 1866 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1867 /* TODO This needs to be multiplied by Nq for non-affine */ 1868 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFECellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1869 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1870 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1871 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1872 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1873 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1874 for (c = cStart; c < cEnd; ++c) { 1875 PetscFECellGeom *cg; 1876 1877 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1878 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1879 ierr = DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v0, cg->J, cg->invJ, &cg->detJ);CHKERRQ(ierr); 1880 if (cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c); 1881 } 1882 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1883 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1884 PetscFunctionReturn(0); 1885 } 1886 1887 /*@ 1888 DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method 1889 1890 Input Parameter: 1891 . dm - The DM 1892 1893 Output Parameters: 1894 + cellgeom - A Vec of PetscFVCellGeom data 1895 . facegeom - A Vec of PetscFVFaceGeom data 1896 1897 Level: developer 1898 1899 .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM() 1900 @*/ 1901 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 1902 { 1903 DM dmFace, dmCell; 1904 DMLabel ghostLabel; 1905 PetscSection sectionFace, sectionCell; 1906 PetscSection coordSection; 1907 Vec coordinates; 1908 PetscScalar *fgeom, *cgeom; 1909 PetscReal minradius, gminradius; 1910 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 1911 PetscErrorCode ierr; 1912 1913 PetscFunctionBegin; 1914 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1915 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1916 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1917 /* Make cell centroids and volumes */ 1918 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1919 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1920 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1921 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1922 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1923 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1924 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1925 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1926 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1927 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1928 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1929 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1930 if (cEndInterior < 0) { 1931 cEndInterior = cEnd; 1932 } 1933 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1934 for (c = cStart; c < cEndInterior; ++c) { 1935 PetscFVCellGeom *cg; 1936 1937 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1938 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1939 ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr); 1940 } 1941 /* Compute face normals and minimum cell radius */ 1942 ierr = DMClone(dm, &dmFace);CHKERRQ(ierr); 1943 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);CHKERRQ(ierr); 1944 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1945 ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr); 1946 for (f = fStart; f < fEnd; ++f) {ierr = PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1947 ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr); 1948 ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr); 1949 ierr = PetscSectionDestroy(§ionFace);CHKERRQ(ierr); 1950 ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr); 1951 ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr); 1952 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1953 minradius = PETSC_MAX_REAL; 1954 for (f = fStart; f < fEnd; ++f) { 1955 PetscFVFaceGeom *fg; 1956 PetscReal area; 1957 PetscInt ghost = -1, d, numChildren; 1958 1959 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1960 ierr = DMPlexGetTreeChildren(dm,f,&numChildren,NULL);CHKERRQ(ierr); 1961 if (ghost >= 0 || numChildren) continue; 1962 ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr); 1963 ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);CHKERRQ(ierr); 1964 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 1965 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 1966 { 1967 PetscFVCellGeom *cL, *cR; 1968 PetscInt ncells; 1969 const PetscInt *cells; 1970 PetscReal *lcentroid, *rcentroid; 1971 PetscReal l[3], r[3], v[3]; 1972 1973 ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr); 1974 ierr = DMPlexGetSupportSize(dm, f, &ncells);CHKERRQ(ierr); 1975 ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr); 1976 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 1977 if (ncells > 1) { 1978 ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr); 1979 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 1980 } 1981 else { 1982 rcentroid = fg->centroid; 1983 } 1984 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);CHKERRQ(ierr); 1985 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);CHKERRQ(ierr); 1986 DMPlex_WaxpyD_Internal(dim, -1, l, r, v); 1987 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 1988 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 1989 } 1990 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 1991 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]); 1992 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]); 1993 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f); 1994 } 1995 if (cells[0] < cEndInterior) { 1996 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 1997 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1998 } 1999 if (ncells > 1 && cells[1] < cEndInterior) { 2000 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 2001 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2002 } 2003 } 2004 } 2005 ierr = MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 2006 ierr = DMPlexSetMinRadius(dm, gminradius);CHKERRQ(ierr); 2007 /* Compute centroids of ghost cells */ 2008 for (c = cEndInterior; c < cEnd; ++c) { 2009 PetscFVFaceGeom *fg; 2010 const PetscInt *cone, *support; 2011 PetscInt coneSize, supportSize, s; 2012 2013 ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr); 2014 if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize); 2015 ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr); 2016 ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr); 2017 if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize); 2018 ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr); 2019 ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr); 2020 for (s = 0; s < 2; ++s) { 2021 /* Reflect ghost centroid across plane of face */ 2022 if (support[s] == c) { 2023 PetscFVCellGeom *ci; 2024 PetscFVCellGeom *cg; 2025 PetscReal c2f[3], a; 2026 2027 ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr); 2028 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 2029 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 2030 ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr); 2031 DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid); 2032 cg->volume = ci->volume; 2033 } 2034 } 2035 } 2036 ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr); 2037 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 2038 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 2039 ierr = DMDestroy(&dmFace);CHKERRQ(ierr); 2040 PetscFunctionReturn(0); 2041 } 2042 2043 /*@C 2044 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 2045 2046 Not collective 2047 2048 Input Argument: 2049 . dm - the DM 2050 2051 Output Argument: 2052 . minradius - the minium cell radius 2053 2054 Level: developer 2055 2056 .seealso: DMGetCoordinates() 2057 @*/ 2058 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 2059 { 2060 PetscFunctionBegin; 2061 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2062 PetscValidPointer(minradius,2); 2063 *minradius = ((DM_Plex*) dm->data)->minradius; 2064 PetscFunctionReturn(0); 2065 } 2066 2067 /*@C 2068 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 2069 2070 Logically collective 2071 2072 Input Arguments: 2073 + dm - the DM 2074 - minradius - the minium cell radius 2075 2076 Level: developer 2077 2078 .seealso: DMSetCoordinates() 2079 @*/ 2080 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 2081 { 2082 PetscFunctionBegin; 2083 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2084 ((DM_Plex*) dm->data)->minradius = minradius; 2085 PetscFunctionReturn(0); 2086 } 2087 2088 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 2089 { 2090 DMLabel ghostLabel; 2091 PetscScalar *dx, *grad, **gref; 2092 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 2093 PetscErrorCode ierr; 2094 2095 PetscFunctionBegin; 2096 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2097 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2098 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2099 ierr = DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);CHKERRQ(ierr); 2100 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 2101 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 2102 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 2103 for (c = cStart; c < cEndInterior; c++) { 2104 const PetscInt *faces; 2105 PetscInt numFaces, usedFaces, f, d; 2106 PetscFVCellGeom *cg; 2107 PetscBool boundary; 2108 PetscInt ghost; 2109 2110 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2111 ierr = DMPlexGetConeSize(dm, c, &numFaces);CHKERRQ(ierr); 2112 ierr = DMPlexGetCone(dm, c, &faces);CHKERRQ(ierr); 2113 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 2114 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 2115 PetscFVCellGeom *cg1; 2116 PetscFVFaceGeom *fg; 2117 const PetscInt *fcells; 2118 PetscInt ncell, side; 2119 2120 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 2121 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 2122 if ((ghost >= 0) || boundary) continue; 2123 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr); 2124 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 2125 ncell = fcells[!side]; /* the neighbor */ 2126 ierr = DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr); 2127 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 2128 for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d]; 2129 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2130 } 2131 if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 2132 ierr = PetscFVComputeGradient(fvm, usedFaces, dx, grad);CHKERRQ(ierr); 2133 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 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 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d]; 2138 ++usedFaces; 2139 } 2140 } 2141 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 2142 PetscFunctionReturn(0); 2143 } 2144 2145 static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 2146 { 2147 DMLabel ghostLabel; 2148 PetscScalar *dx, *grad, **gref; 2149 PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0; 2150 PetscSection neighSec; 2151 PetscInt (*neighbors)[2]; 2152 PetscInt *counter; 2153 PetscErrorCode ierr; 2154 2155 PetscFunctionBegin; 2156 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2157 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2158 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2159 if (cEndInterior < 0) { 2160 cEndInterior = cEnd; 2161 } 2162 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);CHKERRQ(ierr); 2163 ierr = PetscSectionSetChart(neighSec,cStart,cEndInterior);CHKERRQ(ierr); 2164 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2165 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 2166 for (f = fStart; f < fEnd; f++) { 2167 const PetscInt *fcells; 2168 PetscBool boundary; 2169 PetscInt ghost = -1; 2170 PetscInt numChildren, numCells, c; 2171 2172 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 2173 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 2174 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 2175 if ((ghost >= 0) || boundary || numChildren) continue; 2176 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 2177 if (numCells == 2) { 2178 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 2179 for (c = 0; c < 2; c++) { 2180 PetscInt cell = fcells[c]; 2181 2182 if (cell >= cStart && cell < cEndInterior) { 2183 ierr = PetscSectionAddDof(neighSec,cell,1);CHKERRQ(ierr); 2184 } 2185 } 2186 } 2187 } 2188 ierr = PetscSectionSetUp(neighSec);CHKERRQ(ierr); 2189 ierr = PetscSectionGetMaxDof(neighSec,&maxNumFaces);CHKERRQ(ierr); 2190 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 2191 nStart = 0; 2192 ierr = PetscSectionGetStorageSize(neighSec,&nEnd);CHKERRQ(ierr); 2193 ierr = PetscMalloc1((nEnd-nStart),&neighbors);CHKERRQ(ierr); 2194 ierr = PetscCalloc1((cEndInterior-cStart),&counter);CHKERRQ(ierr); 2195 for (f = fStart; f < fEnd; f++) { 2196 const PetscInt *fcells; 2197 PetscBool boundary; 2198 PetscInt ghost = -1; 2199 PetscInt numChildren, numCells, c; 2200 2201 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 2202 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 2203 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 2204 if ((ghost >= 0) || boundary || numChildren) continue; 2205 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 2206 if (numCells == 2) { 2207 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 2208 for (c = 0; c < 2; c++) { 2209 PetscInt cell = fcells[c], off; 2210 2211 if (cell >= cStart && cell < cEndInterior) { 2212 ierr = PetscSectionGetOffset(neighSec,cell,&off);CHKERRQ(ierr); 2213 off += counter[cell - cStart]++; 2214 neighbors[off][0] = f; 2215 neighbors[off][1] = fcells[1 - c]; 2216 } 2217 } 2218 } 2219 } 2220 ierr = PetscFree(counter);CHKERRQ(ierr); 2221 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 2222 for (c = cStart; c < cEndInterior; c++) { 2223 PetscInt numFaces, f, d, off, ghost = -1; 2224 PetscFVCellGeom *cg; 2225 2226 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2227 ierr = PetscSectionGetDof(neighSec, c, &numFaces);CHKERRQ(ierr); 2228 ierr = PetscSectionGetOffset(neighSec, c, &off);CHKERRQ(ierr); 2229 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, c, &ghost);CHKERRQ(ierr);} 2230 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); 2231 for (f = 0; f < numFaces; ++f) { 2232 PetscFVCellGeom *cg1; 2233 PetscFVFaceGeom *fg; 2234 const PetscInt *fcells; 2235 PetscInt ncell, side, nface; 2236 2237 nface = neighbors[off + f][0]; 2238 ncell = neighbors[off + f][1]; 2239 ierr = DMPlexGetSupport(dm,nface,&fcells);CHKERRQ(ierr); 2240 side = (c != fcells[0]); 2241 ierr = DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);CHKERRQ(ierr); 2242 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 2243 for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d]; 2244 gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2245 } 2246 ierr = PetscFVComputeGradient(fvm, numFaces, dx, grad);CHKERRQ(ierr); 2247 for (f = 0; f < numFaces; ++f) { 2248 for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d]; 2249 } 2250 } 2251 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 2252 ierr = PetscSectionDestroy(&neighSec);CHKERRQ(ierr); 2253 ierr = PetscFree(neighbors);CHKERRQ(ierr); 2254 PetscFunctionReturn(0); 2255 } 2256 2257 /*@ 2258 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 2259 2260 Collective on DM 2261 2262 Input Arguments: 2263 + dm - The DM 2264 . fvm - The PetscFV 2265 . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM() 2266 - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM() 2267 2268 Output Parameters: 2269 + faceGeometry - The geometric factors for gradient calculation are inserted 2270 - dmGrad - The DM describing the layout of gradient data 2271 2272 Level: developer 2273 2274 .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM() 2275 @*/ 2276 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 2277 { 2278 DM dmFace, dmCell; 2279 PetscScalar *fgeom, *cgeom; 2280 PetscSection sectionGrad, parentSection; 2281 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 2282 PetscErrorCode ierr; 2283 2284 PetscFunctionBegin; 2285 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2286 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 2287 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2288 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2289 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 2290 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 2291 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 2292 ierr = VecGetArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2293 ierr = VecGetArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2294 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 2295 if (!parentSection) { 2296 ierr = BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2297 } else { 2298 ierr = BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2299 } 2300 ierr = VecRestoreArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2301 ierr = VecRestoreArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2302 /* Create storage for gradients */ 2303 ierr = DMClone(dm, dmGrad);CHKERRQ(ierr); 2304 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 2305 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 2306 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 2307 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 2308 ierr = DMSetDefaultSection(*dmGrad, sectionGrad);CHKERRQ(ierr); 2309 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 2310 PetscFunctionReturn(0); 2311 } 2312 2313 PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM) 2314 { 2315 PetscObject cellgeomobj, facegeomobj; 2316 PetscErrorCode ierr; 2317 2318 PetscFunctionBegin; 2319 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2320 if (!cellgeomobj) { 2321 Vec cellgeomInt, facegeomInt; 2322 2323 ierr = DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);CHKERRQ(ierr); 2324 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);CHKERRQ(ierr); 2325 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);CHKERRQ(ierr); 2326 ierr = VecDestroy(&cellgeomInt);CHKERRQ(ierr); 2327 ierr = VecDestroy(&facegeomInt);CHKERRQ(ierr); 2328 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2329 } 2330 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);CHKERRQ(ierr); 2331 if (cellgeom) *cellgeom = (Vec) cellgeomobj; 2332 if (facegeom) *facegeom = (Vec) facegeomobj; 2333 if (gradDM) { 2334 PetscObject gradobj; 2335 PetscBool computeGradients; 2336 2337 ierr = PetscFVGetComputeGradients(fv,&computeGradients);CHKERRQ(ierr); 2338 if (!computeGradients) { 2339 *gradDM = NULL; 2340 PetscFunctionReturn(0); 2341 } 2342 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2343 if (!gradobj) { 2344 DM dmGradInt; 2345 2346 ierr = DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);CHKERRQ(ierr); 2347 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);CHKERRQ(ierr); 2348 ierr = DMDestroy(&dmGradInt);CHKERRQ(ierr); 2349 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2350 } 2351 *gradDM = (DM) gradobj; 2352 } 2353 PetscFunctionReturn(0); 2354 } 2355 2356 static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess) 2357 { 2358 PetscInt l, m; 2359 2360 PetscFunctionBeginHot; 2361 if (dimC == dimR && dimR <= 3) { 2362 /* invert Jacobian, multiply */ 2363 PetscScalar det, idet; 2364 2365 switch (dimR) { 2366 case 1: 2367 invJ[0] = 1./ J[0]; 2368 break; 2369 case 2: 2370 det = J[0] * J[3] - J[1] * J[2]; 2371 idet = 1./det; 2372 invJ[0] = J[3] * idet; 2373 invJ[1] = -J[1] * idet; 2374 invJ[2] = -J[2] * idet; 2375 invJ[3] = J[0] * idet; 2376 break; 2377 case 3: 2378 { 2379 invJ[0] = J[4] * J[8] - J[5] * J[7]; 2380 invJ[1] = J[2] * J[7] - J[1] * J[8]; 2381 invJ[2] = J[1] * J[5] - J[2] * J[4]; 2382 det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6]; 2383 idet = 1./det; 2384 invJ[0] *= idet; 2385 invJ[1] *= idet; 2386 invJ[2] *= idet; 2387 invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]); 2388 invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]); 2389 invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]); 2390 invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]); 2391 invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]); 2392 invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]); 2393 } 2394 break; 2395 } 2396 for (l = 0; l < dimR; l++) { 2397 for (m = 0; m < dimC; m++) { 2398 guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m]; 2399 } 2400 } 2401 } else { 2402 #if defined(PETSC_USE_COMPLEX) 2403 char transpose = 'C'; 2404 #else 2405 char transpose = 'T'; 2406 #endif 2407 PetscBLASInt m = dimR; 2408 PetscBLASInt n = dimC; 2409 PetscBLASInt one = 1; 2410 PetscBLASInt worksize = dimR * dimC, info; 2411 2412 for (l = 0; l < dimC; l++) {invJ[l] = resNeg[l];} 2413 2414 PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&m,&n,&one,J,&m,invJ,&n,work,&worksize, &info)); 2415 if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS"); 2416 2417 for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);} 2418 } 2419 PetscFunctionReturn(0); 2420 } 2421 2422 static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2423 { 2424 PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR); 2425 PetscScalar *coordsScalar = NULL; 2426 PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg; 2427 PetscScalar *J, *invJ, *work; 2428 PetscErrorCode ierr; 2429 2430 PetscFunctionBegin; 2431 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2432 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2433 if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);CHKERRQ(ierr); 2434 ierr = DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, PETSC_REAL, &cellData);CHKERRQ(ierr); 2435 ierr = DMGetWorkArray(dm, 3 * dimR * dimC, PETSC_SCALAR, &J);CHKERRQ(ierr); 2436 cellCoords = &cellData[0]; 2437 cellCoeffs = &cellData[coordSize]; 2438 extJ = &cellData[2 * coordSize]; 2439 resNeg = &cellData[2 * coordSize + dimR]; 2440 invJ = &J[dimR * dimC]; 2441 work = &J[2 * dimR * dimC]; 2442 if (dimR == 2) { 2443 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 2444 2445 for (i = 0; i < 4; i++) { 2446 PetscInt plexI = zToPlex[i]; 2447 2448 for (j = 0; j < dimC; j++) { 2449 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2450 } 2451 } 2452 } else if (dimR == 3) { 2453 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2454 2455 for (i = 0; i < 8; i++) { 2456 PetscInt plexI = zToPlex[i]; 2457 2458 for (j = 0; j < dimC; j++) { 2459 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2460 } 2461 } 2462 } else { 2463 for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);} 2464 } 2465 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 2466 for (i = 0; i < dimR; i++) { 2467 PetscReal *swap; 2468 2469 for (j = 0; j < (numV / 2); j++) { 2470 for (k = 0; k < dimC; k++) { 2471 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 2472 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 2473 } 2474 } 2475 2476 if (i < dimR - 1) { 2477 swap = cellCoeffs; 2478 cellCoeffs = cellCoords; 2479 cellCoords = swap; 2480 } 2481 } 2482 ierr = PetscMemzero(refCoords,numPoints * dimR * sizeof (PetscReal));CHKERRQ(ierr); 2483 for (j = 0; j < numPoints; j++) { 2484 for (i = 0; i < maxIts; i++) { 2485 PetscReal *guess = &refCoords[dimR * j]; 2486 2487 /* compute -residual and Jacobian */ 2488 for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];} 2489 for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;} 2490 for (k = 0; k < numV; k++) { 2491 PetscReal extCoord = 1.; 2492 for (l = 0; l < dimR; l++) { 2493 PetscReal coord = guess[l]; 2494 PetscInt dep = (k & (1 << l)) >> l; 2495 2496 extCoord *= dep * coord + !dep; 2497 extJ[l] = dep; 2498 2499 for (m = 0; m < dimR; m++) { 2500 PetscReal coord = guess[m]; 2501 PetscInt dep = ((k & (1 << m)) >> m) && (m != l); 2502 PetscReal mult = dep * coord + !dep; 2503 2504 extJ[l] *= mult; 2505 } 2506 } 2507 for (l = 0; l < dimC; l++) { 2508 PetscReal coeff = cellCoeffs[dimC * k + l]; 2509 2510 resNeg[l] -= coeff * extCoord; 2511 for (m = 0; m < dimR; m++) { 2512 J[dimR * l + m] += coeff * extJ[m]; 2513 } 2514 } 2515 } 2516 #if 0 && defined(PETSC_USE_DEBUG) 2517 { 2518 PetscReal maxAbs = 0.; 2519 2520 for (l = 0; l < dimC; l++) { 2521 maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l])); 2522 } 2523 ierr = PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,maxAbs);CHKERRQ(ierr); 2524 } 2525 #endif 2526 2527 ierr = DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);CHKERRQ(ierr); 2528 } 2529 } 2530 ierr = DMRestoreWorkArray(dm, 3 * dimR * dimC, PETSC_SCALAR, &J);CHKERRQ(ierr); 2531 ierr = DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, PETSC_REAL, &cellData);CHKERRQ(ierr); 2532 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2533 PetscFunctionReturn(0); 2534 } 2535 2536 static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2537 { 2538 PetscInt coordSize, i, j, k, l, numV = (1 << dimR); 2539 PetscScalar *coordsScalar = NULL; 2540 PetscReal *cellData, *cellCoords, *cellCoeffs; 2541 PetscErrorCode ierr; 2542 2543 PetscFunctionBegin; 2544 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2545 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2546 if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);CHKERRQ(ierr); 2547 ierr = DMGetWorkArray(dm, 2 * coordSize, PETSC_REAL, &cellData);CHKERRQ(ierr); 2548 cellCoords = &cellData[0]; 2549 cellCoeffs = &cellData[coordSize]; 2550 if (dimR == 2) { 2551 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 2552 2553 for (i = 0; i < 4; i++) { 2554 PetscInt plexI = zToPlex[i]; 2555 2556 for (j = 0; j < dimC; j++) { 2557 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2558 } 2559 } 2560 } else if (dimR == 3) { 2561 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2562 2563 for (i = 0; i < 8; i++) { 2564 PetscInt plexI = zToPlex[i]; 2565 2566 for (j = 0; j < dimC; j++) { 2567 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2568 } 2569 } 2570 } else { 2571 for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);} 2572 } 2573 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 2574 for (i = 0; i < dimR; i++) { 2575 PetscReal *swap; 2576 2577 for (j = 0; j < (numV / 2); j++) { 2578 for (k = 0; k < dimC; k++) { 2579 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 2580 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 2581 } 2582 } 2583 2584 if (i < dimR - 1) { 2585 swap = cellCoeffs; 2586 cellCoeffs = cellCoords; 2587 cellCoords = swap; 2588 } 2589 } 2590 ierr = PetscMemzero(realCoords,numPoints * dimC * sizeof (PetscReal));CHKERRQ(ierr); 2591 for (j = 0; j < numPoints; j++) { 2592 const PetscReal *guess = &refCoords[dimR * j]; 2593 PetscReal *mapped = &realCoords[dimC * j]; 2594 2595 for (k = 0; k < numV; k++) { 2596 PetscReal extCoord = 1.; 2597 for (l = 0; l < dimR; l++) { 2598 PetscReal coord = guess[l]; 2599 PetscInt dep = (k & (1 << l)) >> l; 2600 2601 extCoord *= dep * coord + !dep; 2602 } 2603 for (l = 0; l < dimC; l++) { 2604 PetscReal coeff = cellCoeffs[dimC * k + l]; 2605 2606 mapped[l] += coeff * extCoord; 2607 } 2608 } 2609 } 2610 ierr = DMRestoreWorkArray(dm, 2 * coordSize, PETSC_REAL, &cellData);CHKERRQ(ierr); 2611 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2612 PetscFunctionReturn(0); 2613 } 2614 2615 static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2616 { 2617 PetscInt numComp, numDof, i, j, k, l, m, maxIter = 7, coordSize; 2618 PetscScalar *nodes = NULL; 2619 PetscReal *invV, *modes; 2620 PetscReal *B, *D, *resNeg; 2621 PetscScalar *J, *invJ, *work; 2622 PetscErrorCode ierr; 2623 2624 PetscFunctionBegin; 2625 ierr = PetscFEGetDimension(fe, &numDof);CHKERRQ(ierr); 2626 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 2627 if (numComp != dimC) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,dimC); 2628 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2629 /* convert nodes to values in the stable evaluation basis */ 2630 ierr = DMGetWorkArray(dm,dimC * numDof,PETSC_REAL,&modes);CHKERRQ(ierr); 2631 invV = fe->invV; 2632 for (i = 0; i < numDof; i++) { 2633 for (j = 0; j < dimC; j++) { 2634 modes[i * dimC + j] = 0.; 2635 for (k = 0; k < numDof; k++) { 2636 modes[i * dimC + j] += invV[i * numDof + k] * PetscRealPart(nodes[k * dimC + j]); 2637 } 2638 } 2639 } 2640 ierr = DMGetWorkArray(dm,numDof + numDof * dimR + dimC,PETSC_REAL,&B);CHKERRQ(ierr); 2641 D = &B[numDof]; 2642 resNeg = &D[numDof * dimR]; 2643 ierr = DMGetWorkArray(dm,3 * dimC * dimR,PETSC_SCALAR,&J);CHKERRQ(ierr); 2644 invJ = &J[dimC * dimR]; 2645 work = &invJ[dimC * dimR]; 2646 for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;} 2647 for (j = 0; j < numPoints; j++) { 2648 for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */ 2649 PetscReal *guess = &refCoords[j * dimR]; 2650 ierr = PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);CHKERRQ(ierr); 2651 for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[j * dimC + k];} 2652 for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;} 2653 for (k = 0; k < numDof; k++) { 2654 for (l = 0; l < dimC; l++) { 2655 resNeg[l] -= modes[k * dimC + l] * B[k]; 2656 for (m = 0; m < dimR; m++) { 2657 J[l * dimR + m] += modes[k * dimC + l] * D[k * dimR + m]; 2658 } 2659 } 2660 } 2661 #if 0 && defined(PETSC_USE_DEBUG) 2662 { 2663 PetscReal maxAbs = 0.; 2664 2665 for (l = 0; l < dimC; l++) { 2666 maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l])); 2667 } 2668 ierr = PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,maxAbs);CHKERRQ(ierr); 2669 } 2670 #endif 2671 ierr = DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);CHKERRQ(ierr); 2672 } 2673 } 2674 ierr = DMRestoreWorkArray(dm,3 * dimC * dimR,PETSC_SCALAR,&J);CHKERRQ(ierr); 2675 ierr = DMRestoreWorkArray(dm,numDof + numDof * dimR + dimC,PETSC_REAL,&B);CHKERRQ(ierr); 2676 ierr = DMRestoreWorkArray(dm,dimC * numDof,PETSC_REAL,&modes);CHKERRQ(ierr); 2677 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2678 PetscFunctionReturn(0); 2679 } 2680 2681 static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2682 { 2683 PetscInt numComp, numDof, i, j, k, l, coordSize; 2684 PetscScalar *nodes = NULL; 2685 PetscReal *invV, *modes; 2686 PetscReal *B; 2687 PetscErrorCode ierr; 2688 2689 PetscFunctionBegin; 2690 ierr = PetscFEGetDimension(fe, &numDof);CHKERRQ(ierr); 2691 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 2692 if (numComp != dimC) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,dimC); 2693 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2694 /* convert nodes to values in the stable evaluation basis */ 2695 ierr = DMGetWorkArray(dm,dimC * numDof,PETSC_REAL,&modes);CHKERRQ(ierr); 2696 invV = fe->invV; 2697 for (i = 0; i < numDof; i++) { 2698 for (j = 0; j < dimC; j++) { 2699 modes[i * dimC + j] = 0.; 2700 for (k = 0; k < numDof; k++) { 2701 modes[i * dimC + j] += invV[i * numDof + k] * PetscRealPart(nodes[k * dimC + j]); 2702 } 2703 } 2704 } 2705 ierr = DMGetWorkArray(dm,numDof,PETSC_REAL,&B);CHKERRQ(ierr); 2706 for (i = 0; i < numPoints * dimC; i++) {realCoords[i] = 0.;} 2707 for (j = 0; j < numPoints; j++) { 2708 const PetscReal *guess = &refCoords[j * dimR]; 2709 PetscReal *mapped = &realCoords[j * dimC]; 2710 2711 ierr = PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, NULL, NULL);CHKERRQ(ierr); 2712 for (k = 0; k < numDof; k++) { 2713 for (l = 0; l < dimC; l++) { 2714 mapped[l] += modes[k * dimC + l] * B[k]; 2715 } 2716 } 2717 } 2718 ierr = DMRestoreWorkArray(dm,numDof,PETSC_REAL,&B);CHKERRQ(ierr); 2719 ierr = DMRestoreWorkArray(dm,dimC * numDof,PETSC_REAL,&modes);CHKERRQ(ierr); 2720 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2721 PetscFunctionReturn(0); 2722 } 2723 2724 /*@ 2725 DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element 2726 map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not 2727 extend uniquely outside the reference cell (e.g, most non-affine maps) 2728 2729 Not collective 2730 2731 Input Parameters: 2732 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 2733 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 2734 as a multilinear map for tensor-product elements 2735 . cell - the cell whose map is used. 2736 . numPoints - the number of points to locate 2737 - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 2738 2739 Output Parameters: 2740 . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 2741 2742 Level: intermediate 2743 @*/ 2744 PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[]) 2745 { 2746 PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i; 2747 DM coordDM = NULL; 2748 Vec coords; 2749 PetscFE fe = NULL; 2750 PetscErrorCode ierr; 2751 2752 PetscFunctionBegin; 2753 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2754 ierr = DMGetDimension(dm,&dimR);CHKERRQ(ierr); 2755 ierr = DMGetCoordinateDim(dm,&dimC);CHKERRQ(ierr); 2756 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 2757 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 2758 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2759 ierr = DMGetCoordinateDM(dm,&coordDM);CHKERRQ(ierr); 2760 if (coordDM) { 2761 PetscInt coordFields; 2762 2763 ierr = DMGetNumFields(coordDM,&coordFields);CHKERRQ(ierr); 2764 if (coordFields) { 2765 PetscClassId id; 2766 PetscObject disc; 2767 2768 ierr = DMGetField(coordDM,0,&disc);CHKERRQ(ierr); 2769 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 2770 if (id == PETSCFE_CLASSID) { 2771 fe = (PetscFE) disc; 2772 } 2773 } 2774 } 2775 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 2776 ierr = DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 2777 cEnd = cEndInterior > 0 ? cEndInterior : cEnd; 2778 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); 2779 if (!fe) { /* implicit discretization: affine or multilinear */ 2780 PetscInt coneSize; 2781 PetscBool isSimplex, isTensor; 2782 2783 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 2784 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 2785 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 2786 if (isSimplex) { 2787 PetscReal detJ, *v0, *J, *invJ; 2788 2789 ierr = DMGetWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2790 J = &v0[dimC]; 2791 invJ = &J[dimC * dimC]; 2792 ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);CHKERRQ(ierr); 2793 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 2794 CoordinatesRealToRef(dimC, dimR, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]); 2795 } 2796 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2797 } else if (isTensor) { 2798 ierr = DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);CHKERRQ(ierr); 2799 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize); 2800 } else { 2801 ierr = DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);CHKERRQ(ierr); 2802 } 2803 PetscFunctionReturn(0); 2804 } 2805 2806 /*@ 2807 DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map. 2808 2809 Not collective 2810 2811 Input Parameters: 2812 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 2813 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 2814 as a multilinear map for tensor-product elements 2815 . cell - the cell whose map is used. 2816 . numPoints - the number of points to locate 2817 + refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 2818 2819 Output Parameters: 2820 . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 2821 2822 Level: intermediate 2823 @*/ 2824 PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[]) 2825 { 2826 PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i; 2827 DM coordDM = NULL; 2828 Vec coords; 2829 PetscFE fe = NULL; 2830 PetscErrorCode ierr; 2831 2832 PetscFunctionBegin; 2833 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2834 ierr = DMGetDimension(dm,&dimR);CHKERRQ(ierr); 2835 ierr = DMGetCoordinateDim(dm,&dimC);CHKERRQ(ierr); 2836 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 2837 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 2838 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2839 ierr = DMGetCoordinateDM(dm,&coordDM);CHKERRQ(ierr); 2840 if (coordDM) { 2841 PetscInt coordFields; 2842 2843 ierr = DMGetNumFields(coordDM,&coordFields);CHKERRQ(ierr); 2844 if (coordFields) { 2845 PetscClassId id; 2846 PetscObject disc; 2847 2848 ierr = DMGetField(coordDM,0,&disc);CHKERRQ(ierr); 2849 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 2850 if (id == PETSCFE_CLASSID) { 2851 fe = (PetscFE) disc; 2852 } 2853 } 2854 } 2855 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 2856 ierr = DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 2857 cEnd = cEndInterior > 0 ? cEndInterior : cEnd; 2858 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); 2859 if (!fe) { /* implicit discretization: affine or multilinear */ 2860 PetscInt coneSize; 2861 PetscBool isSimplex, isTensor; 2862 2863 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 2864 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 2865 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 2866 if (isSimplex) { 2867 PetscReal detJ, *v0, *J; 2868 2869 ierr = DMGetWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2870 J = &v0[dimC]; 2871 ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ);CHKERRQ(ierr); 2872 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 2873 CoordinatesRefToReal(dimC, dimR, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]); 2874 } 2875 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2876 } else if (isTensor) { 2877 ierr = DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);CHKERRQ(ierr); 2878 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize); 2879 } else { 2880 ierr = DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);CHKERRQ(ierr); 2881 } 2882 PetscFunctionReturn(0); 2883 } 2884