1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "DMPlexGetLineIntersection_2D_Internal" 5 static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection) 6 { 7 const PetscReal p0_x = segmentA[0*2+0]; 8 const PetscReal p0_y = segmentA[0*2+1]; 9 const PetscReal p1_x = segmentA[1*2+0]; 10 const PetscReal p1_y = segmentA[1*2+1]; 11 const PetscReal p2_x = segmentB[0*2+0]; 12 const PetscReal p2_y = segmentB[0*2+1]; 13 const PetscReal p3_x = segmentB[1*2+0]; 14 const PetscReal p3_y = segmentB[1*2+1]; 15 const PetscReal s1_x = p1_x - p0_x; 16 const PetscReal s1_y = p1_y - p0_y; 17 const PetscReal s2_x = p3_x - p2_x; 18 const PetscReal s2_y = p3_y - p2_y; 19 const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y); 20 21 PetscFunctionBegin; 22 *hasIntersection = PETSC_FALSE; 23 /* Non-parallel lines */ 24 if (denom != 0.0) { 25 const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom; 26 const PetscReal t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom; 27 28 if (s >= 0 && s <= 1 && t >= 0 && t <= 1) { 29 *hasIntersection = PETSC_TRUE; 30 if (intersection) { 31 intersection[0] = p0_x + (t * s1_x); 32 intersection[1] = p0_y + (t * s1_y); 33 } 34 } 35 } 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "DMPlexLocatePoint_Simplex_2D_Internal" 41 static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 42 { 43 const PetscInt embedDim = 2; 44 const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON; 45 PetscReal x = PetscRealPart(point[0]); 46 PetscReal y = PetscRealPart(point[1]); 47 PetscReal v0[2], J[4], invJ[4], detJ; 48 PetscReal xi, eta; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 53 xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]); 54 eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]); 55 56 if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0+eps)) *cell = c; 57 else *cell = -1; 58 PetscFunctionReturn(0); 59 } 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "DMPlexLocatePoint_General_2D_Internal" 63 static PetscErrorCode DMPlexLocatePoint_General_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 64 { 65 PetscSection coordSection; 66 Vec coordsLocal; 67 PetscScalar *coords = NULL; 68 const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0}; 69 PetscReal x = PetscRealPart(point[0]); 70 PetscReal y = PetscRealPart(point[1]); 71 PetscInt crossings = 0, f; 72 PetscErrorCode ierr; 73 74 PetscFunctionBegin; 75 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 76 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 77 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 78 for (f = 0; f < 4; ++f) { 79 PetscReal x_i = PetscRealPart(coords[faces[2*f+0]*2+0]); 80 PetscReal y_i = PetscRealPart(coords[faces[2*f+0]*2+1]); 81 PetscReal x_j = PetscRealPart(coords[faces[2*f+1]*2+0]); 82 PetscReal y_j = PetscRealPart(coords[faces[2*f+1]*2+1]); 83 PetscReal slope = (y_j - y_i) / (x_j - x_i); 84 PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE; 85 PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE; 86 PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE; 87 if ((cond1 || cond2) && above) ++crossings; 88 } 89 if (crossings % 2) *cell = c; 90 else *cell = -1; 91 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 92 PetscFunctionReturn(0); 93 } 94 95 #undef __FUNCT__ 96 #define __FUNCT__ "DMPlexLocatePoint_Simplex_3D_Internal" 97 static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 98 { 99 const PetscInt embedDim = 3; 100 PetscReal v0[3], J[9], invJ[9], detJ; 101 PetscReal x = PetscRealPart(point[0]); 102 PetscReal y = PetscRealPart(point[1]); 103 PetscReal z = PetscRealPart(point[2]); 104 PetscReal xi, eta, zeta; 105 PetscErrorCode ierr; 106 107 PetscFunctionBegin; 108 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 109 xi = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]); 110 eta = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]); 111 zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]); 112 113 if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c; 114 else *cell = -1; 115 PetscFunctionReturn(0); 116 } 117 118 #undef __FUNCT__ 119 #define __FUNCT__ "DMPlexLocatePoint_General_3D_Internal" 120 static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell) 121 { 122 PetscSection coordSection; 123 Vec coordsLocal; 124 PetscScalar *coords; 125 const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5, 126 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4}; 127 PetscBool found = PETSC_TRUE; 128 PetscInt f; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 133 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 134 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 135 for (f = 0; f < 6; ++f) { 136 /* Check the point is under plane */ 137 /* Get face normal */ 138 PetscReal v_i[3]; 139 PetscReal v_j[3]; 140 PetscReal normal[3]; 141 PetscReal pp[3]; 142 PetscReal dot; 143 144 v_i[0] = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]); 145 v_i[1] = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]); 146 v_i[2] = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]); 147 v_j[0] = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]); 148 v_j[1] = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]); 149 v_j[2] = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]); 150 normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1]; 151 normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2]; 152 normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0]; 153 pp[0] = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]); 154 pp[1] = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]); 155 pp[2] = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]); 156 dot = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2]; 157 158 /* Check that projected point is in face (2D location problem) */ 159 if (dot < 0.0) { 160 found = PETSC_FALSE; 161 break; 162 } 163 } 164 if (found) *cell = c; 165 else *cell = -1; 166 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);CHKERRQ(ierr); 167 PetscFunctionReturn(0); 168 } 169 170 #undef __FUNCT__ 171 #define __FUNCT__ "PetscGridHashInitialize_Internal" 172 static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[]) 173 { 174 PetscInt d; 175 176 PetscFunctionBegin; 177 box->dim = dim; 178 for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]); 179 PetscFunctionReturn(0); 180 } 181 182 #undef __FUNCT__ 183 #define __FUNCT__ "PetscGridHashCreate" 184 PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box) 185 { 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 ierr = PetscMalloc1(1, box);CHKERRQ(ierr); 190 ierr = PetscGridHashInitialize_Internal(*box, dim, point);CHKERRQ(ierr); 191 PetscFunctionReturn(0); 192 } 193 194 #undef __FUNCT__ 195 #define __FUNCT__ "PetscGridHashEnlarge" 196 PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[]) 197 { 198 PetscInt d; 199 200 PetscFunctionBegin; 201 for (d = 0; d < box->dim; ++d) { 202 box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d])); 203 box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d])); 204 } 205 PetscFunctionReturn(0); 206 } 207 208 #undef __FUNCT__ 209 #define __FUNCT__ "PetscGridHashSetGrid" 210 PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[]) 211 { 212 PetscInt d; 213 214 PetscFunctionBegin; 215 for (d = 0; d < box->dim; ++d) { 216 box->extent[d] = box->upper[d] - box->lower[d]; 217 if (n[d] == PETSC_DETERMINE) { 218 box->h[d] = h[d]; 219 box->n[d] = PetscCeilReal(box->extent[d]/h[d]); 220 } else { 221 box->n[d] = n[d]; 222 box->h[d] = box->extent[d]/n[d]; 223 } 224 } 225 PetscFunctionReturn(0); 226 } 227 228 #undef __FUNCT__ 229 #define __FUNCT__ "PetscGridHashGetEnclosingBox" 230 PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[]) 231 { 232 const PetscReal *lower = box->lower; 233 const PetscReal *upper = box->upper; 234 const PetscReal *h = box->h; 235 const PetscInt *n = box->n; 236 const PetscInt dim = box->dim; 237 PetscInt d, p; 238 239 PetscFunctionBegin; 240 for (p = 0; p < numPoints; ++p) { 241 for (d = 0; d < dim; ++d) { 242 PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]); 243 244 if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1; 245 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", 246 p, PetscRealPart(points[p*dim+0]), dim > 1 ? PetscRealPart(points[p*dim+1]) : 0.0, dim > 2 ? PetscRealPart(points[p*dim+2]) : 0.0); 247 dboxes[p*dim+d] = dbox; 248 } 249 if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1]; 250 } 251 PetscFunctionReturn(0); 252 } 253 254 #undef __FUNCT__ 255 #define __FUNCT__ "PetscGridHashDestroy" 256 PetscErrorCode PetscGridHashDestroy(PetscGridHash *box) 257 { 258 PetscErrorCode ierr; 259 260 PetscFunctionBegin; 261 if (*box) { 262 ierr = PetscSectionDestroy(&(*box)->cellSection);CHKERRQ(ierr); 263 ierr = ISDestroy(&(*box)->cells);CHKERRQ(ierr); 264 ierr = DMLabelDestroy(&(*box)->cellsSparse);CHKERRQ(ierr); 265 } 266 ierr = PetscFree(*box);CHKERRQ(ierr); 267 PetscFunctionReturn(0); 268 } 269 270 #undef __FUNCT__ 271 #define __FUNCT__ "DMPlexLocatePoint_Internal" 272 PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell) 273 { 274 PetscInt coneSize; 275 PetscErrorCode ierr; 276 277 PetscFunctionBegin; 278 switch (dim) { 279 case 2: 280 ierr = DMPlexGetConeSize(dm, cellStart, &coneSize);CHKERRQ(ierr); 281 switch (coneSize) { 282 case 3: 283 ierr = DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);CHKERRQ(ierr); 284 break; 285 case 4: 286 ierr = DMPlexLocatePoint_General_2D_Internal(dm, point, cellStart, cell);CHKERRQ(ierr); 287 break; 288 default: 289 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize); 290 } 291 break; 292 case 3: 293 ierr = DMPlexGetConeSize(dm, cellStart, &coneSize);CHKERRQ(ierr); 294 switch (coneSize) { 295 case 4: 296 ierr = DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);CHKERRQ(ierr); 297 break; 298 case 6: 299 ierr = DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);CHKERRQ(ierr); 300 break; 301 default: 302 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize); 303 } 304 break; 305 default: 306 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for mesh dimension %D", dim); 307 } 308 PetscFunctionReturn(0); 309 } 310 311 #undef __FUNCT__ 312 #define __FUNCT__ "DMPlexComputeGridHash_Internal" 313 PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox) 314 { 315 MPI_Comm comm; 316 PetscGridHash lbox; 317 Vec coordinates; 318 PetscSection coordSection; 319 Vec coordsLocal; 320 const PetscScalar *coords; 321 PetscInt *dboxes, *boxes; 322 PetscInt n[3] = {10, 10, 10}; 323 PetscInt dim, N, cStart, cEnd, cMax, c, i; 324 PetscErrorCode ierr; 325 326 PetscFunctionBegin; 327 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 328 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 329 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 330 ierr = VecGetLocalSize(coordinates, &N);CHKERRQ(ierr); 331 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 332 ierr = PetscGridHashCreate(comm, dim, coords, &lbox);CHKERRQ(ierr); 333 for (i = 0; i < N; i += dim) {ierr = PetscGridHashEnlarge(lbox, &coords[i]);CHKERRQ(ierr);} 334 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 335 ierr = PetscGridHashSetGrid(lbox, n, NULL);CHKERRQ(ierr); 336 #if 0 337 /* Could define a custom reduction to merge these */ 338 ierr = MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);CHKERRQ(ierr); 339 ierr = MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);CHKERRQ(ierr); 340 #endif 341 /* Is there a reason to snap the local bounding box to a division of the global box? */ 342 /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */ 343 /* Create label */ 344 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 345 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 346 if (cMax >= 0) cEnd = PetscMin(cEnd, cMax); 347 ierr = DMLabelCreate("cells", &lbox->cellsSparse);CHKERRQ(ierr); 348 ierr = DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);CHKERRQ(ierr); 349 /* Compute boxes which overlap each cell: http://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */ 350 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 351 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 352 ierr = PetscCalloc2(16 * dim, &dboxes, 16, &boxes);CHKERRQ(ierr); 353 for (c = cStart; c < cEnd; ++c) { 354 const PetscReal *h = lbox->h; 355 PetscScalar *ccoords = NULL; 356 PetscInt csize = 0; 357 PetscScalar point[3]; 358 PetscInt dlim[6], d, e, i, j, k; 359 360 /* Find boxes enclosing each vertex */ 361 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);CHKERRQ(ierr); 362 ierr = PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);CHKERRQ(ierr); 363 /* Mark cells containing the vertices */ 364 for (e = 0; e < csize/dim; ++e) {ierr = DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);CHKERRQ(ierr);} 365 /* Get grid of boxes containing these */ 366 for (d = 0; d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];} 367 for (d = dim; d < 3; ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;} 368 for (e = 1; e < dim+1; ++e) { 369 for (d = 0; d < dim; ++d) { 370 dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]); 371 dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]); 372 } 373 } 374 /* Check for intersection of box with cell */ 375 for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) { 376 for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) { 377 for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) { 378 const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i; 379 PetscScalar cpoint[3]; 380 PetscInt cell, edge, ii, jj, kk; 381 382 /* Check whether cell contains any vertex of these subboxes TODO vectorize this */ 383 for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) { 384 for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) { 385 for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) { 386 387 ierr = DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);CHKERRQ(ierr); 388 if (cell >= 0) {DMLabelSetValue(lbox->cellsSparse, c, box);CHKERRQ(ierr); ii = jj = kk = 2;} 389 } 390 } 391 } 392 /* Check whether cell edge intersects any edge of these subboxes TODO vectorize this */ 393 for (edge = 0; edge < dim+1; ++edge) { 394 PetscReal segA[6], segB[6]; 395 396 for (d = 0; d < dim; ++d) {segA[d] = PetscRealPart(ccoords[edge*dim+d]); segA[dim+d] = PetscRealPart(ccoords[((edge+1)%(dim+1))*dim+d]);} 397 for (kk = 0; kk < (dim > 2 ? 2 : 1); ++kk) { 398 if (dim > 2) {segB[2] = PetscRealPart(point[2]); 399 segB[dim+2] = PetscRealPart(point[2]) + kk*h[2];} 400 for (jj = 0; jj < (dim > 1 ? 2 : 1); ++jj) { 401 if (dim > 1) {segB[1] = PetscRealPart(point[1]); 402 segB[dim+1] = PetscRealPart(point[1]) + jj*h[1];} 403 for (ii = 0; ii < 2; ++ii) { 404 PetscBool intersects; 405 406 segB[0] = PetscRealPart(point[0]); 407 segB[dim+0] = PetscRealPart(point[0]) + ii*h[0]; 408 ierr = DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);CHKERRQ(ierr); 409 if (intersects) {DMLabelSetValue(lbox->cellsSparse, c, box);CHKERRQ(ierr); edge = ii = jj = kk = dim+1;} 410 } 411 } 412 } 413 } 414 } 415 } 416 } 417 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);CHKERRQ(ierr); 418 } 419 ierr = PetscFree2(dboxes, boxes);CHKERRQ(ierr); 420 ierr = DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);CHKERRQ(ierr); 421 ierr = DMLabelDestroy(&lbox->cellsSparse);CHKERRQ(ierr); 422 *localBox = lbox; 423 PetscFunctionReturn(0); 424 } 425 426 #undef __FUNCT__ 427 #define __FUNCT__ "DMLocatePoints_Plex" 428 PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, PetscSF cellSF) 429 { 430 DM_Plex *mesh = (DM_Plex *) dm->data; 431 PetscBool hash = mesh->useHashLocation; 432 PetscInt bs, numPoints, p, numFound, *found = NULL; 433 PetscInt dim, cStart, cEnd, cMax, numCells, c; 434 const PetscInt *boxCells; 435 PetscSFNode *cells; 436 PetscScalar *a; 437 PetscMPIInt result; 438 PetscErrorCode ierr; 439 440 PetscFunctionBegin; 441 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 442 ierr = VecGetBlockSize(v, &bs);CHKERRQ(ierr); 443 ierr = MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);CHKERRQ(ierr); 444 if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported"); 445 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); 446 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 447 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 448 if (cMax >= 0) cEnd = PetscMin(cEnd, cMax); 449 ierr = VecGetLocalSize(v, &numPoints);CHKERRQ(ierr); 450 ierr = VecGetArray(v, &a);CHKERRQ(ierr); 451 numPoints /= bs; 452 ierr = PetscMalloc1(numPoints, &cells);CHKERRQ(ierr); 453 if (hash) { 454 if (!mesh->lbox) {ierr = PetscInfo(dm, "Initializing grid hashing");CHKERRQ(ierr);ierr = DMPlexComputeGridHash_Internal(dm, &mesh->lbox);CHKERRQ(ierr);} 455 /* Designate the local box for each point */ 456 /* Send points to correct process */ 457 /* Search cells that lie in each subbox */ 458 /* Should we bin points before doing search? */ 459 ierr = ISGetIndices(mesh->lbox->cells, &boxCells);CHKERRQ(ierr); 460 } 461 for (p = 0, numFound = 0; p < numPoints; ++p) { 462 const PetscScalar *point = &a[p*bs]; 463 PetscInt dbin[3], bin, cell = -1, cellOffset; 464 465 cells[p].rank = -1; 466 cells[p].index = -1; 467 if (hash) { 468 ierr = PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);CHKERRQ(ierr); 469 /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */ 470 ierr = PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);CHKERRQ(ierr); 471 ierr = PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);CHKERRQ(ierr); 472 for (c = cellOffset; c < cellOffset + numCells; ++c) { 473 ierr = DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell);CHKERRQ(ierr); 474 if (cell >= 0) { 475 cells[p].rank = 0; 476 cells[p].index = cell; 477 numFound++; 478 break; 479 } 480 } 481 } else { 482 for (c = cStart; c < cEnd; ++c) { 483 ierr = DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);CHKERRQ(ierr); 484 if (cell >= 0) { 485 cells[p].rank = 0; 486 cells[p].index = cell; 487 numFound++; 488 break; 489 } 490 } 491 } 492 } 493 if (hash) {ierr = ISRestoreIndices(mesh->lbox->cells, &boxCells);CHKERRQ(ierr);} 494 /* Check for highest numbered proc that claims a point (do we care?) */ 495 ierr = VecRestoreArray(v, &a);CHKERRQ(ierr); 496 if (numFound < numPoints) { 497 ierr = PetscMalloc1(numFound,&found);CHKERRQ(ierr); 498 for (p = 0, numFound = 0; p < numPoints; p++) { 499 if (cells[p].rank >= 0 && cells[p].index >= 0) { 500 if (numFound < p) { 501 cells[numFound] = cells[p]; 502 } 503 found[numFound++] = p; 504 } 505 } 506 } 507 ierr = PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);CHKERRQ(ierr); 508 PetscFunctionReturn(0); 509 } 510 511 #undef __FUNCT__ 512 #define __FUNCT__ "DMPlexComputeProjection2Dto1D_Internal" 513 /* 514 DMPlexComputeProjection2Dto1D_Internal - Rewrite coordinates to be the 1D projection of the 2D 515 */ 516 PetscErrorCode DMPlexComputeProjection2Dto1D_Internal(PetscScalar coords[], PetscReal R[]) 517 { 518 const PetscReal x = PetscRealPart(coords[2] - coords[0]); 519 const PetscReal y = PetscRealPart(coords[3] - coords[1]); 520 const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r; 521 522 PetscFunctionBegin; 523 R[0] = c; R[1] = -s; 524 R[2] = s; R[3] = c; 525 coords[0] = 0.0; 526 coords[1] = r; 527 PetscFunctionReturn(0); 528 } 529 530 #undef __FUNCT__ 531 #define __FUNCT__ "DMPlexComputeProjection3Dto1D_Internal" 532 /* 533 DMPlexComputeProjection3Dto1D_Internal - Rewrite coordinates to be the 1D projection of the 3D 534 535 This uses the basis completion described by Frisvad, 536 537 http://www.imm.dtu.dk/~jerf/papers/abstracts/onb.html 538 DOI:10.1080/2165347X.2012.689606 539 */ 540 PetscErrorCode DMPlexComputeProjection3Dto1D_Internal(PetscScalar coords[], PetscReal R[]) 541 { 542 PetscReal x = PetscRealPart(coords[3] - coords[0]); 543 PetscReal y = PetscRealPart(coords[4] - coords[1]); 544 PetscReal z = PetscRealPart(coords[5] - coords[2]); 545 PetscReal r = PetscSqrtReal(x*x + y*y + z*z); 546 PetscReal rinv = 1. / r; 547 PetscFunctionBegin; 548 549 x *= rinv; y *= rinv; z *= rinv; 550 if (x > 0.) { 551 PetscReal inv1pX = 1./ (1. + x); 552 553 R[0] = x; R[1] = -y; R[2] = -z; 554 R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] = -y*z*inv1pX; 555 R[6] = z; R[7] = -y*z*inv1pX; R[8] = 1. - z*z*inv1pX; 556 } 557 else { 558 PetscReal inv1mX = 1./ (1. - x); 559 560 R[0] = x; R[1] = z; R[2] = y; 561 R[3] = y; R[4] = -y*z*inv1mX; R[5] = 1. - y*y*inv1mX; 562 R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] = -y*z*inv1mX; 563 } 564 coords[0] = 0.0; 565 coords[1] = r; 566 PetscFunctionReturn(0); 567 } 568 569 #undef __FUNCT__ 570 #define __FUNCT__ "DMPlexComputeProjection3Dto2D_Internal" 571 /* 572 DMPlexComputeProjection3Dto2D_Internal - Rewrite coordinates to be the 2D projection of the 3D 573 */ 574 PetscErrorCode DMPlexComputeProjection3Dto2D_Internal(PetscInt coordSize, PetscScalar coords[], PetscReal R[]) 575 { 576 PetscReal x1[3], x2[3], n[3], norm; 577 PetscReal x1p[3], x2p[3], xnp[3]; 578 PetscReal sqrtz, alpha; 579 const PetscInt dim = 3; 580 PetscInt d, e, p; 581 582 PetscFunctionBegin; 583 /* 0) Calculate normal vector */ 584 for (d = 0; d < dim; ++d) { 585 x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]); 586 x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]); 587 } 588 n[0] = x1[1]*x2[2] - x1[2]*x2[1]; 589 n[1] = x1[2]*x2[0] - x1[0]*x2[2]; 590 n[2] = x1[0]*x2[1] - x1[1]*x2[0]; 591 norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]); 592 n[0] /= norm; 593 n[1] /= norm; 594 n[2] /= norm; 595 /* 1) Take the normal vector and rotate until it is \hat z 596 597 Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then 598 599 R = / alpha nx nz alpha ny nz -1/alpha \ 600 | -alpha ny alpha nx 0 | 601 \ nx ny nz / 602 603 will rotate the normal vector to \hat z 604 */ 605 sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]); 606 /* Check for n = z */ 607 if (sqrtz < 1.0e-10) { 608 const PetscInt s = PetscSign(n[2]); 609 /* If nz < 0, rotate 180 degrees around x-axis */ 610 for (p = 3; p < coordSize/3; ++p) { 611 coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]); 612 coords[p*2+1] = (PetscRealPart(coords[p*dim+1] - coords[0*dim+1])) * s; 613 } 614 coords[0] = 0.0; 615 coords[1] = 0.0; 616 coords[2] = x1[0]; 617 coords[3] = x1[1] * s; 618 coords[4] = x2[0]; 619 coords[5] = x2[1] * s; 620 R[0] = 1.0; R[1] = 0.0; R[2] = 0.0; 621 R[3] = 0.0; R[4] = 1.0 * s; R[5] = 0.0; 622 R[6] = 0.0; R[7] = 0.0; R[8] = 1.0 * s; 623 PetscFunctionReturn(0); 624 } 625 alpha = 1.0/sqrtz; 626 R[0] = alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz; 627 R[3] = -alpha*n[1]; R[4] = alpha*n[0]; R[5] = 0.0; 628 R[6] = n[0]; R[7] = n[1]; R[8] = n[2]; 629 for (d = 0; d < dim; ++d) { 630 x1p[d] = 0.0; 631 x2p[d] = 0.0; 632 for (e = 0; e < dim; ++e) { 633 x1p[d] += R[d*dim+e]*x1[e]; 634 x2p[d] += R[d*dim+e]*x2[e]; 635 } 636 } 637 if (PetscAbsReal(x1p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 638 if (PetscAbsReal(x2p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated"); 639 /* 2) Project to (x, y) */ 640 for (p = 3; p < coordSize/3; ++p) { 641 for (d = 0; d < dim; ++d) { 642 xnp[d] = 0.0; 643 for (e = 0; e < dim; ++e) { 644 xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]); 645 } 646 if (d < dim-1) coords[p*2+d] = xnp[d]; 647 } 648 } 649 coords[0] = 0.0; 650 coords[1] = 0.0; 651 coords[2] = x1p[0]; 652 coords[3] = x1p[1]; 653 coords[4] = x2p[0]; 654 coords[5] = x2p[1]; 655 /* Output R^T which rotates \hat z to the input normal */ 656 for (d = 0; d < dim; ++d) { 657 for (e = d+1; e < dim; ++e) { 658 PetscReal tmp; 659 660 tmp = R[d*dim+e]; 661 R[d*dim+e] = R[e*dim+d]; 662 R[e*dim+d] = tmp; 663 } 664 } 665 PetscFunctionReturn(0); 666 } 667 668 #undef __FUNCT__ 669 #define __FUNCT__ "Volume_Triangle_Internal" 670 PETSC_UNUSED 671 PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[]) 672 { 673 /* Signed volume is 1/2 the determinant 674 675 | 1 1 1 | 676 | x0 x1 x2 | 677 | y0 y1 y2 | 678 679 but if x0,y0 is the origin, we have 680 681 | x1 x2 | 682 | y1 y2 | 683 */ 684 const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1]; 685 const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1]; 686 PetscReal M[4], detM; 687 M[0] = x1; M[1] = x2; 688 M[2] = y1; M[3] = y2; 689 DMPlex_Det2D_Internal(&detM, M); 690 *vol = 0.5*detM; 691 (void)PetscLogFlops(5.0); 692 } 693 694 #undef __FUNCT__ 695 #define __FUNCT__ "Volume_Triangle_Origin_Internal" 696 PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[]) 697 { 698 DMPlex_Det2D_Internal(vol, coords); 699 *vol *= 0.5; 700 } 701 702 #undef __FUNCT__ 703 #define __FUNCT__ "Volume_Tetrahedron_Internal" 704 PETSC_UNUSED 705 PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[]) 706 { 707 /* Signed volume is 1/6th of the determinant 708 709 | 1 1 1 1 | 710 | x0 x1 x2 x3 | 711 | y0 y1 y2 y3 | 712 | z0 z1 z2 z3 | 713 714 but if x0,y0,z0 is the origin, we have 715 716 | x1 x2 x3 | 717 | y1 y2 y3 | 718 | z1 z2 z3 | 719 */ 720 const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2]; 721 const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2]; 722 const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2]; 723 PetscReal M[9], detM; 724 M[0] = x1; M[1] = x2; M[2] = x3; 725 M[3] = y1; M[4] = y2; M[5] = y3; 726 M[6] = z1; M[7] = z2; M[8] = z3; 727 DMPlex_Det3D_Internal(&detM, M); 728 *vol = -0.16666666666666666666666*detM; 729 (void)PetscLogFlops(10.0); 730 } 731 732 #undef __FUNCT__ 733 #define __FUNCT__ "Volume_Tetrahedron_Origin_Internal" 734 PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[]) 735 { 736 DMPlex_Det3D_Internal(vol, coords); 737 *vol *= -0.16666666666666666666666; 738 } 739 740 #undef __FUNCT__ 741 #define __FUNCT__ "DMPlexComputeLineGeometry_Internal" 742 static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 743 { 744 PetscSection coordSection; 745 Vec coordinates; 746 PetscScalar *coords = NULL; 747 PetscInt numCoords, d, pStart, pEnd, numSelfCoords = 0; 748 PetscErrorCode ierr; 749 750 PetscFunctionBegin; 751 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 752 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 753 ierr = PetscSectionGetChart(coordSection,&pStart,&pEnd);CHKERRQ(ierr); 754 if (e >= pStart && e < pEnd) {ierr = PetscSectionGetDof(coordSection,e,&numSelfCoords);CHKERRQ(ierr);} 755 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 756 numCoords = numSelfCoords ? numSelfCoords : numCoords; 757 if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL"); 758 *detJ = 0.0; 759 if (numCoords == 6) { 760 const PetscInt dim = 3; 761 PetscReal R[9], J0; 762 763 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 764 ierr = DMPlexComputeProjection3Dto1D_Internal(coords, R);CHKERRQ(ierr); 765 if (J) { 766 J0 = 0.5*PetscRealPart(coords[1]); 767 J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2]; 768 J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5]; 769 J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8]; 770 DMPlex_Det3D_Internal(detJ, J); 771 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 772 } 773 } else if (numCoords == 4) { 774 const PetscInt dim = 2; 775 PetscReal R[4], J0; 776 777 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 778 ierr = DMPlexComputeProjection2Dto1D_Internal(coords, R);CHKERRQ(ierr); 779 if (J) { 780 J0 = 0.5*PetscRealPart(coords[1]); 781 J[0] = R[0]*J0; J[1] = R[1]; 782 J[2] = R[2]*J0; J[3] = R[3]; 783 DMPlex_Det2D_Internal(detJ, J); 784 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 785 } 786 } else if (numCoords == 2) { 787 const PetscInt dim = 1; 788 789 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 790 if (J) { 791 J[0] = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0])); 792 *detJ = J[0]; 793 ierr = PetscLogFlops(2.0);CHKERRQ(ierr); 794 if (invJ) {invJ[0] = 1.0/J[0]; ierr = PetscLogFlops(1.0);CHKERRQ(ierr);} 795 } 796 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords); 797 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 798 PetscFunctionReturn(0); 799 } 800 801 #undef __FUNCT__ 802 #define __FUNCT__ "DMPlexComputeTriangleGeometry_Internal" 803 static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 804 { 805 PetscSection coordSection; 806 Vec coordinates; 807 PetscScalar *coords = NULL; 808 PetscInt numCoords, d, f, g; 809 PetscErrorCode ierr; 810 811 PetscFunctionBegin; 812 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 813 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 814 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 815 *detJ = 0.0; 816 if (numCoords == 9) { 817 const PetscInt dim = 3; 818 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 819 820 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 821 ierr = DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);CHKERRQ(ierr); 822 if (J) { 823 const PetscInt pdim = 2; 824 825 for (d = 0; d < pdim; d++) { 826 for (f = 0; f < pdim; f++) { 827 J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 828 } 829 } 830 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 831 DMPlex_Det3D_Internal(detJ, J0); 832 for (d = 0; d < dim; d++) { 833 for (f = 0; f < dim; f++) { 834 J[d*dim+f] = 0.0; 835 for (g = 0; g < dim; g++) { 836 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 837 } 838 } 839 } 840 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 841 } 842 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 843 } else if (numCoords == 6) { 844 const PetscInt dim = 2; 845 846 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 847 if (J) { 848 for (d = 0; d < dim; d++) { 849 for (f = 0; f < dim; f++) { 850 J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d])); 851 } 852 } 853 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 854 DMPlex_Det2D_Internal(detJ, J); 855 } 856 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 857 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords); 858 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 859 PetscFunctionReturn(0); 860 } 861 862 #undef __FUNCT__ 863 #define __FUNCT__ "DMPlexComputeRectangleGeometry_Internal" 864 static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 865 { 866 PetscSection coordSection; 867 Vec coordinates; 868 PetscScalar *coords = NULL; 869 PetscInt numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd; 870 PetscErrorCode ierr; 871 872 PetscFunctionBegin; 873 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 874 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 875 ierr = PetscSectionGetChart(coordSection,&pStart,&pEnd);CHKERRQ(ierr); 876 if (e >= pStart && e < pEnd) {ierr = PetscSectionGetDof(coordSection,e,&numSelfCoords);CHKERRQ(ierr);} 877 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 878 numCoords = numSelfCoords ? numSelfCoords : numCoords; 879 *detJ = 0.0; 880 if (numCoords == 12) { 881 const PetscInt dim = 3; 882 PetscReal R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0}; 883 884 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 885 ierr = DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);CHKERRQ(ierr); 886 if (J) { 887 const PetscInt pdim = 2; 888 889 for (d = 0; d < pdim; d++) { 890 J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 891 J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d])); 892 } 893 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 894 DMPlex_Det3D_Internal(detJ, J0); 895 for (d = 0; d < dim; d++) { 896 for (f = 0; f < dim; f++) { 897 J[d*dim+f] = 0.0; 898 for (g = 0; g < dim; g++) { 899 J[d*dim+f] += R[d*dim+g]*J0[g*dim+f]; 900 } 901 } 902 } 903 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 904 } 905 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 906 } else if (numCoords == 8) { 907 const PetscInt dim = 2; 908 909 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 910 if (J) { 911 for (d = 0; d < dim; d++) { 912 J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 913 J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 914 } 915 ierr = PetscLogFlops(8.0);CHKERRQ(ierr); 916 DMPlex_Det2D_Internal(detJ, J); 917 } 918 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 919 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords); 920 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr); 921 PetscFunctionReturn(0); 922 } 923 924 #undef __FUNCT__ 925 #define __FUNCT__ "DMPlexComputeTetrahedronGeometry_Internal" 926 static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 927 { 928 PetscSection coordSection; 929 Vec coordinates; 930 PetscScalar *coords = NULL; 931 const PetscInt dim = 3; 932 PetscInt d; 933 PetscErrorCode ierr; 934 935 PetscFunctionBegin; 936 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 937 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 938 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 939 *detJ = 0.0; 940 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 941 if (J) { 942 for (d = 0; d < dim; d++) { 943 /* I orient with outward face normals */ 944 J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d])); 945 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 946 J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 947 } 948 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 949 DMPlex_Det3D_Internal(detJ, J); 950 } 951 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 952 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 953 PetscFunctionReturn(0); 954 } 955 956 #undef __FUNCT__ 957 #define __FUNCT__ "DMPlexComputeHexahedronGeometry_Internal" 958 static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 959 { 960 PetscSection coordSection; 961 Vec coordinates; 962 PetscScalar *coords = NULL; 963 const PetscInt dim = 3; 964 PetscInt d; 965 PetscErrorCode ierr; 966 967 PetscFunctionBegin; 968 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 969 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 970 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 971 *detJ = 0.0; 972 if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);} 973 if (J) { 974 for (d = 0; d < dim; d++) { 975 J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d])); 976 J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d])); 977 J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d])); 978 } 979 ierr = PetscLogFlops(18.0);CHKERRQ(ierr); 980 DMPlex_Det3D_Internal(detJ, J); 981 } 982 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 983 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr); 984 PetscFunctionReturn(0); 985 } 986 987 #undef __FUNCT__ 988 #define __FUNCT__ "DMPlexComputeCellGeometryAffineFEM" 989 /*@C 990 DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell 991 992 Collective on DM 993 994 Input Arguments: 995 + dm - the DM 996 - cell - the cell 997 998 Output Arguments: 999 + v0 - the translation part of this affine transform 1000 . J - the Jacobian of the transform from the reference element 1001 . invJ - the inverse of the Jacobian 1002 - detJ - the Jacobian determinant 1003 1004 Level: advanced 1005 1006 Fortran Notes: 1007 Since it returns arrays, this routine is only available in Fortran 90, and you must 1008 include petsc.h90 in your code. 1009 1010 .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinateVec() 1011 @*/ 1012 PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1013 { 1014 PetscInt depth, dim, coneSize; 1015 PetscErrorCode ierr; 1016 1017 PetscFunctionBegin; 1018 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1019 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 1020 if (depth == 1) { 1021 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1022 } else { 1023 DMLabel depth; 1024 1025 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 1026 ierr = DMLabelGetValue(depth, cell, &dim);CHKERRQ(ierr); 1027 } 1028 switch (dim) { 1029 case 1: 1030 ierr = DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1031 break; 1032 case 2: 1033 switch (coneSize) { 1034 case 3: 1035 ierr = DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1036 break; 1037 case 4: 1038 ierr = DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1039 break; 1040 default: 1041 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 1042 } 1043 break; 1044 case 3: 1045 switch (coneSize) { 1046 case 4: 1047 ierr = DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1048 break; 1049 case 6: /* Faces */ 1050 case 8: /* Vertices */ 1051 ierr = DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr); 1052 break; 1053 default: 1054 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell); 1055 } 1056 break; 1057 default: 1058 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1059 } 1060 PetscFunctionReturn(0); 1061 } 1062 1063 #undef __FUNCT__ 1064 #define __FUNCT__ "DMPlexComputeIsoparametricGeometry_Internal" 1065 static PetscErrorCode DMPlexComputeIsoparametricGeometry_Internal(DM dm, PetscFE fe, PetscInt point, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 1066 { 1067 PetscQuadrature quad; 1068 PetscSection coordSection; 1069 Vec coordinates; 1070 PetscScalar *coords = NULL; 1071 const PetscReal *quadPoints; 1072 PetscReal *basisDer; 1073 PetscInt dim, cdim, pdim, qdim, Nq, numCoords, d, q; 1074 PetscErrorCode ierr; 1075 1076 PetscFunctionBegin; 1077 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1078 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1079 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1080 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1081 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 1082 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1083 ierr = PetscFEGetDimension(fe, &pdim);CHKERRQ(ierr); 1084 ierr = PetscQuadratureGetData(quad, &qdim, &Nq, &quadPoints, NULL);CHKERRQ(ierr); 1085 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 1086 *detJ = 0.0; 1087 if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim); 1088 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); 1089 if (v0) {for (d = 0; d < cdim; d++) v0[d] = PetscRealPart(coords[d]);} 1090 if (J) { 1091 ierr = PetscMemzero(J, Nq*cdim*dim*sizeof(PetscReal));CHKERRQ(ierr); 1092 for (q = 0; q < Nq; ++q) { 1093 PetscInt i, j, k, c, r; 1094 1095 /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */ 1096 for (k = 0; k < pdim; ++k) 1097 for (j = 0; j < dim; ++j) 1098 for (i = 0; i < cdim; ++i) 1099 J[(q*cdim + i)*dim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]); 1100 ierr = PetscLogFlops(2.0*pdim*dim*cdim);CHKERRQ(ierr); 1101 if (cdim > dim) { 1102 for (c = dim; c < cdim; ++c) 1103 for (r = 0; r < cdim; ++r) 1104 J[r*cdim+c] = r == c ? 1.0 : 0.0; 1105 } 1106 switch (cdim) { 1107 case 3: 1108 DMPlex_Det3D_Internal(detJ, J); 1109 if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);} 1110 break; 1111 case 2: 1112 DMPlex_Det2D_Internal(detJ, J); 1113 if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);} 1114 break; 1115 case 1: 1116 *detJ = J[0]; 1117 if (invJ) invJ[0] = 1.0/J[0]; 1118 } 1119 } 1120 } 1121 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);CHKERRQ(ierr); 1122 PetscFunctionReturn(0); 1123 } 1124 1125 #undef __FUNCT__ 1126 #define __FUNCT__ "DMPlexComputeCellGeometryFEM" 1127 /*@C 1128 DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell 1129 1130 Collective on DM 1131 1132 Input Arguments: 1133 + dm - the DM 1134 . cell - the cell 1135 - fe - the finite element containing the quadrature 1136 1137 Output Arguments: 1138 + v0 - the translation part of this transform 1139 . J - the Jacobian of the transform from the reference element at each quadrature point 1140 . invJ - the inverse of the Jacobian at each quadrature point 1141 - detJ - the Jacobian determinant at each quadrature point 1142 1143 Level: advanced 1144 1145 Fortran Notes: 1146 Since it returns arrays, this routine is only available in Fortran 90, and you must 1147 include petsc.h90 in your code. 1148 1149 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1150 @*/ 1151 PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscFE fe, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) 1152 { 1153 PetscErrorCode ierr; 1154 1155 PetscFunctionBegin; 1156 if (!fe) {ierr = DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, detJ);CHKERRQ(ierr);} 1157 else {ierr = DMPlexComputeIsoparametricGeometry_Internal(dm, fe, cell, v0, J, invJ, detJ);CHKERRQ(ierr);} 1158 PetscFunctionReturn(0); 1159 } 1160 1161 #undef __FUNCT__ 1162 #define __FUNCT__ "DMPlexComputeGeometryFVM_1D_Internal" 1163 static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1164 { 1165 PetscSection coordSection; 1166 Vec coordinates; 1167 PetscScalar *coords = NULL; 1168 PetscScalar tmp[2]; 1169 PetscInt coordSize; 1170 PetscErrorCode ierr; 1171 1172 PetscFunctionBegin; 1173 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1174 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1175 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1176 if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now"); 1177 ierr = DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);CHKERRQ(ierr); 1178 if (centroid) { 1179 centroid[0] = 0.5*PetscRealPart(coords[0] + tmp[0]); 1180 centroid[1] = 0.5*PetscRealPart(coords[1] + tmp[1]); 1181 } 1182 if (normal) { 1183 PetscReal norm; 1184 1185 normal[0] = -PetscRealPart(coords[1] - tmp[1]); 1186 normal[1] = PetscRealPart(coords[0] - tmp[0]); 1187 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]); 1188 normal[0] /= norm; 1189 normal[1] /= norm; 1190 } 1191 if (vol) { 1192 *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - tmp[0])) + PetscSqr(PetscRealPart(coords[1] - tmp[1]))); 1193 } 1194 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1195 PetscFunctionReturn(0); 1196 } 1197 1198 #undef __FUNCT__ 1199 #define __FUNCT__ "DMPlexComputeGeometryFVM_2D_Internal" 1200 /* Centroid_i = (\sum_n A_n Cn_i ) / A */ 1201 static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1202 { 1203 PetscSection coordSection; 1204 Vec coordinates; 1205 PetscScalar *coords = NULL; 1206 PetscReal vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9]; 1207 PetscInt tdim = 2, coordSize, numCorners, p, d, e; 1208 PetscErrorCode ierr; 1209 1210 PetscFunctionBegin; 1211 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1212 ierr = DMPlexGetConeSize(dm, cell, &numCorners);CHKERRQ(ierr); 1213 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1214 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1215 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 1216 if (dim > 2 && centroid) { 1217 v0[0] = PetscRealPart(coords[0]); 1218 v0[1] = PetscRealPart(coords[1]); 1219 v0[2] = PetscRealPart(coords[2]); 1220 } 1221 if (normal) { 1222 if (dim > 2) { 1223 const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]); 1224 const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]); 1225 const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]); 1226 PetscReal norm; 1227 1228 normal[0] = y0*z1 - z0*y1; 1229 normal[1] = z0*x1 - x0*z1; 1230 normal[2] = x0*y1 - y0*x1; 1231 norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]); 1232 normal[0] /= norm; 1233 normal[1] /= norm; 1234 normal[2] /= norm; 1235 } else { 1236 for (d = 0; d < dim; ++d) normal[d] = 0.0; 1237 } 1238 } 1239 if (dim == 3) {ierr = DMPlexComputeProjection3Dto2D_Internal(coordSize, coords, R);CHKERRQ(ierr);} 1240 for (p = 0; p < numCorners; ++p) { 1241 /* Need to do this copy to get types right */ 1242 for (d = 0; d < tdim; ++d) { 1243 ctmp[d] = PetscRealPart(coords[p*tdim+d]); 1244 ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]); 1245 } 1246 Volume_Triangle_Origin_Internal(&vtmp, ctmp); 1247 vsum += vtmp; 1248 for (d = 0; d < tdim; ++d) { 1249 csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp; 1250 } 1251 } 1252 for (d = 0; d < tdim; ++d) { 1253 csum[d] /= (tdim+1)*vsum; 1254 } 1255 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr); 1256 if (vol) *vol = PetscAbsReal(vsum); 1257 if (centroid) { 1258 if (dim > 2) { 1259 for (d = 0; d < dim; ++d) { 1260 centroid[d] = v0[d]; 1261 for (e = 0; e < dim; ++e) { 1262 centroid[d] += R[d*dim+e]*csum[e]; 1263 } 1264 } 1265 } else for (d = 0; d < dim; ++d) centroid[d] = csum[d]; 1266 } 1267 PetscFunctionReturn(0); 1268 } 1269 1270 #undef __FUNCT__ 1271 #define __FUNCT__ "DMPlexComputeGeometryFVM_3D_Internal" 1272 /* Centroid_i = (\sum_n V_n Cn_i ) / V */ 1273 static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1274 { 1275 PetscSection coordSection; 1276 Vec coordinates; 1277 PetscScalar *coords = NULL; 1278 PetscReal vsum = 0.0, vtmp, coordsTmp[3*3]; 1279 const PetscInt *faces, *facesO; 1280 PetscInt numFaces, f, coordSize, numCorners, p, d; 1281 PetscErrorCode ierr; 1282 1283 PetscFunctionBegin; 1284 if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim); 1285 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1286 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1287 1288 if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0; 1289 ierr = DMPlexGetConeSize(dm, cell, &numFaces);CHKERRQ(ierr); 1290 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 1291 ierr = DMPlexGetConeOrientation(dm, cell, &facesO);CHKERRQ(ierr); 1292 for (f = 0; f < numFaces; ++f) { 1293 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 1294 numCorners = coordSize/dim; 1295 switch (numCorners) { 1296 case 3: 1297 for (d = 0; d < dim; ++d) { 1298 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 1299 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 1300 coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]); 1301 } 1302 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1303 if (facesO[f] < 0) vtmp = -vtmp; 1304 vsum += vtmp; 1305 if (centroid) { /* Centroid of OABC = (a+b+c)/4 */ 1306 for (d = 0; d < dim; ++d) { 1307 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1308 } 1309 } 1310 break; 1311 case 4: 1312 /* DO FOR PYRAMID */ 1313 /* First tet */ 1314 for (d = 0; d < dim; ++d) { 1315 coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]); 1316 coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]); 1317 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 1318 } 1319 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1320 if (facesO[f] < 0) vtmp = -vtmp; 1321 vsum += vtmp; 1322 if (centroid) { 1323 for (d = 0; d < dim; ++d) { 1324 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1325 } 1326 } 1327 /* Second tet */ 1328 for (d = 0; d < dim; ++d) { 1329 coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]); 1330 coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]); 1331 coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]); 1332 } 1333 Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp); 1334 if (facesO[f] < 0) vtmp = -vtmp; 1335 vsum += vtmp; 1336 if (centroid) { 1337 for (d = 0; d < dim; ++d) { 1338 for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp; 1339 } 1340 } 1341 break; 1342 default: 1343 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners); 1344 } 1345 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);CHKERRQ(ierr); 1346 } 1347 if (vol) *vol = PetscAbsReal(vsum); 1348 if (normal) for (d = 0; d < dim; ++d) normal[d] = 0.0; 1349 if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4); 1350 PetscFunctionReturn(0); 1351 } 1352 1353 #undef __FUNCT__ 1354 #define __FUNCT__ "DMPlexComputeCellGeometryFVM" 1355 /*@C 1356 DMPlexComputeCellGeometryFVM - Compute the volume for a given cell 1357 1358 Collective on DM 1359 1360 Input Arguments: 1361 + dm - the DM 1362 - cell - the cell 1363 1364 Output Arguments: 1365 + volume - the cell volume 1366 . centroid - the cell centroid 1367 - normal - the cell normal, if appropriate 1368 1369 Level: advanced 1370 1371 Fortran Notes: 1372 Since it returns arrays, this routine is only available in Fortran 90, and you must 1373 include petsc.h90 in your code. 1374 1375 .seealso: DMGetCoordinateSection(), DMGetCoordinateVec() 1376 @*/ 1377 PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[]) 1378 { 1379 PetscInt depth, dim; 1380 PetscErrorCode ierr; 1381 1382 PetscFunctionBegin; 1383 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1384 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1385 if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated"); 1386 /* We need to keep a pointer to the depth label */ 1387 ierr = DMGetLabelValue(dm, "depth", cell, &depth);CHKERRQ(ierr); 1388 /* Cone size is now the number of faces */ 1389 switch (depth) { 1390 case 1: 1391 ierr = DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1392 break; 1393 case 2: 1394 ierr = DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1395 break; 1396 case 3: 1397 ierr = DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);CHKERRQ(ierr); 1398 break; 1399 default: 1400 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim); 1401 } 1402 PetscFunctionReturn(0); 1403 } 1404 1405 #undef __FUNCT__ 1406 #define __FUNCT__ "DMPlexComputeGeometryFEM" 1407 /* This should also take a PetscFE argument I think */ 1408 PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom) 1409 { 1410 DM dmCell; 1411 Vec coordinates; 1412 PetscSection coordSection, sectionCell; 1413 PetscScalar *cgeom; 1414 PetscInt cStart, cEnd, cMax, c; 1415 PetscErrorCode ierr; 1416 1417 PetscFunctionBegin; 1418 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1419 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1420 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1421 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1422 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1423 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1424 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1425 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1426 cEnd = cMax < 0 ? cEnd : cMax; 1427 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1428 /* TODO This needs to be multiplied by Nq for non-affine */ 1429 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFECellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1430 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1431 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1432 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1433 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1434 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1435 for (c = cStart; c < cEnd; ++c) { 1436 PetscFECellGeom *cg; 1437 1438 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1439 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1440 ierr = DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v0, cg->J, cg->invJ, &cg->detJ);CHKERRQ(ierr); 1441 if (cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c); 1442 } 1443 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1444 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1445 PetscFunctionReturn(0); 1446 } 1447 1448 #undef __FUNCT__ 1449 #define __FUNCT__ "DMPlexComputeGeometryFVM" 1450 /*@ 1451 DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method 1452 1453 Input Parameter: 1454 . dm - The DM 1455 1456 Output Parameters: 1457 + cellgeom - A Vec of PetscFVCellGeom data 1458 . facegeom - A Vec of PetscFVFaceGeom data 1459 1460 Level: developer 1461 1462 .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM() 1463 @*/ 1464 PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom) 1465 { 1466 DM dmFace, dmCell; 1467 DMLabel ghostLabel; 1468 PetscSection sectionFace, sectionCell; 1469 PetscSection coordSection; 1470 Vec coordinates; 1471 PetscScalar *fgeom, *cgeom; 1472 PetscReal minradius, gminradius; 1473 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 1474 PetscErrorCode ierr; 1475 1476 PetscFunctionBegin; 1477 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1478 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1479 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1480 /* Make cell centroids and volumes */ 1481 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 1482 ierr = DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);CHKERRQ(ierr); 1483 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 1484 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 1485 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1486 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1487 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 1488 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1489 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 1490 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 1491 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 1492 ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr); 1493 if (cEndInterior < 0) { 1494 cEndInterior = cEnd; 1495 } 1496 ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1497 for (c = cStart; c < cEndInterior; ++c) { 1498 PetscFVCellGeom *cg; 1499 1500 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1501 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 1502 ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr); 1503 } 1504 /* Compute face normals and minimum cell radius */ 1505 ierr = DMClone(dm, &dmFace);CHKERRQ(ierr); 1506 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);CHKERRQ(ierr); 1507 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1508 ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr); 1509 for (f = fStart; f < fEnd; ++f) {ierr = PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));CHKERRQ(ierr);} 1510 ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr); 1511 ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr); 1512 ierr = PetscSectionDestroy(§ionFace);CHKERRQ(ierr); 1513 ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr); 1514 ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr); 1515 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1516 minradius = PETSC_MAX_REAL; 1517 for (f = fStart; f < fEnd; ++f) { 1518 PetscFVFaceGeom *fg; 1519 PetscReal area; 1520 PetscInt ghost = -1, d, numChildren; 1521 1522 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1523 ierr = DMPlexGetTreeChildren(dm,f,&numChildren,NULL);CHKERRQ(ierr); 1524 if (ghost >= 0 || numChildren) continue; 1525 ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr); 1526 ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);CHKERRQ(ierr); 1527 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 1528 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 1529 { 1530 PetscFVCellGeom *cL, *cR; 1531 PetscInt ncells; 1532 const PetscInt *cells; 1533 PetscReal *lcentroid, *rcentroid; 1534 PetscReal l[3], r[3], v[3]; 1535 1536 ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr); 1537 ierr = DMPlexGetSupportSize(dm, f, &ncells);CHKERRQ(ierr); 1538 ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr); 1539 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 1540 if (ncells > 1) { 1541 ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr); 1542 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 1543 } 1544 else { 1545 rcentroid = fg->centroid; 1546 } 1547 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);CHKERRQ(ierr); 1548 ierr = DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);CHKERRQ(ierr); 1549 DMPlex_WaxpyD_Internal(dim, -1, l, r, v); 1550 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) { 1551 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 1552 } 1553 if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) { 1554 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]); 1555 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]); 1556 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f); 1557 } 1558 if (cells[0] < cEndInterior) { 1559 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v); 1560 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1561 } 1562 if (ncells > 1 && cells[1] < cEndInterior) { 1563 DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v); 1564 minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v)); 1565 } 1566 } 1567 } 1568 ierr = MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1569 ierr = DMPlexSetMinRadius(dm, gminradius);CHKERRQ(ierr); 1570 /* Compute centroids of ghost cells */ 1571 for (c = cEndInterior; c < cEnd; ++c) { 1572 PetscFVFaceGeom *fg; 1573 const PetscInt *cone, *support; 1574 PetscInt coneSize, supportSize, s; 1575 1576 ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr); 1577 if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize); 1578 ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr); 1579 ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr); 1580 if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize); 1581 ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr); 1582 ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr); 1583 for (s = 0; s < 2; ++s) { 1584 /* Reflect ghost centroid across plane of face */ 1585 if (support[s] == c) { 1586 PetscFVCellGeom *ci; 1587 PetscFVCellGeom *cg; 1588 PetscReal c2f[3], a; 1589 1590 ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr); 1591 DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 1592 a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal); 1593 ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr); 1594 DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid); 1595 cg->volume = ci->volume; 1596 } 1597 } 1598 } 1599 ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr); 1600 ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr); 1601 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 1602 ierr = DMDestroy(&dmFace);CHKERRQ(ierr); 1603 PetscFunctionReturn(0); 1604 } 1605 1606 #undef __FUNCT__ 1607 #define __FUNCT__ "DMPlexGetMinRadius" 1608 /*@C 1609 DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face 1610 1611 Not collective 1612 1613 Input Argument: 1614 . dm - the DM 1615 1616 Output Argument: 1617 . minradius - the minium cell radius 1618 1619 Level: developer 1620 1621 .seealso: DMGetCoordinates() 1622 @*/ 1623 PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius) 1624 { 1625 PetscFunctionBegin; 1626 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1627 PetscValidPointer(minradius,2); 1628 *minradius = ((DM_Plex*) dm->data)->minradius; 1629 PetscFunctionReturn(0); 1630 } 1631 1632 #undef __FUNCT__ 1633 #define __FUNCT__ "DMPlexSetMinRadius" 1634 /*@C 1635 DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face 1636 1637 Logically collective 1638 1639 Input Arguments: 1640 + dm - the DM 1641 - minradius - the minium cell radius 1642 1643 Level: developer 1644 1645 .seealso: DMSetCoordinates() 1646 @*/ 1647 PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius) 1648 { 1649 PetscFunctionBegin; 1650 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1651 ((DM_Plex*) dm->data)->minradius = minradius; 1652 PetscFunctionReturn(0); 1653 } 1654 1655 #undef __FUNCT__ 1656 #define __FUNCT__ "BuildGradientReconstruction_Internal" 1657 static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 1658 { 1659 DMLabel ghostLabel; 1660 PetscScalar *dx, *grad, **gref; 1661 PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces; 1662 PetscErrorCode ierr; 1663 1664 PetscFunctionBegin; 1665 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1666 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1667 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1668 ierr = DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);CHKERRQ(ierr); 1669 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 1670 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1671 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 1672 for (c = cStart; c < cEndInterior; c++) { 1673 const PetscInt *faces; 1674 PetscInt numFaces, usedFaces, f, d; 1675 PetscFVCellGeom *cg; 1676 PetscBool boundary; 1677 PetscInt ghost; 1678 1679 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1680 ierr = DMPlexGetConeSize(dm, c, &numFaces);CHKERRQ(ierr); 1681 ierr = DMPlexGetCone(dm, c, &faces);CHKERRQ(ierr); 1682 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces); 1683 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 1684 PetscFVCellGeom *cg1; 1685 PetscFVFaceGeom *fg; 1686 const PetscInt *fcells; 1687 PetscInt ncell, side; 1688 1689 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 1690 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 1691 if ((ghost >= 0) || boundary) continue; 1692 ierr = DMPlexGetSupport(dm, faces[f], &fcells);CHKERRQ(ierr); 1693 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 1694 ncell = fcells[!side]; /* the neighbor */ 1695 ierr = DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);CHKERRQ(ierr); 1696 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 1697 for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d]; 1698 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 1699 } 1700 if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?"); 1701 ierr = PetscFVComputeGradient(fvm, usedFaces, dx, grad);CHKERRQ(ierr); 1702 for (f = 0, usedFaces = 0; f < numFaces; ++f) { 1703 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 1704 ierr = DMIsBoundaryPoint(dm, faces[f], &boundary);CHKERRQ(ierr); 1705 if ((ghost >= 0) || boundary) continue; 1706 for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d]; 1707 ++usedFaces; 1708 } 1709 } 1710 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 1711 PetscFunctionReturn(0); 1712 } 1713 1714 #undef __FUNCT__ 1715 #define __FUNCT__ "BuildGradientReconstruction_Internal_Tree" 1716 static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom) 1717 { 1718 DMLabel ghostLabel; 1719 PetscScalar *dx, *grad, **gref; 1720 PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0; 1721 PetscSection neighSec; 1722 PetscInt (*neighbors)[2]; 1723 PetscInt *counter; 1724 PetscErrorCode ierr; 1725 1726 PetscFunctionBegin; 1727 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1728 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1729 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1730 if (cEndInterior < 0) { 1731 cEndInterior = cEnd; 1732 } 1733 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);CHKERRQ(ierr); 1734 ierr = PetscSectionSetChart(neighSec,cStart,cEndInterior);CHKERRQ(ierr); 1735 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1736 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 1737 for (f = fStart; f < fEnd; f++) { 1738 const PetscInt *fcells; 1739 PetscBool boundary; 1740 PetscInt ghost = -1; 1741 PetscInt numChildren, numCells, c; 1742 1743 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1744 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 1745 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 1746 if ((ghost >= 0) || boundary || numChildren) continue; 1747 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 1748 if (numCells == 2) { 1749 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 1750 for (c = 0; c < 2; c++) { 1751 PetscInt cell = fcells[c]; 1752 1753 if (cell >= cStart && cell < cEndInterior) { 1754 ierr = PetscSectionAddDof(neighSec,cell,1);CHKERRQ(ierr); 1755 } 1756 } 1757 } 1758 } 1759 ierr = PetscSectionSetUp(neighSec);CHKERRQ(ierr); 1760 ierr = PetscSectionGetMaxDof(neighSec,&maxNumFaces);CHKERRQ(ierr); 1761 ierr = PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);CHKERRQ(ierr); 1762 nStart = 0; 1763 ierr = PetscSectionGetStorageSize(neighSec,&nEnd);CHKERRQ(ierr); 1764 ierr = PetscMalloc1((nEnd-nStart),&neighbors);CHKERRQ(ierr); 1765 ierr = PetscCalloc1((cEndInterior-cStart),&counter);CHKERRQ(ierr); 1766 for (f = fStart; f < fEnd; f++) { 1767 const PetscInt *fcells; 1768 PetscBool boundary; 1769 PetscInt ghost = -1; 1770 PetscInt numChildren, numCells, c; 1771 1772 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);} 1773 ierr = DMIsBoundaryPoint(dm, f, &boundary);CHKERRQ(ierr); 1774 ierr = DMPlexGetTreeChildren(dm, f, &numChildren, NULL);CHKERRQ(ierr); 1775 if ((ghost >= 0) || boundary || numChildren) continue; 1776 ierr = DMPlexGetSupportSize(dm, f, &numCells);CHKERRQ(ierr); 1777 if (numCells == 2) { 1778 ierr = DMPlexGetSupport(dm, f, &fcells);CHKERRQ(ierr); 1779 for (c = 0; c < 2; c++) { 1780 PetscInt cell = fcells[c], off; 1781 1782 if (cell >= cStart && cell < cEndInterior) { 1783 ierr = PetscSectionGetOffset(neighSec,cell,&off);CHKERRQ(ierr); 1784 off += counter[cell - cStart]++; 1785 neighbors[off][0] = f; 1786 neighbors[off][1] = fcells[1 - c]; 1787 } 1788 } 1789 } 1790 } 1791 ierr = PetscFree(counter);CHKERRQ(ierr); 1792 ierr = PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);CHKERRQ(ierr); 1793 for (c = cStart; c < cEndInterior; c++) { 1794 PetscInt numFaces, f, d, off, ghost = -1; 1795 PetscFVCellGeom *cg; 1796 1797 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 1798 ierr = PetscSectionGetDof(neighSec, c, &numFaces);CHKERRQ(ierr); 1799 ierr = PetscSectionGetOffset(neighSec, c, &off);CHKERRQ(ierr); 1800 if (ghostLabel) {ierr = DMLabelGetValue(ghostLabel, c, &ghost);CHKERRQ(ierr);} 1801 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); 1802 for (f = 0; f < numFaces; ++f) { 1803 PetscFVCellGeom *cg1; 1804 PetscFVFaceGeom *fg; 1805 const PetscInt *fcells; 1806 PetscInt ncell, side, nface; 1807 1808 nface = neighbors[off + f][0]; 1809 ncell = neighbors[off + f][1]; 1810 ierr = DMPlexGetSupport(dm,nface,&fcells);CHKERRQ(ierr); 1811 side = (c != fcells[0]); 1812 ierr = DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);CHKERRQ(ierr); 1813 ierr = DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);CHKERRQ(ierr); 1814 for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d]; 1815 gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */ 1816 } 1817 ierr = PetscFVComputeGradient(fvm, numFaces, dx, grad);CHKERRQ(ierr); 1818 for (f = 0; f < numFaces; ++f) { 1819 for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d]; 1820 } 1821 } 1822 ierr = PetscFree3(dx, grad, gref);CHKERRQ(ierr); 1823 ierr = PetscSectionDestroy(&neighSec);CHKERRQ(ierr); 1824 ierr = PetscFree(neighbors);CHKERRQ(ierr); 1825 PetscFunctionReturn(0); 1826 } 1827 1828 #undef __FUNCT__ 1829 #define __FUNCT__ "DMPlexComputeGradientFVM" 1830 /*@ 1831 DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data 1832 1833 Collective on DM 1834 1835 Input Arguments: 1836 + dm - The DM 1837 . fvm - The PetscFV 1838 . faceGeometry - The face geometry from DMPlexGetFaceGeometryFVM() 1839 - cellGeometry - The face geometry from DMPlexGetCellGeometryFVM() 1840 1841 Output Parameters: 1842 + faceGeometry - The geometric factors for gradient calculation are inserted 1843 - dmGrad - The DM describing the layout of gradient data 1844 1845 Level: developer 1846 1847 .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM() 1848 @*/ 1849 PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad) 1850 { 1851 DM dmFace, dmCell; 1852 PetscScalar *fgeom, *cgeom; 1853 PetscSection sectionGrad, parentSection; 1854 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 1855 PetscErrorCode ierr; 1856 1857 PetscFunctionBegin; 1858 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1859 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 1860 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1861 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1862 /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */ 1863 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 1864 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 1865 ierr = VecGetArray(faceGeometry, &fgeom);CHKERRQ(ierr); 1866 ierr = VecGetArray(cellGeometry, &cgeom);CHKERRQ(ierr); 1867 ierr = DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 1868 if (!parentSection) { 1869 ierr = BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 1870 } 1871 else { 1872 ierr = BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 1873 } 1874 ierr = VecRestoreArray(faceGeometry, &fgeom);CHKERRQ(ierr); 1875 ierr = VecRestoreArray(cellGeometry, &cgeom);CHKERRQ(ierr); 1876 /* Create storage for gradients */ 1877 ierr = DMClone(dm, dmGrad);CHKERRQ(ierr); 1878 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 1879 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 1880 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 1881 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 1882 ierr = DMSetDefaultSection(*dmGrad, sectionGrad);CHKERRQ(ierr); 1883 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 1884 PetscFunctionReturn(0); 1885 } 1886 1887 #undef __FUNCT__ 1888 #define __FUNCT__ "DMPlexGetDataFVM" 1889 PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM) 1890 { 1891 PetscObject cellgeomobj, facegeomobj; 1892 PetscErrorCode ierr; 1893 1894 PetscFunctionBegin; 1895 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 1896 if (!cellgeomobj) { 1897 Vec cellgeomInt, facegeomInt; 1898 1899 ierr = DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);CHKERRQ(ierr); 1900 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);CHKERRQ(ierr); 1901 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);CHKERRQ(ierr); 1902 ierr = VecDestroy(&cellgeomInt);CHKERRQ(ierr); 1903 ierr = VecDestroy(&facegeomInt);CHKERRQ(ierr); 1904 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);CHKERRQ(ierr); 1905 } 1906 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);CHKERRQ(ierr); 1907 if (cellgeom) *cellgeom = (Vec) cellgeomobj; 1908 if (facegeom) *facegeom = (Vec) facegeomobj; 1909 if (gradDM) { 1910 PetscObject gradobj; 1911 PetscBool computeGradients; 1912 1913 ierr = PetscFVGetComputeGradients(fv,&computeGradients);CHKERRQ(ierr); 1914 if (!computeGradients) { 1915 *gradDM = NULL; 1916 PetscFunctionReturn(0); 1917 } 1918 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 1919 if (!gradobj) { 1920 DM dmGradInt; 1921 1922 ierr = DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);CHKERRQ(ierr); 1923 ierr = PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);CHKERRQ(ierr); 1924 ierr = DMDestroy(&dmGradInt);CHKERRQ(ierr); 1925 ierr = PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);CHKERRQ(ierr); 1926 } 1927 *gradDM = (DM) gradobj; 1928 } 1929 PetscFunctionReturn(0); 1930 } 1931