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]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 829 if (PetscAbsReal(x2p[2]) > 1.0e-9) 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, PetscReal v0[], 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 *detJ = 0.0; 1102 if (numCoords == 12) { 1103 const PetscInt dim = 3; 1104 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 1105 1106 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1107 ierr = DMPlexComputeProjection3Dto2D(numCoords, coords, R);CHKERRQ(ierr); 1108 if (J) { 1109 const PetscInt pdim = 2; 1110 1111 for (d = 0; d < pdim; d++) { 1112 J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 1113 J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 1114 } 1115 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1116 DMPlex_Det3D_Internal(detJ, J0); 1117 for (d = 0; d < dim; d++) { 1118 for (f = 0; f < dim; f++) { 1119 J[d*dim+f] = 0.0; 1120 for (g = 0; g < dim; g++) { 1121 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 1122 } 1123 } 1124 } 1125 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1126 } 1127 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1128 } else if (numCoords == 8) { 1129 const PetscInt dim = 2; 1130 1131 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1132 if (J) { 1133 for (d = 0; d < dim; d++) { 1134 J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1135 J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1136 } 1137 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 1138 DMPlex_Det2D_Internal(detJ, J); 1139 } 1140 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 1141 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords); 1142 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 1143 PetscFunctionReturn(0); 1144 } 1145 1146 #undef __FUNCT__ 1147 #define __FUNCT__ "DMPlexComputeTetrahedronGeometry_Internal" 1148 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1149 { 1150 PetscSection coordSection; 1151 Vec coordinates; 1152 PetscScalar *coords = NULL; 1153 const PetscInt dim = 3; 1154 PetscInt d; 1155 PetscErrorCode ierr; 1156 1157 PetscFunctionBegin; 1158 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1159 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1160 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1161 *detJ = 0.0; 1162 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1163 if (J) { 1164 for (d = 0; d < dim; d++) { 1165 /* I orient with outward face normals */ 1166 J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d])); 1167 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1168 J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1169 } 1170 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1171 DMPlex_Det3D_Internal(detJ, J); 1172 } 1173 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1174 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1175 PetscFunctionReturn(0); 1176 } 1177 1178 #undef __FUNCT__ 1179 #define __FUNCT__ "DMPlexComputeHexahedronGeometry_Internal" 1180 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1181 { 1182 PetscSection coordSection; 1183 Vec coordinates; 1184 PetscScalar *coords = NULL; 1185 const PetscInt dim = 3; 1186 PetscInt d; 1187 PetscErrorCode ierr; 1188 1189 PetscFunctionBegin; 1190 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1191 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1192 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1193 *detJ = 0.0; 1194 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 1195 if (J) { 1196 for (d = 0; d < dim; d++) { 1197 J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 1198 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 1199 J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d])); 1200 } 1201 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 1202 DMPlex_Det3D_Internal(detJ, J); 1203 } 1204 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1205 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 1206 PetscFunctionReturn(0); 1207 } 1208 1209 #undef __FUNCT__ 1210 #define __FUNCT__ "DMPlexComputeCellGeometryAffineFEM" 1211 /*@C 1212 DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell 1213 1214 Collective on DM 1215 1216 Input Arguments: 1217 + dm - the DM 1218 - cell - the cell 1219 1220 Output Arguments: 1221 + v0 - the translation part of this affine transform 1222 . J - the Jacobian of the transform from the reference element 1223 . invJ - the inverse of the Jacobian 1224 - detJ - the Jacobian determinant 1225 1226 Level: advanced 1227 1228 Fortran Notes: 1229 Since it returns arrays, this routine is only available in Fortran 90, and you must 1230 include petsc.h90 in your code. 1231 1232 .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinateVec() 1233 @*/ 1234 PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1235 { 1236 PetscInt depth, dim, coneSize; 1237 DMLabel depthLabel; 1238 PetscErrorCode ierr; 1239 1240 PetscFunctionBegin; 1241 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1242 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 1243 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 1244 ierr = DMLabelGetValue(depthLabel, cell, &dim);CHKERRQ(ierr); 1245 if (depth == 1 && dim == 1) { 1246 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1247 } 1248 switch (dim) { 1249 case 0: 1250 ierr = DMPlexComputePointGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1251 break; 1252 case 1: 1253 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1254 break; 1255 case 2: 1256 switch (coneSize) { 1257 case 3: 1258 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1259 break; 1260 case 4: 1261 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1262 break; 1263 default: 1264 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 1265 } 1266 break; 1267 case 3: 1268 switch (coneSize) { 1269 case 4: 1270 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1271 break; 1272 case 6: /* Faces */ 1273 case 8: /* Vertices */ 1274 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1275 break; 1276 default: 1277 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 1278 } 1279 break; 1280 default: 1281 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1282 } 1283 PetscFunctionReturn(0); 1284 } 1285 1286 #undef __FUNCT__ 1287 #define __FUNCT__ "DMPlexComputeIsoparametricGeometry_Internal" 1288 static PetscErrorCode DMPlexComputeIsoparametricGeometry_Internal(DM dm, PetscFE fe, PetscInt point, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1289 { 1290 PetscQuadrature quad; 1291 PetscSection coordSection; 1292 Vec coordinates; 1293 PetscScalar *coords = NULL; 1294 const PetscReal *quadPoints; 1295 PetscReal *basisDer; 1296 PetscInt dim, cdim, pdim, qdim, Nq, numCoords, d, q; 1297 PetscErrorCode ierr; 1298 1299 PetscFunctionBegin; 1300 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1301 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1302 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1303 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1304 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 1305 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1306 ierr = PetscFEGetDimension(fe, &pdim);CHKERRQ(ierr); 1307 ierr = PetscQuadratureGetData(quad, &qdim, &Nq, &quadPoints, NULL);CHKERRQ(ierr); 1308 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 1309 *detJ = 0.0; 1310 if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim); 1311 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); 1312 if (v0) {for (d = 0; d < cdim; d++) v0[d] = PetscRealPart(coords[d]);} 1313 if (J) { 1314 ierr = PetscMemzero(J, Nq*cdim*dim*sizeof(PetscReal));CHKERRQ(ierr); 1315 for (q = 0; q < Nq; ++q) { 1316 PetscInt i, j, k, c, r; 1317 1318 /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */ 1319 for (k = 0; k < pdim; ++k) 1320 for (j = 0; j < dim; ++j) 1321 for (i = 0; i < cdim; ++i) 1322 J[(q*cdim + i)*dim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]); 1323 ierr = PetscLogFlops(2.0*pdim*dim*cdim);CHKERRQ(ierr); 1324 if (cdim > dim) { 1325 for (c = dim; c < cdim; ++c) 1326 for (r = 0; r < cdim; ++r) 1327 J[r*cdim+c] = r == c ? 1.0 : 0.0; 1328 } 1329 switch (cdim) { 1330 case 3: 1331 DMPlex_Det3D_Internal(detJ, J); 1332 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1333 break; 1334 case 2: 1335 DMPlex_Det2D_Internal(detJ, J); 1336 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 1337 break; 1338 case 1: 1339 *detJ = J[0]; 1340 if (invJ) invJ[0] = 1.0/J[0]; 1341 } 1342 } 1343 } 1344 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1345 PetscFunctionReturn(0); 1346 } 1347 1348 #undef __FUNCT__ 1349 #define __FUNCT__ "DMPlexComputeCellGeometryFEM" 1350 /*@C 1351 DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell 1352 1353 Collective on DM 1354 1355 Input Arguments: 1356 + dm - the DM 1357 . cell - the cell 1358 - fe - the finite element containing the quadrature 1359 1360 Output Arguments: 1361 + v0 - the translation part of this transform 1362 . J - the Jacobian of the transform from the reference element at each quadrature point 1363 . invJ - the inverse of the Jacobian at each quadrature point 1364 - detJ - the Jacobian determinant at each quadrature point 1365 1366 Level: advanced 1367 1368 Fortran Notes: 1369 Since it returns arrays, this routine is only available in Fortran 90, and you must 1370 include petsc.h90 in your code. 1371 1372 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1373 @*/ 1374 PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscFE fe, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1375 { 1376 PetscErrorCode ierr; 1377 1378 PetscFunctionBegin; 1379 if (!fe) {ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);} 1380 else {ierr = DMPlexComputeIsoparametricGeometry_Internal(dm, fe, cell, v0, J, invJ, detJ);CHKERRQ(ierr);} 1381 PetscFunctionReturn(0); 1382 } 1383 1384 #undef __FUNCT__ 1385 #define __FUNCT__ "DMPlexComputeGeometryFVM_1D_Internal" 1386 static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1387 { 1388 PetscSection coordSection; 1389 Vec coordinates; 1390 PetscScalar *coords = NULL; 1391 PetscScalar tmp[2]; 1392 PetscInt coordSize; 1393 PetscErrorCode ierr; 1394 1395 PetscFunctionBegin; 1396 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1397 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1398 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1399 if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now"); 1400 ierr = DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);CHKERRQ(ierr); 1401 if (centroid) { 1402 centroid[0] = 0.5*PetscRealPart(coords[0] + tmp[0]); 1403 centroid[1] = 0.5*PetscRealPart(coords[1] + tmp[1]); 1404 } 1405 if (normal) { 1406 PetscReal norm; 1407 1408 normal[0] = -PetscRealPart(coords[1] - tmp[1]); 1409 normal[1] = PetscRealPart(coords[0] - tmp[0]); 1410 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]); 1411 normal[0] /= norm; 1412 normal[1] /= norm; 1413 } 1414 if (vol) { 1415 *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - tmp[0])) + PetscSqr(PetscRealPart(coords[1] - tmp[1]))); 1416 } 1417 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1418 PetscFunctionReturn(0); 1419 } 1420 1421 #undef __FUNCT__ 1422 #define __FUNCT__ "DMPlexComputeGeometryFVM_2D_Internal" 1423 /* Centroid_i = (\sum_n A_n Cn_i ) / A */ 1424 static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1425 { 1426 PetscSection coordSection; 1427 Vec coordinates; 1428 PetscScalar *coords = NULL; 1429 PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9]; 1430 PetscInt tdim = 2, coordSize, numCorners, p, d, e; 1431 PetscErrorCode ierr; 1432 1433 PetscFunctionBegin; 1434 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1435 ierr = DMPlexGetConeSize(dm, cell, &numCorners);CHKERRQ(ierr); 1436 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1437 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1438 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 1439 if (dim > 2 && centroid) { 1440 v0[0] = PetscRealPart(coords[0]); 1441 v0[1] = PetscRealPart(coords[1]); 1442 v0[2] = PetscRealPart(coords[2]); 1443 } 1444 if (normal) { 1445 if (dim > 2) { 1446 const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]); 1447 const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]); 1448 const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]); 1449 PetscReal norm; 1450 1451 normal[0] = y0*z1 - z0*y1; 1452 normal[1] = z0*x1 - x0*z1; 1453 normal[2] = x0*y1 - y0*x1; 1454 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]); 1455 normal[0] /= norm; 1456 normal[1] /= norm; 1457 normal[2] /= norm; 1458 } else { 1459 for (d = 0; d < dim; ++d) normal[d] = 0.0; 1460 } 1461 } 1462 if (dim == 3) {ierr = DMPlexComputeProjection3Dto2D(coordSize, coords, R);CHKERRQ(ierr);} 1463 for (p = 0; p < numCorners; ++p) { 1464 /* Need to do this copy to get types right */ 1465 for (d = 0; d < tdim; ++d) { 1466 ctmp[d] = PetscRealPart(coords[p*tdim+d]); 1467 ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]); 1468 } 1469 Volume_Triangle_Origin_Internal(&vtmp, ctmp); 1470 vsum += vtmp; 1471 for (d = 0; d < tdim; ++d) { 1472 csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp; 1473 } 1474 } 1475 for (d = 0; d < tdim; ++d) { 1476 csum[d] /= (tdim+1)*vsum; 1477 } 1478 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1479 if (vol) *vol = PetscAbsReal(vsum); 1480 if (centroid) { 1481 if (dim > 2) { 1482 for (d = 0; d < dim; ++d) { 1483 centroid[d] = v0[d]; 1484 for (e = 0; e < dim; ++e) { 1485 centroid[d] += R[d*dim+e]*csum[e]; 1486 } 1487 } 1488 } else for (d = 0; d < dim; ++d) centroid[d] = csum[d]; 1489 } 1490 PetscFunctionReturn(0); 1491 } 1492 1493 #undef __FUNCT__ 1494 #define __FUNCT__ "DMPlexComputeGeometryFVM_3D_Internal" 1495 /* Centroid_i = (\sum_n V_n Cn_i ) / V */ 1496 static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1497 { 1498 PetscSection coordSection; 1499 Vec coordinates; 1500 PetscScalar *coords = NULL; 1501 PetscReal vsum = 0.0, vtmp, coordsTmp[3*3]; 1502 const PetscInt *faces, *facesO; 1503 PetscInt numFaces, f, coordSize, numCorners, p, d; 1504 PetscErrorCode ierr; 1505 1506 PetscFunctionBegin; 1507 if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim); 1508 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1509 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1510 1511 if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0; 1512 ierr = DMPlexGetConeSize(dm, cell, &numFaces);CHKERRQ(ierr); 1513 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 1514 ierr = DMPlexGetConeOrientation(dm, cell, &facesO);CHKERRQ(ierr); 1515 for (f = 0; f < numFaces; ++f) { 1516 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 1517 numCorners = coordSize/dim; 1518 switch (numCorners) { 1519 case 3: 1520 for (d = 0; d < dim; ++d) { 1521 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 1522 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 1523 coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]); 1524 } 1525 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1526 if (facesO[f] < 0) vtmp = -vtmp; 1527 vsum += vtmp; 1528 if (centroid) { /* Centroid of OABC = (a+b+c)/4 */ 1529 for (d = 0; d < dim; ++d) { 1530 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1531 } 1532 } 1533 break; 1534 case 4: 1535 /* DO FOR PYRAMID */ 1536 /* First tet */ 1537 for (d = 0; d < dim; ++d) { 1538 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 1539 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 1540 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 1541 } 1542 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1543 if (facesO[f] < 0) vtmp = -vtmp; 1544 vsum += vtmp; 1545 if (centroid) { 1546 for (d = 0; d < dim; ++d) { 1547 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1548 } 1549 } 1550 /* Second tet */ 1551 for (d = 0; d < dim; ++d) { 1552 coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]); 1553 coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]); 1554 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 1555 } 1556 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1557 if (facesO[f] < 0) vtmp = -vtmp; 1558 vsum += vtmp; 1559 if (centroid) { 1560 for (d = 0; d < dim; ++d) { 1561 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1562 } 1563 } 1564 break; 1565 default: 1566 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners); 1567 } 1568 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 1569 } 1570 if (vol) *vol = PetscAbsReal(vsum); 1571 if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0; 1572 if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4); 1573 PetscFunctionReturn(0); 1574 } 1575 1576 #undef __FUNCT__ 1577 #define __FUNCT__ "DMPlexComputeCellGeometryFVM" 1578 /*@C 1579 DMPlexComputeCellGeometryFVM - Compute the volume for a given cell 1580 1581 Collective on DM 1582 1583 Input Arguments: 1584 + dm - the DM 1585 - cell - the cell 1586 1587 Output Arguments: 1588 + volume - the cell volume 1589 . centroid - the cell centroid 1590 - normal - the cell normal, if appropriate 1591 1592 Level: advanced 1593 1594 Fortran Notes: 1595 Since it returns arrays, this routine is only available in Fortran 90, and you must 1596 include petsc.h90 in your code. 1597 1598 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1599 @*/ 1600 PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1601 { 1602 PetscInt depth, dim; 1603 PetscErrorCode ierr; 1604 1605 PetscFunctionBegin; 1606 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1607 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1608 if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 1609 /* We need to keep a pointer to the depth label */ 1610 ierr = DMGetLabelValue(dm, "depth", cell, &depth);CHKERRQ(ierr); 1611 /* Cone size is now the number of faces */ 1612 switch (depth) { 1613 case 1: 1614 ierr = DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1615 break; 1616 case 2: 1617 ierr = DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1618 break; 1619 case 3: 1620 ierr = DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1621 break; 1622 default: 1623 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1624 } 1625 PetscFunctionReturn(0); 1626 } 1627 1628 #undef __FUNCT__ 1629 #define __FUNCT__ "DMPlexComputeGeometryFEM" 1630 /* This should also take a PetscFE argument I think */ 1631 PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom) 1632 { 1633 DM dmCell; 1634 Vec coordinates; 1635 PetscSection coordSection, sectionCell; 1636 PetscScalar *cgeom; 1637 PetscInt cStart, cEnd, cMax, c; 1638 PetscErrorCode ierr; 1639 1640 PetscFunctionBegin; 1641 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1642 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1643 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1644 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1645 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1646 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1647 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1648 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1649 cEnd = cMax < 0 ? cEnd : cMax; 1650 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1651 /* TODO This needs to be multiplied by Nq for non-affine */ 1652 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFECellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1653 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1654 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1655 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1656 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1657 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1658 for (c = cStart; c < cEnd; ++c) { 1659 PetscFECellGeom *cg; 1660 1661 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1662 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1663 ierr = DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v0, cg->J, cg->invJ, &cg->detJ);CHKERRQ(ierr); 1664 if (cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c); 1665 } 1666 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1667 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1668 PetscFunctionReturn(0); 1669 } 1670 1671 #undef __FUNCT__ 1672 #define __FUNCT__ "DMPlexComputeGeometryFVM" 1673 /*@ 1674 DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method 1675 1676 Input Parameter: 1677 . dm - The DM 1678 1679 Output Parameters: 1680 + cellgeom - A Vec of PetscFVCellGeom data 1681 . facegeom - A Vec of PetscFVFaceGeom data 1682 1683 Level: developer 1684 1685 .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM() 1686 @*/ 1687 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 1688 { 1689 DM dmFace, dmCell; 1690 DMLabel ghostLabel; 1691 PetscSection sectionFace, sectionCell; 1692 PetscSection coordSection; 1693 Vec coordinates; 1694 PetscScalar *fgeom, *cgeom; 1695 PetscReal minradius, gminradius; 1696 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 1697 PetscErrorCode ierr; 1698 1699 PetscFunctionBegin; 1700 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1701 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1702 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1703 /* Make cell centroids and volumes */ 1704 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1705 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1706 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1707 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1708 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1709 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1710 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1711 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1712 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1713 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1714 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1715 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1716 if (cEndInterior < 0) { 1717 cEndInterior = cEnd; 1718 } 1719 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1720 for (c = cStart; c < cEndInterior; ++c) { 1721 PetscFVCellGeom *cg; 1722 1723 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1724 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1725 ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr); 1726 } 1727 /* Compute face normals and minimum cell radius */ 1728 ierr = DMClone(dm, &dmFace);CHKERRQ(ierr); 1729 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);CHKERRQ(ierr); 1730 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1731 ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr); 1732 for (f = fStart; f < fEnd; ++f) {ierr = PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1733 ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr); 1734 ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr); 1735 ierr = PetscSectionDestroy(§ionFace);CHKERRQ(ierr); 1736 ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr); 1737 ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr); 1738 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1739 minradius = PETSC_MAX_REAL; 1740 for (f = fStart; f < fEnd; ++f) { 1741 PetscFVFaceGeom *fg; 1742 PetscReal area; 1743 PetscInt ghost = -1, d, numChildren; 1744 1745 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1746 ierr = DMPlexGetTreeChildren(dm,f,&numChildren,NULL);CHKERRQ(ierr); 1747 if (ghost >= 0 || numChildren) continue; 1748 ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr); 1749 ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);CHKERRQ(ierr); 1750 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 1751 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 1752 { 1753 PetscFVCellGeom *cL, *cR; 1754 PetscInt ncells; 1755 const PetscInt *cells; 1756 PetscReal *lcentroid, *rcentroid; 1757 PetscReal l[3], r[3], v[3]; 1758 1759 ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr); 1760 ierr = DMPlexGetSupportSize(dm, f, &ncells);CHKERRQ(ierr); 1761 ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr); 1762 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 1763 if (ncells > 1) { 1764 ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr); 1765 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 1766 } 1767 else { 1768 rcentroid = fg->centroid; 1769 } 1770 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);CHKERRQ(ierr); 1771 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);CHKERRQ(ierr); 1772 DMPlex_WaxpyD_Internal(dim, -1, l, r, v); 1773 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 1774 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 1775 } 1776 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 1777 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]); 1778 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]); 1779 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f); 1780 } 1781 if (cells[0] < cEndInterior) { 1782 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 1783 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1784 } 1785 if (ncells > 1 && cells[1] < cEndInterior) { 1786 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 1787 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1788 } 1789 } 1790 } 1791 ierr = MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1792 ierr = DMPlexSetMinRadius(dm, gminradius);CHKERRQ(ierr); 1793 /* Compute centroids of ghost cells */ 1794 for (c = cEndInterior; c < cEnd; ++c) { 1795 PetscFVFaceGeom *fg; 1796 const PetscInt *cone, *support; 1797 PetscInt coneSize, supportSize, s; 1798 1799 ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr); 1800 if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize); 1801 ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr); 1802 ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr); 1803 if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize); 1804 ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr); 1805 ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr); 1806 for (s = 0; s < 2; ++s) { 1807 /* Reflect ghost centroid across plane of face */ 1808 if (support[s] == c) { 1809 PetscFVCellGeom *ci; 1810 PetscFVCellGeom *cg; 1811 PetscReal c2f[3], a; 1812 1813 ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr); 1814 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 1815 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 1816 ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr); 1817 DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid); 1818 cg->volume = ci->volume; 1819 } 1820 } 1821 } 1822 ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr); 1823 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1824 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1825 ierr = DMDestroy(&dmFace);CHKERRQ(ierr); 1826 PetscFunctionReturn(0); 1827 } 1828 1829 #undef __FUNCT__ 1830 #define __FUNCT__ "DMPlexGetMinRadius" 1831 /*@C 1832 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 1833 1834 Not collective 1835 1836 Input Argument: 1837 . dm - the DM 1838 1839 Output Argument: 1840 . minradius - the minium cell radius 1841 1842 Level: developer 1843 1844 .seealso: DMGetCoordinates() 1845 @*/ 1846 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 1847 { 1848 PetscFunctionBegin; 1849 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1850 PetscValidPointer(minradius,2); 1851 *minradius = ((DM_Plex*) dm->data)->minradius; 1852 PetscFunctionReturn(0); 1853 } 1854 1855 #undef __FUNCT__ 1856 #define __FUNCT__ "DMPlexSetMinRadius" 1857 /*@C 1858 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 1859 1860 Logically collective 1861 1862 Input Arguments: 1863 + dm - the DM 1864 - minradius - the minium cell radius 1865 1866 Level: developer 1867 1868 .seealso: DMSetCoordinates() 1869 @*/ 1870 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 1871 { 1872 PetscFunctionBegin; 1873 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1874 ((DM_Plex*) dm->data)->minradius = minradius; 1875 PetscFunctionReturn(0); 1876 } 1877 1878 #undef __FUNCT__ 1879 #define __FUNCT__ "BuildGradientReconstruction_Internal" 1880 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 1881 { 1882 DMLabel ghostLabel; 1883 PetscScalar *dx, *grad, **gref; 1884 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 1885 PetscErrorCode ierr; 1886 1887 PetscFunctionBegin; 1888 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1889 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1890 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1891 ierr = DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);CHKERRQ(ierr); 1892 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 1893 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1894 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 1895 for (c = cStart; c < cEndInterior; c++) { 1896 const PetscInt *faces; 1897 PetscInt numFaces, usedFaces, f, d; 1898 PetscFVCellGeom *cg; 1899 PetscBool boundary; 1900 PetscInt ghost; 1901 1902 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1903 ierr = DMPlexGetConeSize(dm, c, &numFaces);CHKERRQ(ierr); 1904 ierr = DMPlexGetCone(dm, c, &faces);CHKERRQ(ierr); 1905 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 1906 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 1907 PetscFVCellGeom *cg1; 1908 PetscFVFaceGeom *fg; 1909 const PetscInt *fcells; 1910 PetscInt ncell, side; 1911 1912 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 1913 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 1914 if ((ghost >= 0) || boundary) continue; 1915 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr); 1916 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 1917 ncell = fcells[!side]; /* the neighbor */ 1918 ierr = DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr); 1919 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 1920 for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d]; 1921 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 1922 } 1923 if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 1924 ierr = PetscFVComputeGradient(fvm, usedFaces, dx, grad);CHKERRQ(ierr); 1925 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 1926 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 1927 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 1928 if ((ghost >= 0) || boundary) continue; 1929 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d]; 1930 ++usedFaces; 1931 } 1932 } 1933 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 1934 PetscFunctionReturn(0); 1935 } 1936 1937 #undef __FUNCT__ 1938 #define __FUNCT__ "BuildGradientReconstruction_Internal_Tree" 1939 static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 1940 { 1941 DMLabel ghostLabel; 1942 PetscScalar *dx, *grad, **gref; 1943 PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0; 1944 PetscSection neighSec; 1945 PetscInt (*neighbors)[2]; 1946 PetscInt *counter; 1947 PetscErrorCode ierr; 1948 1949 PetscFunctionBegin; 1950 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1951 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1952 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1953 if (cEndInterior < 0) { 1954 cEndInterior = cEnd; 1955 } 1956 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);CHKERRQ(ierr); 1957 ierr = PetscSectionSetChart(neighSec,cStart,cEndInterior);CHKERRQ(ierr); 1958 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1959 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1960 for (f = fStart; f < fEnd; f++) { 1961 const PetscInt *fcells; 1962 PetscBool boundary; 1963 PetscInt ghost = -1; 1964 PetscInt numChildren, numCells, c; 1965 1966 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1967 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 1968 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 1969 if ((ghost >= 0) || boundary || numChildren) continue; 1970 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 1971 if (numCells == 2) { 1972 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 1973 for (c = 0; c < 2; c++) { 1974 PetscInt cell = fcells[c]; 1975 1976 if (cell >= cStart && cell < cEndInterior) { 1977 ierr = PetscSectionAddDof(neighSec,cell,1);CHKERRQ(ierr); 1978 } 1979 } 1980 } 1981 } 1982 ierr = PetscSectionSetUp(neighSec);CHKERRQ(ierr); 1983 ierr = PetscSectionGetMaxDof(neighSec,&maxNumFaces);CHKERRQ(ierr); 1984 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 1985 nStart = 0; 1986 ierr = PetscSectionGetStorageSize(neighSec,&nEnd);CHKERRQ(ierr); 1987 ierr = PetscMalloc1((nEnd-nStart),&neighbors);CHKERRQ(ierr); 1988 ierr = PetscCalloc1((cEndInterior-cStart),&counter);CHKERRQ(ierr); 1989 for (f = fStart; f < fEnd; f++) { 1990 const PetscInt *fcells; 1991 PetscBool boundary; 1992 PetscInt ghost = -1; 1993 PetscInt numChildren, numCells, c; 1994 1995 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1996 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 1997 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 1998 if ((ghost >= 0) || boundary || numChildren) continue; 1999 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 2000 if (numCells == 2) { 2001 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 2002 for (c = 0; c < 2; c++) { 2003 PetscInt cell = fcells[c], off; 2004 2005 if (cell >= cStart && cell < cEndInterior) { 2006 ierr = PetscSectionGetOffset(neighSec,cell,&off);CHKERRQ(ierr); 2007 off += counter[cell - cStart]++; 2008 neighbors[off][0] = f; 2009 neighbors[off][1] = fcells[1 - c]; 2010 } 2011 } 2012 } 2013 } 2014 ierr = PetscFree(counter);CHKERRQ(ierr); 2015 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 2016 for (c = cStart; c < cEndInterior; c++) { 2017 PetscInt numFaces, f, d, off, ghost = -1; 2018 PetscFVCellGeom *cg; 2019 2020 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 2021 ierr = PetscSectionGetDof(neighSec, c, &numFaces);CHKERRQ(ierr); 2022 ierr = PetscSectionGetOffset(neighSec, c, &off);CHKERRQ(ierr); 2023 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, c, &ghost);CHKERRQ(ierr);} 2024 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); 2025 for (f = 0; f < numFaces; ++f) { 2026 PetscFVCellGeom *cg1; 2027 PetscFVFaceGeom *fg; 2028 const PetscInt *fcells; 2029 PetscInt ncell, side, nface; 2030 2031 nface = neighbors[off + f][0]; 2032 ncell = neighbors[off + f][1]; 2033 ierr = DMPlexGetSupport(dm,nface,&fcells);CHKERRQ(ierr); 2034 side = (c != fcells[0]); 2035 ierr = DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);CHKERRQ(ierr); 2036 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 2037 for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d]; 2038 gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */ 2039 } 2040 ierr = PetscFVComputeGradient(fvm, numFaces, dx, grad);CHKERRQ(ierr); 2041 for (f = 0; f < numFaces; ++f) { 2042 for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d]; 2043 } 2044 } 2045 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 2046 ierr = PetscSectionDestroy(&neighSec);CHKERRQ(ierr); 2047 ierr = PetscFree(neighbors);CHKERRQ(ierr); 2048 PetscFunctionReturn(0); 2049 } 2050 2051 #undef __FUNCT__ 2052 #define __FUNCT__ "DMPlexComputeGradientFVM" 2053 /*@ 2054 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 2055 2056 Collective on DM 2057 2058 Input Arguments: 2059 + dm - The DM 2060 . fvm - The PetscFV 2061 . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM() 2062 - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM() 2063 2064 Output Parameters: 2065 + faceGeometry - The geometric factors for gradient calculation are inserted 2066 - dmGrad - The DM describing the layout of gradient data 2067 2068 Level: developer 2069 2070 .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM() 2071 @*/ 2072 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 2073 { 2074 DM dmFace, dmCell; 2075 PetscScalar *fgeom, *cgeom; 2076 PetscSection sectionGrad, parentSection; 2077 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 2078 PetscErrorCode ierr; 2079 2080 PetscFunctionBegin; 2081 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2082 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 2083 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2084 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2085 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 2086 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 2087 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 2088 ierr = VecGetArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2089 ierr = VecGetArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2090 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 2091 if (!parentSection) { 2092 ierr = BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2093 } else { 2094 ierr = BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 2095 } 2096 ierr = VecRestoreArray(faceGeometry, &fgeom);CHKERRQ(ierr); 2097 ierr = VecRestoreArray(cellGeometry, &cgeom);CHKERRQ(ierr); 2098 /* Create storage for gradients */ 2099 ierr = DMClone(dm, dmGrad);CHKERRQ(ierr); 2100 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 2101 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 2102 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 2103 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 2104 ierr = DMSetDefaultSection(*dmGrad, sectionGrad);CHKERRQ(ierr); 2105 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 2106 PetscFunctionReturn(0); 2107 } 2108 2109 #undef __FUNCT__ 2110 #define __FUNCT__ "DMPlexGetDataFVM" 2111 PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM) 2112 { 2113 PetscObject cellgeomobj, facegeomobj; 2114 PetscErrorCode ierr; 2115 2116 PetscFunctionBegin; 2117 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2118 if (!cellgeomobj) { 2119 Vec cellgeomInt, facegeomInt; 2120 2121 ierr = DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);CHKERRQ(ierr); 2122 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);CHKERRQ(ierr); 2123 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);CHKERRQ(ierr); 2124 ierr = VecDestroy(&cellgeomInt);CHKERRQ(ierr); 2125 ierr = VecDestroy(&facegeomInt);CHKERRQ(ierr); 2126 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 2127 } 2128 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);CHKERRQ(ierr); 2129 if (cellgeom) *cellgeom = (Vec) cellgeomobj; 2130 if (facegeom) *facegeom = (Vec) facegeomobj; 2131 if (gradDM) { 2132 PetscObject gradobj; 2133 PetscBool computeGradients; 2134 2135 ierr = PetscFVGetComputeGradients(fv,&computeGradients);CHKERRQ(ierr); 2136 if (!computeGradients) { 2137 *gradDM = NULL; 2138 PetscFunctionReturn(0); 2139 } 2140 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2141 if (!gradobj) { 2142 DM dmGradInt; 2143 2144 ierr = DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);CHKERRQ(ierr); 2145 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);CHKERRQ(ierr); 2146 ierr = DMDestroy(&dmGradInt);CHKERRQ(ierr); 2147 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 2148 } 2149 *gradDM = (DM) gradobj; 2150 } 2151 PetscFunctionReturn(0); 2152 } 2153 2154 #undef __FUNCT__ 2155 #define __FUNCT__ "DMPlexCoordinatesToReference_NewtonUpdate" 2156 static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess) 2157 { 2158 PetscFunctionBeginHot; 2159 2160 PetscInt l, m; 2161 2162 if (dimC == dimR && dimR <= 3) { 2163 /* invert Jacobian, multiply */ 2164 PetscScalar det, idet; 2165 2166 switch (dimR) { 2167 case 1: 2168 invJ[0] = 1./ J[0]; 2169 break; 2170 case 2: 2171 det = J[0] * J[3] - J[1] * J[2]; 2172 idet = 1./det; 2173 invJ[0] = J[3] * idet; 2174 invJ[1] = -J[1] * idet; 2175 invJ[2] = -J[2] * idet; 2176 invJ[3] = J[0] * idet; 2177 break; 2178 case 3: 2179 { 2180 invJ[0] = J[4] * J[8] - J[5] * J[7]; 2181 invJ[1] = J[2] * J[7] - J[1] * J[8]; 2182 invJ[2] = J[1] * J[5] - J[2] * J[4]; 2183 det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6]; 2184 idet = 1./det; 2185 invJ[0] *= idet; 2186 invJ[1] *= idet; 2187 invJ[2] *= idet; 2188 invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]); 2189 invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]); 2190 invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]); 2191 invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]); 2192 invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]); 2193 invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]); 2194 } 2195 break; 2196 } 2197 for (l = 0; l < dimR; l++) { 2198 for (m = 0; m < dimC; m++) { 2199 guess[m] += invJ[l * dimC + m] * resNeg[m]; 2200 } 2201 } 2202 } else { 2203 #if defined(PETSC_USE_COMPLEX) 2204 char transpose = 'C'; 2205 #else 2206 char transpose = 'T'; 2207 #endif 2208 PetscBLASInt m = dimR; 2209 PetscBLASInt n = dimC; 2210 PetscBLASInt one = 1; 2211 PetscBLASInt worksize = dimR * dimC, info; 2212 2213 for (l = 0; l < dimC; l++) {invJ[l] = resNeg[l];} 2214 2215 PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&m,&n,&one,J,&m,invJ,&n,work,&worksize, &info)); 2216 if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS"); 2217 2218 for (l = 0; l < dimR; l++) {guess[l] += invJ[l];} 2219 } 2220 PetscFunctionReturn(0); 2221 } 2222 2223 #undef __FUNCT__ 2224 #define __FUNCT__ "DMPlexCoordinatesToReference_Tensor" 2225 static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2226 { 2227 PetscInt coordSize, i, j, k, l, m, maxIts = 10, numV = (1 << dimR); 2228 PetscScalar *coordsScalar = NULL; 2229 PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg; 2230 PetscScalar *J, *invJ, *work; 2231 PetscErrorCode ierr; 2232 2233 PetscFunctionBegin; 2234 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2235 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2236 if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);CHKERRQ(ierr); 2237 ierr = DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, PETSC_REAL, &cellData);CHKERRQ(ierr); 2238 ierr = DMGetWorkArray(dm, 3 * dimR * dimC, PETSC_SCALAR, &J);CHKERRQ(ierr); 2239 cellCoords = &cellData[0]; 2240 cellCoeffs = &cellData[coordSize]; 2241 extJ = &cellData[2 * coordSize]; 2242 resNeg = &cellData[2 * coordSize + dimR]; 2243 invJ = &J[dimR * dimC]; 2244 work = &J[2 * dimR * dimC]; 2245 if (dimR == 2) { 2246 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 2247 2248 for (i = 0; i < 4; i++) { 2249 PetscInt plexI = zToPlex[i]; 2250 2251 for (j = 0; j < dimC; j++) { 2252 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2253 } 2254 } 2255 } else if (dimR == 3) { 2256 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2257 2258 for (i = 0; i < 8; i++) { 2259 PetscInt plexI = zToPlex[i]; 2260 2261 for (j = 0; j < dimC; j++) { 2262 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2263 } 2264 } 2265 } else { 2266 for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);} 2267 } 2268 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 2269 for (i = 0; i < dimR; i++) { 2270 PetscReal *swap; 2271 2272 for (j = 0; j < (numV / 2); j++) { 2273 for (k = 0; k < dimC; k++) { 2274 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 2275 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 2276 } 2277 } 2278 2279 if (i < dimR - 1) { 2280 swap = cellCoeffs; 2281 cellCoeffs = cellCoords; 2282 cellCoords = swap; 2283 } 2284 } 2285 ierr = PetscMemzero(refCoords,numPoints * dimR * sizeof (PetscReal));CHKERRQ(ierr); 2286 for (j = 0; j < numPoints; j++) { 2287 for (i = 0; i < maxIts; i++) { 2288 PetscReal *guess = &refCoords[dimR * j]; 2289 2290 /* compute -residual and Jacobian */ 2291 for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];} 2292 for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;} 2293 for (k = 0; k < numV; k++) { 2294 PetscReal extCoord = 1.; 2295 for (l = 0; l < dimR; l++) { 2296 PetscReal coord = guess[l]; 2297 PetscInt dep = (k & (1 << l)) >> l; 2298 2299 extCoord *= dep * coord + !dep; 2300 extJ[l] = dep; 2301 2302 for (m = 0; m < dimR; m++) { 2303 PetscReal coord = guess[m]; 2304 PetscInt dep = ((k & (1 << m)) >> m) && (m != l); 2305 PetscReal mult = dep * coord + !dep; 2306 2307 extJ[l] *= mult; 2308 } 2309 } 2310 for (l = 0; l < dimC; l++) { 2311 PetscReal coeff = cellCoeffs[dimC * k + l]; 2312 2313 resNeg[l] -= coeff * extCoord; 2314 for (m = 0; m < dimR; m++) { 2315 J[dimR * l + m] += coeff * extJ[m]; 2316 } 2317 } 2318 } 2319 2320 ierr = DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);CHKERRQ(ierr); 2321 } 2322 } 2323 ierr = DMRestoreWorkArray(dm, 3 * dimR * dimC, PETSC_SCALAR, &J);CHKERRQ(ierr); 2324 ierr = DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, PETSC_REAL, &cellData);CHKERRQ(ierr); 2325 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2326 PetscFunctionReturn(0); 2327 } 2328 2329 #undef __FUNCT__ 2330 #define __FUNCT__ "DMPlexReferenceToCoordinates_Tensor" 2331 static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2332 { 2333 PetscInt coordSize, i, j, k, l, numV = (1 << dimR); 2334 PetscScalar *coordsScalar = NULL; 2335 PetscReal *cellData, *cellCoords, *cellCoeffs; 2336 PetscErrorCode ierr; 2337 2338 PetscFunctionBegin; 2339 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2340 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2341 if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);CHKERRQ(ierr); 2342 ierr = DMGetWorkArray(dm, 2 * coordSize, PETSC_REAL, &cellData);CHKERRQ(ierr); 2343 cellCoords = &cellData[0]; 2344 cellCoeffs = &cellData[coordSize]; 2345 if (dimR == 2) { 2346 const PetscInt zToPlex[4] = {0, 1, 3, 2}; 2347 2348 for (i = 0; i < 4; i++) { 2349 PetscInt plexI = zToPlex[i]; 2350 2351 for (j = 0; j < dimC; j++) { 2352 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2353 } 2354 } 2355 } else if (dimR == 3) { 2356 const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6}; 2357 2358 for (i = 0; i < 8; i++) { 2359 PetscInt plexI = zToPlex[i]; 2360 2361 for (j = 0; j < dimC; j++) { 2362 cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]); 2363 } 2364 } 2365 } else { 2366 for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);} 2367 } 2368 /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */ 2369 for (i = 0; i < dimR; i++) { 2370 PetscReal *swap; 2371 2372 for (j = 0; j < (numV / 2); j++) { 2373 for (k = 0; k < dimC; k++) { 2374 cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]); 2375 cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]); 2376 } 2377 } 2378 2379 if (i < dimR - 1) { 2380 swap = cellCoeffs; 2381 cellCoeffs = cellCoords; 2382 cellCoords = swap; 2383 } 2384 } 2385 ierr = PetscMemzero(realCoords,numPoints * dimC * sizeof (PetscReal));CHKERRQ(ierr); 2386 for (j = 0; j < numPoints; j++) { 2387 const PetscReal *guess = &refCoords[dimR * j]; 2388 PetscReal *mapped = &realCoords[dimC * j]; 2389 2390 for (k = 0; k < numV; k++) { 2391 PetscReal extCoord = 1.; 2392 for (l = 0; l < dimR; l++) { 2393 PetscReal coord = guess[l]; 2394 PetscInt dep = (k & (1 << l)) >> l; 2395 2396 extCoord *= dep * coord + !dep; 2397 } 2398 for (l = 0; l < dimC; l++) { 2399 PetscReal coeff = cellCoeffs[dimC * k + l]; 2400 2401 mapped[l] += coeff * extCoord; 2402 } 2403 } 2404 } 2405 ierr = DMRestoreWorkArray(dm, 2 * coordSize, PETSC_REAL, &cellData);CHKERRQ(ierr); 2406 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);CHKERRQ(ierr); 2407 PetscFunctionReturn(0); 2408 } 2409 2410 #undef __FUNCT__ 2411 #define __FUNCT__ "DMPlexCoordinatesToReference_FE" 2412 static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2413 { 2414 PetscInt numComp, numDof, i, j, k, l, m, maxIter = 10, coordSize; 2415 PetscScalar *nodes, *modes, *invV; 2416 PetscReal *B; 2417 PetscReal *D; 2418 PetscScalar *resNeg, *J, *invJ, *work; 2419 PetscErrorCode ierr; 2420 2421 PetscFunctionBegin; 2422 ierr = PetscFEGetDimension(fe, &numDof);CHKERRQ(ierr); 2423 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 2424 if (numComp != dimC) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,dimC); 2425 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2426 /* convert nodes to values in the stable evaluation basis */ 2427 ierr = DMGetWorkArray(dm,dimC * numDof,PETSC_SCALAR,&modes);CHKERRQ(ierr); 2428 invV = fe->invV; 2429 for (i = 0; i < numDof; i++) { 2430 for (j = 0; j < dimC; j++) { 2431 modes[i * dimC + j] = 0.; 2432 for (k = 0; k < numDof; k++) { 2433 modes[i * dimC + j] += invV[i * numDof + k] * nodes[k * dimC + j]; 2434 } 2435 } 2436 } 2437 ierr = DMGetWorkArray(dm,numDof * numDof,PETSC_REAL,&B);CHKERRQ(ierr); 2438 ierr = DMGetWorkArray(dm,numDof * dimR,PETSC_REAL,&D);CHKERRQ(ierr); 2439 ierr = DMGetWorkArray(dm,dimC + 3 * dimC * dimR,PETSC_SCALAR,&resNeg);CHKERRQ(ierr); 2440 J = &resNeg[dimC]; 2441 invJ = &J[dimC * dimR]; 2442 work = &invJ[dimC * dimR]; 2443 for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;} 2444 for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */ 2445 for (j = 0; j < numPoints; j++) { 2446 PetscReal *guess = &refCoords[j * dimR]; 2447 ierr = PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);CHKERRQ(ierr); 2448 for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[j * dimC + k];} 2449 for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;} 2450 for (k = 0; k < numDof; k++) { 2451 for (l = 0; l < dimC; l++) { 2452 resNeg[l] -= modes[k * dimC + l] * B[k]; 2453 for (m = 0; m < dimR; m++) { 2454 J[l * dimR + m] += modes[k * dimC + l] * D[k * dimR + m]; 2455 } 2456 } 2457 } 2458 ierr = DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);CHKERRQ(ierr); 2459 } 2460 } 2461 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC + dimR,PETSC_SCALAR,&resNeg);CHKERRQ(ierr); 2462 ierr = DMRestoreWorkArray(dm,numDof * numDof * dimR,PETSC_REAL,&D);CHKERRQ(ierr); 2463 ierr = DMRestoreWorkArray(dm,numDof * numDof,PETSC_REAL,&B);CHKERRQ(ierr); 2464 ierr = DMRestoreWorkArray(dm,dimC * numDof,PETSC_SCALAR,&modes);CHKERRQ(ierr); 2465 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2466 PetscFunctionReturn(0); 2467 } 2468 2469 #undef __FUNCT__ 2470 #define __FUNCT__ "DMPlexReferenceToCoordinates_FE" 2471 static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR) 2472 { 2473 PetscInt numComp, numDof, i, j, k, l, coordSize; 2474 PetscScalar *nodes, *modes, *invV; 2475 PetscReal *B; 2476 PetscErrorCode ierr; 2477 2478 PetscFunctionBegin; 2479 ierr = PetscFEGetDimension(fe, &numDof);CHKERRQ(ierr); 2480 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 2481 if (numComp != dimC) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,dimC); 2482 ierr = DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2483 /* convert nodes to values in the stable evaluation basis */ 2484 ierr = DMGetWorkArray(dm,dimC * numDof,PETSC_SCALAR,&modes);CHKERRQ(ierr); 2485 invV = fe->invV; 2486 for (i = 0; i < numDof; i++) { 2487 for (j = 0; j < dimC; j++) { 2488 modes[i * dimC + j] = 0.; 2489 for (k = 0; k < numDof; k++) { 2490 modes[i * dimC + j] += invV[i * numDof + k] * nodes[k * dimC + j]; 2491 } 2492 } 2493 } 2494 ierr = DMGetWorkArray(dm,numDof * numDof,PETSC_REAL,&B);CHKERRQ(ierr); 2495 for (i = 0; i < numPoints * dimC; i++) {realCoords[i] = 0.;} 2496 for (j = 0; j < numPoints; j++) { 2497 const PetscReal *guess = &refCoords[j * dimR]; 2498 PetscReal *mapped = &realCoords[j * dimC]; 2499 2500 ierr = PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, NULL, NULL);CHKERRQ(ierr); 2501 for (k = 0; k < numDof; k++) { 2502 for (l = 0; l < dimC; l++) { 2503 mapped[l] += modes[k * dimC + l] * B[k]; 2504 } 2505 } 2506 } 2507 ierr = DMRestoreWorkArray(dm,numDof * numDof,PETSC_REAL,&B);CHKERRQ(ierr); 2508 ierr = DMRestoreWorkArray(dm,dimC * numDof,PETSC_SCALAR,&modes);CHKERRQ(ierr); 2509 ierr = DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);CHKERRQ(ierr); 2510 PetscFunctionReturn(0); 2511 } 2512 2513 #undef __FUNCT__ 2514 #define __FUNCT__ "DMPlexCoordinatesToReference" 2515 /*@ 2516 DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element 2517 map. This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not 2518 extend uniquely outside the reference cell (e.g, most non-affine maps) 2519 2520 Not collective 2521 2522 Input Parameters: 2523 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 2524 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 2525 as a multilinear map for tensor-product elements 2526 . cell - the cell whose map is used. 2527 . numPoints - the number of points to locate 2528 + realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 2529 2530 Output Parameters: 2531 . refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 2532 @*/ 2533 PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[]) 2534 { 2535 PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i; 2536 DM coordDM = NULL; 2537 Vec coords; 2538 PetscFE fe = NULL; 2539 PetscErrorCode ierr; 2540 2541 PetscFunctionBegin; 2542 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2543 ierr = DMGetDimension(dm,&dimR);CHKERRQ(ierr); 2544 ierr = DMGetCoordinateDim(dm,&dimC);CHKERRQ(ierr); 2545 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 2546 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 2547 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2548 ierr = DMGetCoordinateDM(dm,&coordDM);CHKERRQ(ierr); 2549 if (coordDM) { 2550 PetscInt coordFields; 2551 2552 ierr = DMGetNumFields(coordDM,&coordFields);CHKERRQ(ierr); 2553 if (coordFields) { 2554 PetscClassId id; 2555 PetscObject disc; 2556 2557 ierr = DMGetField(coordDM,0,&disc);CHKERRQ(ierr); 2558 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 2559 if (id == PETSCFE_CLASSID) { 2560 fe = (PetscFE) disc; 2561 } 2562 } 2563 } 2564 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 2565 ierr = DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 2566 cEnd = cEndInterior > 0 ? cEndInterior : cEnd; 2567 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); 2568 if (!fe) { /* implicit discretization: affine or multilinear */ 2569 PetscInt coneSize; 2570 PetscBool isSimplex, isTensor; 2571 2572 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 2573 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 2574 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 2575 if (isSimplex) { 2576 PetscReal detJ, *v0, *J, *invJ; 2577 2578 ierr = DMGetWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2579 J = &v0[dimC]; 2580 invJ = &J[dimC * dimC]; 2581 ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);CHKERRQ(ierr); 2582 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 2583 CoordinatesRealToRef(dimC, dimR, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]); 2584 } 2585 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2586 } else if (isTensor) { 2587 ierr = DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);CHKERRQ(ierr); 2588 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize); 2589 } else { 2590 ierr = DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);CHKERRQ(ierr); 2591 } 2592 PetscFunctionReturn(0); 2593 } 2594 2595 #undef __FUNCT__ 2596 #define __FUNCT__ "DMPlexReferenceToCoordinates" 2597 /*@ 2598 DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the the mesh for a single element map. 2599 2600 Not collective 2601 2602 Input Parameters: 2603 + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or 2604 implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or 2605 as a multilinear map for tensor-product elements 2606 . cell - the cell whose map is used. 2607 . numPoints - the number of points to locate 2608 + refCoords - (numPoints x dimension) array of reference coordinates (see DMGetDimension()) 2609 2610 Output Parameters: 2611 . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim()) 2612 @*/ 2613 PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[]) 2614 { 2615 PetscInt dimC, dimR, depth, cStart, cEnd, cEndInterior, i; 2616 DM coordDM = NULL; 2617 Vec coords; 2618 PetscFE fe = NULL; 2619 PetscErrorCode ierr; 2620 2621 PetscFunctionBegin; 2622 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2623 ierr = DMGetDimension(dm,&dimR);CHKERRQ(ierr); 2624 ierr = DMGetCoordinateDim(dm,&dimC);CHKERRQ(ierr); 2625 if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(0); 2626 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 2627 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2628 ierr = DMGetCoordinateDM(dm,&coordDM);CHKERRQ(ierr); 2629 if (coordDM) { 2630 PetscInt coordFields; 2631 2632 ierr = DMGetNumFields(coordDM,&coordFields);CHKERRQ(ierr); 2633 if (coordFields) { 2634 PetscClassId id; 2635 PetscObject disc; 2636 2637 ierr = DMGetField(coordDM,0,&disc);CHKERRQ(ierr); 2638 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 2639 if (id == PETSCFE_CLASSID) { 2640 fe = (PetscFE) disc; 2641 } 2642 } 2643 } 2644 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 2645 ierr = DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 2646 cEnd = cEndInterior > 0 ? cEndInterior : cEnd; 2647 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); 2648 if (!fe) { /* implicit discretization: affine or multilinear */ 2649 PetscInt coneSize; 2650 PetscBool isSimplex, isTensor; 2651 2652 ierr = DMPlexGetConeSize(dm,cell,&coneSize);CHKERRQ(ierr); 2653 isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE; 2654 isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE; 2655 if (isSimplex) { 2656 PetscReal detJ, *v0, *J; 2657 2658 ierr = DMGetWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2659 J = &v0[dimC]; 2660 ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ);CHKERRQ(ierr); 2661 for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */ 2662 CoordinatesRefToReal(dimC, dimR, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]); 2663 } 2664 ierr = DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, PETSC_REAL, &v0);CHKERRQ(ierr); 2665 } else if (isTensor) { 2666 ierr = DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);CHKERRQ(ierr); 2667 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize); 2668 } else { 2669 ierr = DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);CHKERRQ(ierr); 2670 } 2671 PetscFunctionReturn(0); 2672 } 2673