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 const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.); 872 PetscReal M[9], detM; 873 M[0] = x1; M[1] = x2; M[2] = x3; 874 M[3] = y1; M[4] = y2; M[5] = y3; 875 M[6] = z1; M[7] = z2; M[8] = z3; 876 DMPlex_Det3D_Internal(&detM, M); 877 *vol = -onesixth*detM; 878 (void)PetscLogFlops(10.0); 879 } 880 881 PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[]) 882 { 883 const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.); 884 DMPlex_Det3D_Internal(vol, coords); 885 *vol *= -onesixth; 886 } 887 888 static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 889 { 890 PetscSection coordSection; 891 Vec coordinates; 892 const PetscScalar *coords; 893 PetscInt dim, d, off; 894 PetscErrorCode ierr; 895 896 PetscFunctionBegin; 897 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 898 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 899 ierr = PetscSectionGetDof(coordSection,e,&dim);CHKERRQ(ierr); 900 if (!dim) PetscFunctionReturn(0); 901 ierr = PetscSectionGetOffset(coordSection,e,&off);CHKERRQ(ierr); 902 ierr = VecGetArrayRead(coordinates,&coords);CHKERRQ(ierr); 903 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);} 904 ierr = VecRestoreArrayRead(coordinates,&coords);CHKERRQ(ierr); 905 *detJ = 1.; 906 if (J) { 907 for (d = 0; d < dim * dim; d++) J[d] = 0.; 908 for (d = 0; d < dim; d++) J[d * dim + d] = 1.; 909 if (invJ) { 910 for (d = 0; d < dim * dim; d++) invJ[d] = 0.; 911 for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.; 912 } 913 } 914 PetscFunctionReturn(0); 915 } 916 917 static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 918 { 919 PetscSection coordSection; 920 Vec coordinates; 921 PetscScalar *coords = NULL; 922 PetscInt numCoords, d, pStart, pEnd, numSelfCoords = 0; 923 PetscErrorCode ierr; 924 925 PetscFunctionBegin; 926 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 927 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 928 ierr = PetscSectionGetChart(coordSection,&pStart,&pEnd);CHKERRQ(ierr); 929 if (e >= pStart && e < pEnd) {ierr = PetscSectionGetDof(coordSection,e,&numSelfCoords);CHKERRQ(ierr);} 930 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 931 numCoords = numSelfCoords ? numSelfCoords : numCoords; 932 if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL"); 933 *detJ = 0.0; 934 if (numCoords == 6) { 935 const PetscInt dim = 3; 936 PetscReal R[9], J0; 937 938 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 939 ierr = DMPlexComputeProjection3Dto1D(coords, R);CHKERRQ(ierr); 940 if (J) { 941 J0 = 0.5*PetscRealPart(coords[1]); 942 J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2]; 943 J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5]; 944 J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8]; 945 DMPlex_Det3D_Internal(detJ, J); 946 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 947 } 948 } else if (numCoords == 4) { 949 const PetscInt dim = 2; 950 PetscReal R[4], J0; 951 952 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 953 ierr = DMPlexComputeProjection2Dto1D(coords, R);CHKERRQ(ierr); 954 if (J) { 955 J0 = 0.5*PetscRealPart(coords[1]); 956 J[0] = R[0]*J0; J[1] = R[1]; 957 J[2] = R[2]*J0; J[3] = R[3]; 958 DMPlex_Det2D_Internal(detJ, J); 959 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 960 } 961 } else if (numCoords == 2) { 962 const PetscInt dim = 1; 963 964 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 965 if (J) { 966 J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0])); 967 *detJ = J[0]; 968 ierr = PetscLogFlops(2.0);CHKERRQ(ierr); 969 if (invJ) {invJ[0] = 1.0/J[0]; ierr = PetscLogFlops(1.0);CHKERRQ(ierr);} 970 } 971 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords); 972 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 973 PetscFunctionReturn(0); 974 } 975 976 static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 977 { 978 PetscSection coordSection; 979 Vec coordinates; 980 PetscScalar *coords = NULL; 981 PetscInt numCoords, d, f, g; 982 PetscErrorCode ierr; 983 984 PetscFunctionBegin; 985 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 986 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 987 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 988 *detJ = 0.0; 989 if (numCoords == 9) { 990 const PetscInt dim = 3; 991 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 992 993 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 994 ierr = DMPlexComputeProjection3Dto2D(numCoords, coords, R);CHKERRQ(ierr); 995 if (J) { 996 const PetscInt pdim = 2; 997 998 for (d = 0; d < pdim; d++) { 999 for (f = 0; f < pdim; f++) { 1000 J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 1001 } 1002 } 1003 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1004 DMPlex_Det3D_Internal(detJ, J0); 1005 for (d = 0; d < dim; d++) { 1006 for (f = 0; f < dim; f++) { 1007 J[d*dim+f] = 0.0; 1008 for (g = 0; g < dim; g++) { 1009 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 1010 } 1011 } 1012 } 1013 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1014 } 1015 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1016 } else if (numCoords == 6) { 1017 const PetscInt dim = 2; 1018 1019 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1020 if (J) { 1021 for (d = 0; d < dim; d++) { 1022 for (f = 0; f < dim; f++) { 1023 J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 1024 } 1025 } 1026 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1027 DMPlex_Det2D_Internal(detJ, J); 1028 } 1029 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 1030 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords); 1031 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1032 PetscFunctionReturn(0); 1033 } 1034 1035 static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1036 { 1037 PetscSection coordSection; 1038 Vec coordinates; 1039 PetscScalar *coords = NULL; 1040 PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd; 1041 PetscErrorCode ierr; 1042 1043 PetscFunctionBegin; 1044 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1045 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1046 ierr = PetscSectionGetChart(coordSection,&pStart,&pEnd);CHKERRQ(ierr); 1047 if (e >= pStart && e < pEnd) {ierr = PetscSectionGetDof(coordSection,e,&numSelfCoords);CHKERRQ(ierr);} 1048 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1049 numCoords = numSelfCoords ? numSelfCoords : numCoords; 1050 if (!Nq) { 1051 *detJ = 0.0; 1052 if (numCoords == 12) { 1053 const PetscInt dim = 3; 1054 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 1055 1056 if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);} 1057 ierr = DMPlexComputeProjection3Dto2D(numCoords, coords, R);CHKERRQ(ierr); 1058 if (J) { 1059 const PetscInt pdim = 2; 1060 1061 for (d = 0; d < pdim; d++) { 1062 J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 1063 J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 1064 } 1065 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1066 DMPlex_Det3D_Internal(detJ, J0); 1067 for (d = 0; d < dim; d++) { 1068 for (f = 0; f < dim; f++) { 1069 J[d*dim+f] = 0.0; 1070 for (g = 0; g < dim; g++) { 1071 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 1072 } 1073 } 1074 } 1075 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1076 } 1077 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1078 } else if (numCoords == 8) { 1079 const PetscInt dim = 2; 1080 1081 if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);} 1082 if (J) { 1083 for (d = 0; d < dim; d++) { 1084 J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1085 J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1086 } 1087 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1088 DMPlex_Det2D_Internal(detJ, J); 1089 } 1090 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 1091 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords); 1092 } else { 1093 const PetscInt Nv = 4; 1094 const PetscInt dimR = 2; 1095 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 1096 PetscReal zOrder[12]; 1097 PetscReal zCoeff[12]; 1098 PetscInt i, j, k, l, dim; 1099 1100 if (numCoords == 12) { 1101 dim = 3; 1102 } else if (numCoords == 8) { 1103 dim = 2; 1104 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords); 1105 for (i = 0; i < Nv; i++) { 1106 PetscInt zi = zToPlex[i]; 1107 1108 for (j = 0; j < dim; j++) { 1109 zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]); 1110 } 1111 } 1112 for (j = 0; j < dim; j++) { 1113 zCoeff[dim * 0 + j] = 0.25 * ( zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1114 zCoeff[dim * 1 + j] = 0.25 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1115 zCoeff[dim * 2 + j] = 0.25 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1116 zCoeff[dim * 3 + j] = 0.25 * ( zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]); 1117 } 1118 for (i = 0; i < Nq; i++) { 1119 PetscReal xi = points[dimR * i], eta = points[dimR * i + 1]; 1120 1121 if (v) { 1122 PetscReal extPoint[4]; 1123 1124 extPoint[0] = 1.; 1125 extPoint[1] = xi; 1126 extPoint[2] = eta; 1127 extPoint[3] = xi * eta; 1128 for (j = 0; j < dim; j++) { 1129 PetscReal val = 0.; 1130 1131 for (k = 0; k < Nv; k++) { 1132 val += extPoint[k] * zCoeff[dim * k + j]; 1133 } 1134 v[i * dim + j] = val; 1135 } 1136 } 1137 if (J) { 1138 PetscReal extJ[8]; 1139 1140 extJ[0] = 0.; 1141 extJ[1] = 0.; 1142 extJ[2] = 1.; 1143 extJ[3] = 0.; 1144 extJ[4] = 0.; 1145 extJ[5] = 1.; 1146 extJ[6] = eta; 1147 extJ[7] = xi; 1148 for (j = 0; j < dim; j++) { 1149 for (k = 0; k < dimR; k++) { 1150 PetscReal val = 0.; 1151 1152 for (l = 0; l < Nv; l++) { 1153 val += zCoeff[dim * l + j] * extJ[dimR * l + k]; 1154 } 1155 J[i * dim * dim + dim * j + k] = val; 1156 } 1157 } 1158 if (dim == 3) { /* put the cross product in the third component of the Jacobian */ 1159 PetscReal x, y, z; 1160 PetscReal *iJ = &J[i * dim * dim]; 1161 PetscReal norm; 1162 1163 x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0]; 1164 y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1]; 1165 z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0]; 1166 norm = PetscSqrtReal(x * x + y * y + z * z); 1167 iJ[2] = x / norm; 1168 iJ[5] = y / norm; 1169 iJ[8] = z / norm; 1170 DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]); 1171 if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);} 1172 } else { 1173 DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]); 1174 if (invJ) {DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);} 1175 } 1176 } 1177 } 1178 } 1179 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1180 PetscFunctionReturn(0); 1181 } 1182 1183 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1184 { 1185 PetscSection coordSection; 1186 Vec coordinates; 1187 PetscScalar *coords = NULL; 1188 const PetscInt dim = 3; 1189 PetscInt d; 1190 PetscErrorCode ierr; 1191 1192 PetscFunctionBegin; 1193 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1194 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1195 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1196 *detJ = 0.0; 1197 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1198 if (J) { 1199 for (d = 0; d < dim; d++) { 1200 /* I orient with outward face normals */ 1201 J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d])); 1202 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1203 J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1204 } 1205 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1206 DMPlex_Det3D_Internal(detJ, J); 1207 } 1208 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1209 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1210 PetscFunctionReturn(0); 1211 } 1212 1213 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1214 { 1215 PetscSection coordSection; 1216 Vec coordinates; 1217 PetscScalar *coords = NULL; 1218 const PetscInt dim = 3; 1219 PetscInt d; 1220 PetscErrorCode ierr; 1221 1222 PetscFunctionBegin; 1223 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1224 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1225 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1226 if (!Nq) { 1227 *detJ = 0.0; 1228 if (v) {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);} 1229 if (J) { 1230 for (d = 0; d < dim; d++) { 1231 J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1232 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1233 J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d])); 1234 } 1235 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1236 DMPlex_Det3D_Internal(detJ, J); 1237 } 1238 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1239 } else { 1240 const PetscInt Nv = 8; 1241 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 1242 const PetscInt dim = 3; 1243 const PetscInt dimR = 3; 1244 PetscReal zOrder[24]; 1245 PetscReal zCoeff[24]; 1246 PetscInt i, j, k, l; 1247 1248 for (i = 0; i < Nv; i++) { 1249 PetscInt zi = zToPlex[i]; 1250 1251 for (j = 0; j < dim; j++) { 1252 zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]); 1253 } 1254 } 1255 for (j = 0; j < dim; j++) { 1256 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]); 1257 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]); 1258 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]); 1259 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]); 1260 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]); 1261 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]); 1262 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]); 1263 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]); 1264 } 1265 for (i = 0; i < Nq; i++) { 1266 PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2]; 1267 1268 if (v) { 1269 PetscReal extPoint[8]; 1270 1271 extPoint[0] = 1.; 1272 extPoint[1] = xi; 1273 extPoint[2] = eta; 1274 extPoint[3] = xi * eta; 1275 extPoint[4] = theta; 1276 extPoint[5] = theta * xi; 1277 extPoint[6] = theta * eta; 1278 extPoint[7] = theta * eta * xi; 1279 for (j = 0; j < dim; j++) { 1280 PetscReal val = 0.; 1281 1282 for (k = 0; k < Nv; k++) { 1283 val += extPoint[k] * zCoeff[dim * k + j]; 1284 } 1285 v[i * dim + j] = val; 1286 } 1287 } 1288 if (J) { 1289 PetscReal extJ[24]; 1290 1291 extJ[0] = 0. ; extJ[1] = 0. ; extJ[2] = 0. ; 1292 extJ[3] = 1. ; extJ[4] = 0. ; extJ[5] = 0. ; 1293 extJ[6] = 0. ; extJ[7] = 1. ; extJ[8] = 0. ; 1294 extJ[9] = eta ; extJ[10] = xi ; extJ[11] = 0. ; 1295 extJ[12] = 0. ; extJ[13] = 0. ; extJ[14] = 1. ; 1296 extJ[15] = theta ; extJ[16] = 0. ; extJ[17] = xi ; 1297 extJ[18] = 0. ; extJ[19] = theta ; extJ[20] = eta ; 1298 extJ[21] = theta * eta; extJ[22] = theta * xi; extJ[23] = eta * xi; 1299 1300 for (j = 0; j < dim; j++) { 1301 for (k = 0; k < dimR; k++) { 1302 PetscReal val = 0.; 1303 1304 for (l = 0; l < Nv; l++) { 1305 val += zCoeff[dim * l + j] * extJ[dimR * l + k]; 1306 } 1307 J[i * dim * dim + dim * j + k] = val; 1308 } 1309 } 1310 DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]); 1311 if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);} 1312 } 1313 } 1314 } 1315 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1316 PetscFunctionReturn(0); 1317 } 1318 1319 static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1320 { 1321 PetscInt depth, dim, coordDim, coneSize, i; 1322 PetscInt Nq = 0; 1323 const PetscReal *points = NULL; 1324 DMLabel depthLabel; 1325 PetscReal v0[3] = {-1.}, J0[9] = {-1.}, detJ0 = -1.; 1326 PetscBool isAffine = PETSC_TRUE; 1327 PetscErrorCode ierr; 1328 1329 PetscFunctionBegin; 1330 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1331 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 1332 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 1333 ierr = DMLabelGetValue(depthLabel, cell, &dim);CHKERRQ(ierr); 1334 if (depth == 1 && dim == 1) { 1335 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1336 } 1337 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1338 if (coordDim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %D > 3", coordDim); 1339 if (quad) {ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL);CHKERRQ(ierr);} 1340 switch (dim) { 1341 case 0: 1342 ierr = DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);CHKERRQ(ierr); 1343 isAffine = PETSC_FALSE; 1344 break; 1345 case 1: 1346 if (Nq) { 1347 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);CHKERRQ(ierr); 1348 } else { 1349 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ);CHKERRQ(ierr); 1350 } 1351 break; 1352 case 2: 1353 switch (coneSize) { 1354 case 3: 1355 if (Nq) { 1356 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);CHKERRQ(ierr); 1357 } else { 1358 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ);CHKERRQ(ierr); 1359 } 1360 break; 1361 case 4: 1362 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);CHKERRQ(ierr); 1363 isAffine = PETSC_FALSE; 1364 break; 1365 default: 1366 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 1367 } 1368 break; 1369 case 3: 1370 switch (coneSize) { 1371 case 4: 1372 if (Nq) { 1373 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);CHKERRQ(ierr); 1374 } else { 1375 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ);CHKERRQ(ierr); 1376 } 1377 break; 1378 case 6: /* Faces */ 1379 case 8: /* Vertices */ 1380 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);CHKERRQ(ierr); 1381 isAffine = PETSC_FALSE; 1382 break; 1383 default: 1384 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 1385 } 1386 break; 1387 default: 1388 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1389 } 1390 if (isAffine && Nq) { 1391 if (v) { 1392 for (i = 0; i < Nq; i++) { 1393 CoordinatesRefToReal(coordDim,dim,v0, J0, &points[dim * i], &v[coordDim * i]); 1394 } 1395 } 1396 if (detJ) { 1397 for (i = 0; i < Nq; i++) { 1398 detJ[i] = detJ0; 1399 } 1400 } 1401 if (J) { 1402 PetscInt k; 1403 1404 for (i = 0, k = 0; i < Nq; i++) { 1405 PetscInt j; 1406 1407 for (j = 0; j < coordDim * coordDim; j++, k++) { 1408 J[k] = J0[j]; 1409 } 1410 } 1411 } 1412 if (invJ) { 1413 PetscInt k; 1414 switch (coordDim) { 1415 case 0: 1416 break; 1417 case 1: 1418 invJ[0] = 1./J0[0]; 1419 break; 1420 case 2: 1421 DMPlex_Invert2D_Internal(invJ, J0, detJ0); 1422 break; 1423 case 3: 1424 DMPlex_Invert3D_Internal(invJ, J0, detJ0); 1425 break; 1426 } 1427 for (i = 1, k = coordDim * coordDim; i < Nq; i++) { 1428 PetscInt j; 1429 1430 for (j = 0; j < coordDim * coordDim; j++, k++) { 1431 invJ[k] = invJ[j]; 1432 } 1433 } 1434 } 1435 } 1436 PetscFunctionReturn(0); 1437 } 1438 1439 /*@C 1440 DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell 1441 1442 Collective on DM 1443 1444 Input Arguments: 1445 + dm - the DM 1446 - cell - the cell 1447 1448 Output Arguments: 1449 + v0 - the translation part of this affine transform 1450 . J - the Jacobian of the transform from the reference element 1451 . invJ - the inverse of the Jacobian 1452 - detJ - the Jacobian determinant 1453 1454 Level: advanced 1455 1456 Fortran Notes: 1457 Since it returns arrays, this routine is only available in Fortran 90, and you must 1458 include petsc.h90 in your code. 1459 1460 .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinateVec() 1461 @*/ 1462 PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1463 { 1464 PetscErrorCode ierr; 1465 1466 PetscFunctionBegin; 1467 ierr = DMPlexComputeCellGeometryFEM_Implicit(dm,cell,NULL,v0,J,invJ,detJ);CHKERRQ(ierr); 1468 PetscFunctionReturn(0); 1469 } 1470 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, NULL, &Nq, &quadPoints, NULL);CHKERRQ(ierr); 1494 Nq = 1; 1495 } else { 1496 ierr = PetscQuadratureGetData(quad, &qdim, NULL, &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 /*@ 1848 DMPlexComputeGeometryFEM - Precompute cell geometry for the entire mesh 1849 1850 Collective on dm 1851 1852 Input Parameter: 1853 . dm - The DMPlex 1854 1855 Output Parameter: 1856 . cellgeom - A vector with the cell geometry data for each cell 1857 1858 Level: beginner 1859 1860 .keywords: DMPlexComputeCellGeometryFEM() 1861 @*/ 1862 PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom) 1863 { 1864 DM dmCell; 1865 Vec coordinates; 1866 PetscSection coordSection, sectionCell; 1867 PetscScalar *cgeom; 1868 PetscInt cStart, cEnd, cMax, c; 1869 PetscErrorCode ierr; 1870 1871 PetscFunctionBegin; 1872 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1873 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1874 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1875 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1876 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1877 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1878 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1879 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1880 cEnd = cMax < 0 ? cEnd : cMax; 1881 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1882 /* TODO This needs to be multiplied by Nq for non-affine */ 1883 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFECellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1884 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1885 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1886 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1887 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1888 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1889 for (c = cStart; c < cEnd; ++c) { 1890 PetscFECellGeom *cg; 1891 1892 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1893 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1894 ierr = DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v0, cg->J, cg->invJ, &cg->detJ);CHKERRQ(ierr); 1895 if (cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c); 1896 } 1897 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1898 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1899 PetscFunctionReturn(0); 1900 } 1901 1902 /*@ 1903 DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method 1904 1905 Input Parameter: 1906 . dm - The DM 1907 1908 Output Parameters: 1909 + cellgeom - A Vec of PetscFVCellGeom data 1910 . facegeom - A Vec of PetscFVFaceGeom data 1911 1912 Level: developer 1913 1914 .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM() 1915 @*/ 1916 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 1917 { 1918 DM dmFace, dmCell; 1919 DMLabel ghostLabel; 1920 PetscSection sectionFace, sectionCell; 1921 PetscSection coordSection; 1922 Vec coordinates; 1923 PetscScalar *fgeom, *cgeom; 1924 PetscReal minradius, gminradius; 1925 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 1926 PetscErrorCode ierr; 1927 1928 PetscFunctionBegin; 1929 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1930 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1931 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1932 /* Make cell centroids and volumes */ 1933 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1934 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1935 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1936 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1937 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1938 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1939 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1940 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1941 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1942 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1943 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1944 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1945 if (cEndInterior < 0) { 1946 cEndInterior = cEnd; 1947 } 1948 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1949 for (c = cStart; c < cEndInterior; ++c) { 1950 PetscFVCellGeom *cg; 1951 1952 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1953 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1954 ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr); 1955 } 1956 /* Compute face normals and minimum cell radius */ 1957 ierr = DMClone(dm, &dmFace);CHKERRQ(ierr); 1958 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);CHKERRQ(ierr); 1959 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1960 ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr); 1961 for (f = fStart; f < fEnd; ++f) {ierr = PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1962 ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr); 1963 ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr); 1964 ierr = PetscSectionDestroy(§ionFace);CHKERRQ(ierr); 1965 ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr); 1966 ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr); 1967 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1968 minradius = PETSC_MAX_REAL; 1969 for (f = fStart; f < fEnd; ++f) { 1970 PetscFVFaceGeom *fg; 1971 PetscReal area; 1972 PetscInt ghost = -1, d, numChildren; 1973 1974 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1975 ierr = DMPlexGetTreeChildren(dm,f,&numChildren,NULL);CHKERRQ(ierr); 1976 if (ghost >= 0 || numChildren) continue; 1977 ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr); 1978 ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);CHKERRQ(ierr); 1979 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 1980 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 1981 { 1982 PetscFVCellGeom *cL, *cR; 1983 PetscInt ncells; 1984 const PetscInt *cells; 1985 PetscReal *lcentroid, *rcentroid; 1986 PetscReal l[3], r[3], v[3]; 1987 1988 ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr); 1989 ierr = DMPlexGetSupportSize(dm, f, &ncells);CHKERRQ(ierr); 1990 ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr); 1991 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 1992 if (ncells > 1) { 1993 ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr); 1994 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 1995 } 1996 else { 1997 rcentroid = fg->centroid; 1998 } 1999 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);CHKERRQ(ierr); 2000 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);CHKERRQ(ierr); 2001 DMPlex_WaxpyD_Internal(dim, -1, l, r, v); 2002 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 2003 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 2004 } 2005 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 2006 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]); 2007 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]); 2008 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f); 2009 } 2010 if (cells[0] < cEndInterior) { 2011 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 2012 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2013 } 2014 if (ncells > 1 && cells[1] < cEndInterior) { 2015 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 2016 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 2017 } 2018 } 2019 } 2020 ierr = MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 2021 ierr = DMPlexSetMinRadius(dm, gminradius);CHKERRQ(ierr); 2022 /* Compute centroids of ghost cells */ 2023 for (c = cEndInterior; c < cEnd; ++c) { 2024 PetscFVFaceGeom *fg; 2025 const PetscInt *cone, *support; 2026 PetscInt coneSize, supportSize, s; 2027 2028 ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr); 2029 if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize); 2030 ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr); 2031 ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr); 2032 if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize); 2033 ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr); 2034 ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr); 2035 for (s = 0; s < 2; ++s) { 2036 /* Reflect ghost centroid across plane of face */ 2037 if (support[s] == c) { 2038 PetscFVCellGeom *ci; 2039 PetscFVCellGeom *cg; 2040 PetscReal c2f[3], a; 2041 2042 ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr); 2043 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 2044 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 2045 ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr); 2046 DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid); 2047 cg->volume = ci->volume; 2048 } 2049 } 2050 } 2051 ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr); 2052 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 2053 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 2054 ierr = DMDestroy(&dmFace);CHKERRQ(ierr); 2055 PetscFunctionReturn(0); 2056 } 2057 2058 /*@C 2059 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 2060 2061 Not collective 2062 2063 Input Argument: 2064 . dm - the DM 2065 2066 Output Argument: 2067 . minradius - the minium cell radius 2068 2069 Level: developer 2070 2071 .seealso: DMGetCoordinates() 2072 @*/ 2073 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 2074 { 2075 PetscFunctionBegin; 2076 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2077 PetscValidPointer(minradius,2); 2078 *minradius = ((DM_Plex*) dm->data)->minradius; 2079 PetscFunctionReturn(0); 2080 } 2081 2082 /*@C 2083 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 2084 2085 Logically collective 2086 2087 Input Arguments: 2088 + dm - the DM 2089 - minradius - the minium cell radius 2090 2091 Level: developer 2092 2093 .seealso: DMSetCoordinates() 2094 @*/ 2095 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 2096 { 2097 PetscFunctionBegin; 2098 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2099 ((DM_Plex*) dm->data)->minradius = minradius; 2100 PetscFunctionReturn(0); 2101 } 2102 2103 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 2104 { 2105 DMLabel ghostLabel; 2106 PetscScalar *dx, *grad, **gref; 2107 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 2108 PetscErrorCode ierr; 2109 2110 PetscFunctionBegin; 2111 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2112 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2113 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2114 ierr = DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);CHKERRQ(ierr); 2115 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 2116 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 2117 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 2118 for (c = cStart; c < cEndInterior; c++) { 2119 const PetscInt *faces; 2120 PetscInt numFaces, usedFaces, f, d; 2121 PetscFVCellGeom *cg; 2122 PetscBool boundary; 2123 PetscInt ghost; 2124 2125 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2126 ierr = DMPlexGetConeSize(dm, c, &numFaces);CHKERRQ(ierr); 2127 ierr = DMPlexGetCone(dm, c, &faces);CHKERRQ(ierr); 2128 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 2129 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 2130 PetscFVCellGeom *cg1; 2131 PetscFVFaceGeom *fg; 2132 const PetscInt *fcells; 2133 PetscInt ncell, side; 2134 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 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr); 2139 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 2140 ncell = fcells[!side]; /* the neighbor */ 2141 ierr = DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr); 2142 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 2143 for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d]; 2144 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2145 } 2146 if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 2147 ierr = PetscFVComputeGradient(fvm, usedFaces, dx, grad);CHKERRQ(ierr); 2148 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 2149 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 2150 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 2151 if ((ghost >= 0) || boundary) continue; 2152 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d]; 2153 ++usedFaces; 2154 } 2155 } 2156 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 2157 PetscFunctionReturn(0); 2158 } 2159 2160 static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 2161 { 2162 DMLabel ghostLabel; 2163 PetscScalar *dx, *grad, **gref; 2164 PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0; 2165 PetscSection neighSec; 2166 PetscInt (*neighbors)[2]; 2167 PetscInt *counter; 2168 PetscErrorCode ierr; 2169 2170 PetscFunctionBegin; 2171 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2172 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2173 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2174 if (cEndInterior < 0) { 2175 cEndInterior = cEnd; 2176 } 2177 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);CHKERRQ(ierr); 2178 ierr = PetscSectionSetChart(neighSec,cStart,cEndInterior);CHKERRQ(ierr); 2179 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 2180 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 2181 for (f = fStart; f < fEnd; f++) { 2182 const PetscInt *fcells; 2183 PetscBool boundary; 2184 PetscInt ghost = -1; 2185 PetscInt numChildren, numCells, c; 2186 2187 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 2188 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 2189 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 2190 if ((ghost >= 0) || boundary || numChildren) continue; 2191 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 2192 if (numCells == 2) { 2193 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 2194 for (c = 0; c < 2; c++) { 2195 PetscInt cell = fcells[c]; 2196 2197 if (cell >= cStart && cell < cEndInterior) { 2198 ierr = PetscSectionAddDof(neighSec,cell,1);CHKERRQ(ierr); 2199 } 2200 } 2201 } 2202 } 2203 ierr = PetscSectionSetUp(neighSec);CHKERRQ(ierr); 2204 ierr = PetscSectionGetMaxDof(neighSec,&maxNumFaces);CHKERRQ(ierr); 2205 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 2206 nStart = 0; 2207 ierr = PetscSectionGetStorageSize(neighSec,&nEnd);CHKERRQ(ierr); 2208 ierr = PetscMalloc1((nEnd-nStart),&neighbors);CHKERRQ(ierr); 2209 ierr = PetscCalloc1((cEndInterior-cStart),&counter);CHKERRQ(ierr); 2210 for (f = fStart; f < fEnd; f++) { 2211 const PetscInt *fcells; 2212 PetscBool boundary; 2213 PetscInt ghost = -1; 2214 PetscInt numChildren, numCells, c; 2215 2216 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 2217 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 2218 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 2219 if ((ghost >= 0) || boundary || numChildren) continue; 2220 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 2221 if (numCells == 2) { 2222 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 2223 for (c = 0; c < 2; c++) { 2224 PetscInt cell = fcells[c], off; 2225 2226 if (cell >= cStart && cell < cEndInterior) { 2227 ierr = PetscSectionGetOffset(neighSec,cell,&off);CHKERRQ(ierr); 2228 off += counter[cell - cStart]++; 2229 neighbors[off][0] = f; 2230 neighbors[off][1] = fcells[1 - c]; 2231 } 2232 } 2233 } 2234 } 2235 ierr = PetscFree(counter);CHKERRQ(ierr); 2236 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 2237 for (c = cStart; c < cEndInterior; c++) { 2238 PetscInt numFaces, f, d, off, ghost = -1; 2239 PetscFVCellGeom *cg; 2240 2241 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2242 ierr = PetscSectionGetDof(neighSec, c, &numFaces);CHKERRQ(ierr); 2243 ierr = PetscSectionGetOffset(neighSec, c, &off);CHKERRQ(ierr); 2244 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, c, &ghost);CHKERRQ(ierr);} 2245 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); 2246 for (f = 0; f < numFaces; ++f) { 2247 PetscFVCellGeom *cg1; 2248 PetscFVFaceGeom *fg; 2249 const PetscInt *fcells; 2250 PetscInt ncell, side, nface; 2251 2252 nface = neighbors[off + f][0]; 2253 ncell = neighbors[off + f][1]; 2254 ierr = DMPlexGetSupport(dm,nface,&fcells);CHKERRQ(ierr); 2255 side = (c != fcells[0]); 2256 ierr = DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);CHKERRQ(ierr); 2257 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 2258 for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d]; 2259 gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2260 } 2261 ierr = PetscFVComputeGradient(fvm, numFaces, dx, grad);CHKERRQ(ierr); 2262 for (f = 0; f < numFaces; ++f) { 2263 for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d]; 2264 } 2265 } 2266 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 2267 ierr = PetscSectionDestroy(&neighSec);CHKERRQ(ierr); 2268 ierr = PetscFree(neighbors);CHKERRQ(ierr); 2269 PetscFunctionReturn(0); 2270 } 2271 2272 /*@ 2273 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 2274 2275 Collective on DM 2276 2277 Input Arguments: 2278 + dm - The DM 2279 . fvm - The PetscFV 2280 . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM() 2281 - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM() 2282 2283 Output Parameters: 2284 + faceGeometry - The geometric factors for gradient calculation are inserted 2285 - dmGrad - The DM describing the layout of gradient data 2286 2287 Level: developer 2288 2289 .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM() 2290 @*/ 2291 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 2292 { 2293 DM dmFace, dmCell; 2294 PetscScalar *fgeom, *cgeom; 2295 PetscSection sectionGrad, parentSection; 2296 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 2297 PetscErrorCode ierr; 2298 2299 PetscFunctionBegin; 2300 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2301 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 2302 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2303 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2304 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 2305 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 2306 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 2307 ierr = VecGetArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2308 ierr = VecGetArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2309 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 2310 if (!parentSection) { 2311 ierr = BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2312 } else { 2313 ierr = BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2314 } 2315 ierr = VecRestoreArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2316 ierr = VecRestoreArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2317 /* Create storage for gradients */ 2318 ierr = DMClone(dm, dmGrad);CHKERRQ(ierr); 2319 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 2320 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 2321 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 2322 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 2323 ierr = DMSetDefaultSection(*dmGrad, sectionGrad);CHKERRQ(ierr); 2324 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 2325 PetscFunctionReturn(0); 2326 } 2327 2328 /*@ 2329 DMPlexGetDataFVM - Retrieve precomputed cell geometry 2330 2331 Collective on DM 2332 2333 Input Arguments: 2334 + dm - The DM 2335 - fvm - The PetscFV 2336 2337 Output Parameters: 2338 + cellGeometry - The cell geometry 2339 . faceGeometry - The face geometry 2340 - dmGrad - The gradient matrices 2341 2342 Level: developer 2343 2344 .seealso: DMPlexComputeGeometryFVM() 2345 @*/ 2346 PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM) 2347 { 2348 PetscObject cellgeomobj, facegeomobj; 2349 PetscErrorCode ierr; 2350 2351 PetscFunctionBegin; 2352 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2353 if (!cellgeomobj) { 2354 Vec cellgeomInt, facegeomInt; 2355 2356 ierr = DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);CHKERRQ(ierr); 2357 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);CHKERRQ(ierr); 2358 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);CHKERRQ(ierr); 2359 ierr = VecDestroy(&cellgeomInt);CHKERRQ(ierr); 2360 ierr = VecDestroy(&facegeomInt);CHKERRQ(ierr); 2361 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2362 } 2363 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);CHKERRQ(ierr); 2364 if (cellgeom) *cellgeom = (Vec) cellgeomobj; 2365 if (facegeom) *facegeom = (Vec) facegeomobj; 2366 if (gradDM) { 2367 PetscObject gradobj; 2368 PetscBool computeGradients; 2369 2370 ierr = PetscFVGetComputeGradients(fv,&computeGradients);CHKERRQ(ierr); 2371 if (!computeGradients) { 2372 *gradDM = NULL; 2373 PetscFunctionReturn(0); 2374 } 2375 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2376 if (!gradobj) { 2377 DM dmGradInt; 2378 2379 ierr = DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);CHKERRQ(ierr); 2380 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);CHKERRQ(ierr); 2381 ierr = DMDestroy(&dmGradInt);CHKERRQ(ierr); 2382 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2383 } 2384 *gradDM = (DM) gradobj; 2385 } 2386 PetscFunctionReturn(0); 2387 } 2388 2389 static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess) 2390 { 2391 PetscInt l, m; 2392 2393 PetscFunctionBeginHot; 2394 if (dimC == dimR && dimR <= 3) { 2395 /* invert Jacobian, multiply */ 2396 PetscScalar det, idet; 2397 2398 switch (dimR) { 2399 case 1: 2400 invJ[0] = 1./ J[0]; 2401 break; 2402 case 2: 2403 det = J[0] * J[3] - J[1] * J[2]; 2404 idet = 1./det; 2405 invJ[0] = J[3] * idet; 2406 invJ[1] = -J[1] * idet; 2407 invJ[2] = -J[2] * idet; 2408 invJ[3] = J[0] * idet; 2409 break; 2410 case 3: 2411 { 2412 invJ[0] = J[4] * J[8] - J[5] * J[7]; 2413 invJ[1] = J[2] * J[7] - J[1] * J[8]; 2414 invJ[2] = J[1] * J[5] - J[2] * J[4]; 2415 det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6]; 2416 idet = 1./det; 2417 invJ[0] *= idet; 2418 invJ[1] *= idet; 2419 invJ[2] *= idet; 2420 invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]); 2421 invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]); 2422 invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]); 2423 invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]); 2424 invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]); 2425 invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]); 2426 } 2427 break; 2428 } 2429 for (l = 0; l < dimR; l++) { 2430 for (m = 0; m < dimC; m++) { 2431 guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m]; 2432 } 2433 } 2434 } else { 2435 #if defined(PETSC_USE_COMPLEX) 2436 char transpose = 'C'; 2437 #else 2438 char transpose = 'T'; 2439 #endif 2440 PetscBLASInt m = dimR; 2441 PetscBLASInt n = dimC; 2442 PetscBLASInt one = 1; 2443 PetscBLASInt worksize = dimR * dimC, info; 2444 2445 for (l = 0; l < dimC; l++) {invJ[l] = resNeg[l];} 2446 2447 PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&m,&n,&one,J,&m,invJ,&n,work,&worksize, &info)); 2448 if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS"); 2449 2450 for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);} 2451 } 2452 PetscFunctionReturn(0); 2453 } 2454 2455 static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2456 { 2457 PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR); 2458 PetscScalar *coordsScalar = NULL; 2459 PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg; 2460 PetscScalar *J, *invJ, *work; 2461 PetscErrorCode ierr; 2462 2463 PetscFunctionBegin; 2464 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2465 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2466 if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);CHKERRQ(ierr); 2467 ierr = DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, PETSC_REAL, &cellData);CHKERRQ(ierr); 2468 ierr = DMGetWorkArray(dm, 3 * dimR * dimC, PETSC_SCALAR, &J);CHKERRQ(ierr); 2469 cellCoords = &cellData[0]; 2470 cellCoeffs = &cellData[coordSize]; 2471 extJ = &cellData[2 * coordSize]; 2472 resNeg = &cellData[2 * coordSize + dimR]; 2473 invJ = &J[dimR * dimC]; 2474 work = &J[2 * dimR * dimC]; 2475 if (dimR == 2) { 2476 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 2477 2478 for (i = 0; i < 4; i++) { 2479 PetscInt plexI = zToPlex[i]; 2480 2481 for (j = 0; j < dimC; j++) { 2482 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2483 } 2484 } 2485 } else if (dimR == 3) { 2486 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2487 2488 for (i = 0; i < 8; i++) { 2489 PetscInt plexI = zToPlex[i]; 2490 2491 for (j = 0; j < dimC; j++) { 2492 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2493 } 2494 } 2495 } else { 2496 for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);} 2497 } 2498 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 2499 for (i = 0; i < dimR; i++) { 2500 PetscReal *swap; 2501 2502 for (j = 0; j < (numV / 2); j++) { 2503 for (k = 0; k < dimC; k++) { 2504 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 2505 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 2506 } 2507 } 2508 2509 if (i < dimR - 1) { 2510 swap = cellCoeffs; 2511 cellCoeffs = cellCoords; 2512 cellCoords = swap; 2513 } 2514 } 2515 ierr = PetscMemzero(refCoords,numPoints * dimR * sizeof (PetscReal));CHKERRQ(ierr); 2516 for (j = 0; j < numPoints; j++) { 2517 for (i = 0; i < maxIts; i++) { 2518 PetscReal *guess = &refCoords[dimR * j]; 2519 2520 /* compute -residual and Jacobian */ 2521 for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];} 2522 for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;} 2523 for (k = 0; k < numV; k++) { 2524 PetscReal extCoord = 1.; 2525 for (l = 0; l < dimR; l++) { 2526 PetscReal coord = guess[l]; 2527 PetscInt dep = (k & (1 << l)) >> l; 2528 2529 extCoord *= dep * coord + !dep; 2530 extJ[l] = dep; 2531 2532 for (m = 0; m < dimR; m++) { 2533 PetscReal coord = guess[m]; 2534 PetscInt dep = ((k & (1 << m)) >> m) && (m != l); 2535 PetscReal mult = dep * coord + !dep; 2536 2537 extJ[l] *= mult; 2538 } 2539 } 2540 for (l = 0; l < dimC; l++) { 2541 PetscReal coeff = cellCoeffs[dimC * k + l]; 2542 2543 resNeg[l] -= coeff * extCoord; 2544 for (m = 0; m < dimR; m++) { 2545 J[dimR * l + m] += coeff * extJ[m]; 2546 } 2547 } 2548 } 2549 #if 0 && defined(PETSC_USE_DEBUG) 2550 { 2551 PetscReal maxAbs = 0.; 2552 2553 for (l = 0; l < dimC; l++) { 2554 maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l])); 2555 } 2556 ierr = PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,maxAbs);CHKERRQ(ierr); 2557 } 2558 #endif 2559 2560 ierr = DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);CHKERRQ(ierr); 2561 } 2562 } 2563 ierr = DMRestoreWorkArray(dm, 3 * dimR * dimC, PETSC_SCALAR, &J);CHKERRQ(ierr); 2564 ierr = DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, PETSC_REAL, &cellData);CHKERRQ(ierr); 2565 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2566 PetscFunctionReturn(0); 2567 } 2568 2569 static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2570 { 2571 PetscInt coordSize, i, j, k, l, numV = (1 << dimR); 2572 PetscScalar *coordsScalar = NULL; 2573 PetscReal *cellData, *cellCoords, *cellCoeffs; 2574 PetscErrorCode ierr; 2575 2576 PetscFunctionBegin; 2577 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2578 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2579 if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);CHKERRQ(ierr); 2580 ierr = DMGetWorkArray(dm, 2 * coordSize, PETSC_REAL, &cellData);CHKERRQ(ierr); 2581 cellCoords = &cellData[0]; 2582 cellCoeffs = &cellData[coordSize]; 2583 if (dimR == 2) { 2584 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 2585 2586 for (i = 0; i < 4; i++) { 2587 PetscInt plexI = zToPlex[i]; 2588 2589 for (j = 0; j < dimC; j++) { 2590 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2591 } 2592 } 2593 } else if (dimR == 3) { 2594 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2595 2596 for (i = 0; i < 8; i++) { 2597 PetscInt plexI = zToPlex[i]; 2598 2599 for (j = 0; j < dimC; j++) { 2600 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2601 } 2602 } 2603 } else { 2604 for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);} 2605 } 2606 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 2607 for (i = 0; i < dimR; i++) { 2608 PetscReal *swap; 2609 2610 for (j = 0; j < (numV / 2); j++) { 2611 for (k = 0; k < dimC; k++) { 2612 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 2613 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 2614 } 2615 } 2616 2617 if (i < dimR - 1) { 2618 swap = cellCoeffs; 2619 cellCoeffs = cellCoords; 2620 cellCoords = swap; 2621 } 2622 } 2623 ierr = PetscMemzero(realCoords,numPoints * dimC * sizeof (PetscReal));CHKERRQ(ierr); 2624 for (j = 0; j < numPoints; j++) { 2625 const PetscReal *guess = &refCoords[dimR * j]; 2626 PetscReal *mapped = &realCoords[dimC * j]; 2627 2628 for (k = 0; k < numV; k++) { 2629 PetscReal extCoord = 1.; 2630 for (l = 0; l < dimR; l++) { 2631 PetscReal coord = guess[l]; 2632 PetscInt dep = (k & (1 << l)) >> l; 2633 2634 extCoord *= dep * coord + !dep; 2635 } 2636 for (l = 0; l < dimC; l++) { 2637 PetscReal coeff = cellCoeffs[dimC * k + l]; 2638 2639 mapped[l] += coeff * extCoord; 2640 } 2641 } 2642 } 2643 ierr = DMRestoreWorkArray(dm, 2 * coordSize, PETSC_REAL, &cellData);CHKERRQ(ierr); 2644 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2645 PetscFunctionReturn(0); 2646 } 2647 2648 /* TODO: TOBY please fix this for Nc > 1 */ 2649 static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR) 2650 { 2651 PetscInt numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize; 2652 PetscScalar *nodes = NULL; 2653 PetscReal *invV, *modes; 2654 PetscReal *B, *D, *resNeg; 2655 PetscScalar *J, *invJ, *work; 2656 PetscErrorCode ierr; 2657 2658 PetscFunctionBegin; 2659 ierr = PetscFEGetDimension(fe, &pdim);CHKERRQ(ierr); 2660 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 2661 if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc); 2662 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2663 /* convert nodes to values in the stable evaluation basis */ 2664 ierr = DMGetWorkArray(dm,pdim,PETSC_REAL,&modes);CHKERRQ(ierr); 2665 invV = fe->invV; 2666 for (i = 0; i < pdim; ++i) { 2667 modes[i] = 0.; 2668 for (j = 0; j < pdim; ++j) { 2669 modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]); 2670 } 2671 } 2672 ierr = DMGetWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,PETSC_REAL,&B);CHKERRQ(ierr); 2673 D = &B[pdim*Nc]; 2674 resNeg = &D[pdim*Nc * dimR]; 2675 ierr = DMGetWorkArray(dm,3 * Nc * dimR,PETSC_SCALAR,&J);CHKERRQ(ierr); 2676 invJ = &J[Nc * dimR]; 2677 work = &invJ[Nc * dimR]; 2678 for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;} 2679 for (j = 0; j < numPoints; j++) { 2680 for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */ 2681 PetscReal *guess = &refCoords[j * dimR]; 2682 ierr = PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);CHKERRQ(ierr); 2683 for (k = 0; k < Nc; k++) {resNeg[k] = realCoords[j * Nc + k];} 2684 for (k = 0; k < Nc * dimR; k++) {J[k] = 0.;} 2685 for (k = 0; k < pdim; k++) { 2686 for (l = 0; l < Nc; l++) { 2687 resNeg[l] -= modes[k] * B[k * Nc + l]; 2688 for (m = 0; m < dimR; m++) { 2689 J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m]; 2690 } 2691 } 2692 } 2693 #if 0 && defined(PETSC_USE_DEBUG) 2694 { 2695 PetscReal maxAbs = 0.; 2696 2697 for (l = 0; l < Nc; l++) { 2698 maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l])); 2699 } 2700 ierr = PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,maxAbs);CHKERRQ(ierr); 2701 } 2702 #endif 2703 ierr = DMPlexCoordinatesToReference_NewtonUpdate(Nc,dimR,J,invJ,work,resNeg,guess);CHKERRQ(ierr); 2704 } 2705 } 2706 ierr = DMRestoreWorkArray(dm,3 * Nc * dimR,PETSC_SCALAR,&J);CHKERRQ(ierr); 2707 ierr = DMRestoreWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,PETSC_REAL,&B);CHKERRQ(ierr); 2708 ierr = DMRestoreWorkArray(dm,pdim,PETSC_REAL,&modes);CHKERRQ(ierr); 2709 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2710 PetscFunctionReturn(0); 2711 } 2712 2713 /* TODO: TOBY please fix this for Nc > 1 */ 2714 static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR) 2715 { 2716 PetscInt numComp, pdim, i, j, k, l, coordSize; 2717 PetscScalar *nodes = NULL; 2718 PetscReal *invV, *modes; 2719 PetscReal *B; 2720 PetscErrorCode ierr; 2721 2722 PetscFunctionBegin; 2723 ierr = PetscFEGetDimension(fe, &pdim);CHKERRQ(ierr); 2724 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 2725 if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc); 2726 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2727 /* convert nodes to values in the stable evaluation basis */ 2728 ierr = DMGetWorkArray(dm,pdim,PETSC_REAL,&modes);CHKERRQ(ierr); 2729 invV = fe->invV; 2730 for (i = 0; i < pdim; ++i) { 2731 modes[i] = 0.; 2732 for (j = 0; j < pdim; ++j) { 2733 modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]); 2734 } 2735 } 2736 ierr = DMGetWorkArray(dm,numPoints * pdim * Nc,PETSC_REAL,&B);CHKERRQ(ierr); 2737 ierr = PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL);CHKERRQ(ierr); 2738 for (i = 0; i < numPoints * Nc; i++) {realCoords[i] = 0.;} 2739 for (j = 0; j < numPoints; j++) { 2740 PetscReal *mapped = &realCoords[j * Nc]; 2741 2742 for (k = 0; k < pdim; k++) { 2743 for (l = 0; l < Nc; l++) { 2744 mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l]; 2745 } 2746 } 2747 } 2748 ierr = DMRestoreWorkArray(dm,numPoints * pdim * Nc,PETSC_REAL,&B);CHKERRQ(ierr); 2749 ierr = DMRestoreWorkArray(dm,pdim,PETSC_REAL,&modes);CHKERRQ(ierr); 2750 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2751 PetscFunctionReturn(0); 2752 } 2753 2754 /*@ 2755 DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element 2756 map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not 2757 extend uniquely outside the reference cell (e.g, most non-affine maps) 2758 2759 Not collective 2760 2761 Input Parameters: 2762 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 2763 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 2764 as a multilinear map for tensor-product elements 2765 . cell - the cell whose map is used. 2766 . numPoints - the number of points to locate 2767 - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 2768 2769 Output Parameters: 2770 . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 2771 2772 Level: intermediate 2773 @*/ 2774 PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[]) 2775 { 2776 PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i; 2777 DM coordDM = NULL; 2778 Vec coords; 2779 PetscFE fe = NULL; 2780 PetscErrorCode ierr; 2781 2782 PetscFunctionBegin; 2783 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2784 ierr = DMGetDimension(dm,&dimR);CHKERRQ(ierr); 2785 ierr = DMGetCoordinateDim(dm,&dimC);CHKERRQ(ierr); 2786 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 2787 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 2788 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2789 ierr = DMGetCoordinateDM(dm,&coordDM);CHKERRQ(ierr); 2790 if (coordDM) { 2791 PetscInt coordFields; 2792 2793 ierr = DMGetNumFields(coordDM,&coordFields);CHKERRQ(ierr); 2794 if (coordFields) { 2795 PetscClassId id; 2796 PetscObject disc; 2797 2798 ierr = DMGetField(coordDM,0,&disc);CHKERRQ(ierr); 2799 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 2800 if (id == PETSCFE_CLASSID) { 2801 fe = (PetscFE) disc; 2802 } 2803 } 2804 } 2805 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 2806 ierr = DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 2807 cEnd = cEndInterior > 0 ? cEndInterior : cEnd; 2808 if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);CHKERRQ(ierr); 2809 if (!fe) { /* implicit discretization: affine or multilinear */ 2810 PetscInt coneSize; 2811 PetscBool isSimplex, isTensor; 2812 2813 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 2814 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 2815 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 2816 if (isSimplex) { 2817 PetscReal detJ, *v0, *J, *invJ; 2818 2819 ierr = DMGetWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2820 J = &v0[dimC]; 2821 invJ = &J[dimC * dimC]; 2822 ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);CHKERRQ(ierr); 2823 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 2824 CoordinatesRealToRef(dimC, dimR, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]); 2825 } 2826 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2827 } else if (isTensor) { 2828 ierr = DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);CHKERRQ(ierr); 2829 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize); 2830 } else { 2831 ierr = DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);CHKERRQ(ierr); 2832 } 2833 PetscFunctionReturn(0); 2834 } 2835 2836 /*@ 2837 DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map. 2838 2839 Not collective 2840 2841 Input Parameters: 2842 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 2843 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 2844 as a multilinear map for tensor-product elements 2845 . cell - the cell whose map is used. 2846 . numPoints - the number of points to locate 2847 + refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 2848 2849 Output Parameters: 2850 . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 2851 2852 Level: intermediate 2853 @*/ 2854 PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[]) 2855 { 2856 PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i; 2857 DM coordDM = NULL; 2858 Vec coords; 2859 PetscFE fe = NULL; 2860 PetscErrorCode ierr; 2861 2862 PetscFunctionBegin; 2863 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2864 ierr = DMGetDimension(dm,&dimR);CHKERRQ(ierr); 2865 ierr = DMGetCoordinateDim(dm,&dimC);CHKERRQ(ierr); 2866 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 2867 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 2868 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2869 ierr = DMGetCoordinateDM(dm,&coordDM);CHKERRQ(ierr); 2870 if (coordDM) { 2871 PetscInt coordFields; 2872 2873 ierr = DMGetNumFields(coordDM,&coordFields);CHKERRQ(ierr); 2874 if (coordFields) { 2875 PetscClassId id; 2876 PetscObject disc; 2877 2878 ierr = DMGetField(coordDM,0,&disc);CHKERRQ(ierr); 2879 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 2880 if (id == PETSCFE_CLASSID) { 2881 fe = (PetscFE) disc; 2882 } 2883 } 2884 } 2885 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 2886 ierr = DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 2887 cEnd = cEndInterior > 0 ? cEndInterior : cEnd; 2888 if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);CHKERRQ(ierr); 2889 if (!fe) { /* implicit discretization: affine or multilinear */ 2890 PetscInt coneSize; 2891 PetscBool isSimplex, isTensor; 2892 2893 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 2894 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 2895 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 2896 if (isSimplex) { 2897 PetscReal detJ, *v0, *J; 2898 2899 ierr = DMGetWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2900 J = &v0[dimC]; 2901 ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ);CHKERRQ(ierr); 2902 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 2903 CoordinatesRefToReal(dimC, dimR, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]); 2904 } 2905 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2906 } else if (isTensor) { 2907 ierr = DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);CHKERRQ(ierr); 2908 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize); 2909 } else { 2910 ierr = DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);CHKERRQ(ierr); 2911 } 2912 PetscFunctionReturn(0); 2913 } 2914