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