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