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