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