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