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