Skip to content

Commit

Permalink
Updated some petsc support functions to work with ghost cells (#518)
Browse files Browse the repository at this point in the history 
* Updated some petsc support functions to work with ghost cells

* Updated formatting

* Version bump
  • Loading branch information
@dsalac
dsalac authored Feb 23, 2024
1 parent 14f6541 commit 33e8421
Show file tree
Hide file tree
Showing 2 changed files with 38 additions and 243 deletions.
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ cmake_minimum_required(VERSION 3.18.4)
include(config/petscCompilers.cmake)

# Set the project details
project(ablateLibrary VERSION 0.12.23)
project(ablateLibrary VERSION 0.12.24)

# Load the Required 3rd Party Libaries
pkg_check_modules(PETSc REQUIRED IMPORTED_TARGET GLOBAL PETSc)
Expand Down
279 changes: 37 additions & 242 deletions src/utilities/petscSupport.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include "petscSupport.hpp"
#include <petsc/private/vecimpl.h>
#include <petscdm.h> // For DMPolytopeTypeGetNumVertices

/**
* Return the cell containing the location xyz
Expand Down Expand Up @@ -258,7 +259,7 @@ PetscErrorCode DMPlexRestoreNeighbors(DM dm, PetscInt p, PetscInt maxLevels, Pet

PetscErrorCode DMPlexGetNeighbors(DM dm, PetscInt p, PetscInt maxLevels, PetscReal maxDist, PetscInt numberCells, PetscBool useCells, PetscBool returnNeighborVertices, PetscInt *nCells,
PetscInt **cells) {
const PetscInt maxLevelListSize = 10000;
const PetscInt maxLevelListSize = 100000;
const PetscInt maxListSize = 100000;
PetscInt numNew, nLevelList[2];
PetscInt *addList = NULL, levelList[2][maxLevelListSize], currentLevelLoc, prevLevelLoc;
Expand Down Expand Up @@ -381,7 +382,7 @@ PetscErrorCode DMPlexGetNeighbors(DM dm, PetscInt p, PetscInt maxLevels, PetscRe
PetscCall(PetscSortRealWithArrayInt(n, dist, list));
PetscCall(PetscFree(dist));
} else if (type == 0 && cte == 1) {
// Now only include the the numberCells closest cells
// Now only include the the numberCells closest vertices
PetscReal *dist;
PetscInt j, dim, i_x, vStart;
Vec coords;
Expand Down Expand Up @@ -419,36 +420,7 @@ PetscErrorCode DMPlexCellGetNumVertices(DM dm, const PetscInt p, PetscInt *nv) {

PetscCall(DMPlexGetCellType(dm, p, &ct));

switch (ct) {
case DM_POLYTOPE_POINT:
*nv = 1;
break;
case DM_POLYTOPE_SEGMENT:
case DM_POLYTOPE_POINT_PRISM_TENSOR:
*nv = 2;
break;
case DM_POLYTOPE_TRIANGLE:
*nv = 3;
break;
case DM_POLYTOPE_QUADRILATERAL:
case DM_POLYTOPE_SEG_PRISM_TENSOR:
case DM_POLYTOPE_TETRAHEDRON:
*nv = 4;
break;
case DM_POLYTOPE_PYRAMID:
*nv = 5;
break;
case DM_POLYTOPE_TRI_PRISM:
case DM_POLYTOPE_TRI_PRISM_TENSOR:
*nv = 6;
break;
case DM_POLYTOPE_HEXAHEDRON:
case DM_POLYTOPE_QUAD_PRISM_TENSOR:
*nv = 8;
break;
default:
SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cannot determine number of vertices for cell type %s", DMPolytopeTypes[ct]);
}
*nv = DMPolytopeTypeGetNumVertices(ct);

PetscFunctionReturn(PETSC_SUCCESS);
}
Expand Down Expand Up @@ -601,7 +573,7 @@ PetscErrorCode xDMPlexPointLocalRead(DM dm, PetscInt p, PetscInt fID, const Pets
* @param centroid - Centroid of the face
* @param n - Outward facing surface area normal
*/
static PetscErrorCode DMPlexFaceCentroidOutwardAreaNormal(DM dm, PetscInt cell, PetscInt face, PetscReal *centroid, PetscReal *n) {
PetscErrorCode DMPlexFaceCentroidOutwardAreaNormal(DM dm, PetscInt cell, PetscInt face, PetscReal *centroid, PetscReal *n) {
PetscFunctionBegin;

// Get the cell center
Expand Down Expand Up @@ -892,6 +864,8 @@ PetscErrorCode DMPlexRestoreCommonPoints(DM dm, const PetscInt p1, const PetscIn
PetscFunctionReturn(PETSC_SUCCESS);
}

// Compute the corner surface area normal as defined in Morgan and Waltz with respect to a given vertex and an edge center
// NOTE: This does NOT check if the vertex and cell are actually associated with each other.
PetscErrorCode DMPlexCornerSurfaceAreaNormal(DM dm, const PetscInt v, const PetscInt c, PetscReal N[]) {
PetscFunctionBegin;

Expand Down Expand Up @@ -1156,160 +1130,7 @@ PetscErrorCode DMPlexCellGradFromVertex(DM dm, const PetscInt c, Vec data, Petsc
PetscFunctionReturn(PETSC_SUCCESS);
}

// Copied from petsc/src/dm/impls/plex/plexgeometry
static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection) {
const PetscReal p0_x = segmentA[0 * 2 + 0];
const PetscReal p0_y = segmentA[0 * 2 + 1];
const PetscReal p1_x = segmentA[1 * 2 + 0];
const PetscReal p1_y = segmentA[1 * 2 + 1];
const PetscReal p2_x = segmentB[0 * 2 + 0];
const PetscReal p2_y = segmentB[0 * 2 + 1];
const PetscReal p3_x = segmentB[1 * 2 + 0];
const PetscReal p3_y = segmentB[1 * 2 + 1];
const PetscReal s1_x = p1_x - p0_x;
const PetscReal s1_y = p1_y - p0_y;
const PetscReal s2_x = p3_x - p2_x;
const PetscReal s2_y = p3_y - p2_y;
const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);

PetscFunctionBegin;
*hasIntersection = PETSC_FALSE;
/* Non-parallel lines */
if (denom != 0.0) {
const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
const PetscReal t = (s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;

if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
*hasIntersection = PETSC_TRUE;
if (intersection) {
intersection[0] = p0_x + (t * s1_x);
intersection[1] = p0_y + (t * s1_y);
}
}
}
PetscFunctionReturn(PETSC_SUCCESS);
}

/* The plane is segmentB x segmentC: https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection */
static PetscErrorCode DMPlexGetLineTriangleIntersection_3D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], const PetscReal segmentC[], PetscReal intersection[],
PetscBool *hasIntersection) {
const PetscReal p0_x = segmentA[0 * 3 + 0];
const PetscReal p0_y = segmentA[0 * 3 + 1];
const PetscReal p0_z = segmentA[0 * 3 + 2];
const PetscReal p1_x = segmentA[1 * 3 + 0];
const PetscReal p1_y = segmentA[1 * 3 + 1];
const PetscReal p1_z = segmentA[1 * 3 + 2];

const PetscReal q0_x = segmentB[0 * 3 + 0];
const PetscReal q0_y = segmentB[0 * 3 + 1];
const PetscReal q0_z = segmentB[0 * 3 + 2];
const PetscReal q1_x = segmentB[1 * 3 + 0];
const PetscReal q1_y = segmentB[1 * 3 + 1];
const PetscReal q1_z = segmentB[1 * 3 + 2];

const PetscReal r0_x = segmentC[0 * 3 + 0];
const PetscReal r0_y = segmentC[0 * 3 + 1];
const PetscReal r0_z = segmentC[0 * 3 + 2];
const PetscReal r1_x = segmentC[1 * 3 + 0];
const PetscReal r1_y = segmentC[1 * 3 + 1];
const PetscReal r1_z = segmentC[1 * 3 + 2];

const PetscReal s0_x = p1_x - p0_x;
const PetscReal s0_y = p1_y - p0_y;
const PetscReal s0_z = p1_z - p0_z;
const PetscReal s1_x = q1_x - q0_x;
const PetscReal s1_y = q1_y - q0_y;
const PetscReal s1_z = q1_z - q0_z;
const PetscReal s2_x = r1_x - r0_x;
const PetscReal s2_y = r1_y - r0_y;
const PetscReal s2_z = r1_z - r0_z;
const PetscReal s3_x = s1_y * s2_z - s1_z * s2_y; /* s1 x s2 */
const PetscReal s3_y = s1_z * s2_x - s1_x * s2_z;
const PetscReal s3_z = s1_x * s2_y - s1_y * s2_x;
const PetscReal s4_x = s0_y * s2_z - s0_z * s2_y; /* s0 x s2 */
const PetscReal s4_y = s0_z * s2_x - s0_x * s2_z;
const PetscReal s4_z = s0_x * s2_y - s0_y * s2_x;
const PetscReal s5_x = s1_y * s0_z - s1_z * s0_y; /* s1 x s0 */
const PetscReal s5_y = s1_z * s0_x - s1_x * s0_z;
const PetscReal s5_z = s1_x * s0_y - s1_y * s0_x;
const PetscReal denom = -(s0_x * s3_x + s0_y * s3_y + s0_z * s3_z); /* -s0 . (s1 x s2) */

PetscFunctionBegin;
*hasIntersection = PETSC_FALSE;
/* Line not parallel to plane */
if (denom != 0.0) {
const PetscReal t = (s3_x * (p0_x - q0_x) + s3_y * (p0_y - q0_y) + s3_z * (p0_z - q0_z)) / denom;
const PetscReal u = (s4_x * (p0_x - q0_x) + s4_y * (p0_y - q0_y) + s4_z * (p0_z - q0_z)) / denom;
const PetscReal v = (s5_x * (p0_x - q0_x) + s5_y * (p0_y - q0_y) + s5_z * (p0_z - q0_z)) / denom;

if (t >= 0 && t <= 1 && u >= 0 && u <= 1 && v >= 0 && v <= 1 && (u + v) <= 1) {
*hasIntersection = PETSC_TRUE;
if (intersection) {
intersection[0] = p0_x + (t * s0_x);
intersection[1] = p0_y + (t * s0_y);
intersection[2] = p0_z + (t * s0_z);
}
}
// printf("%+f\t%+f\t%+f\n", t, u, v);
}
PetscFunctionReturn(PETSC_SUCCESS);
}

// Determine the intersection of a face and a line
PetscErrorCode DMPlexFaceLineIntersection(DM dm, const PetscInt f, const PetscReal segment[], PetscReal intersection[], PetscBool *hasIntersection) {
PetscInt dim;
const PetscScalar *array;
PetscScalar *coords = NULL;
PetscInt numCoords;
PetscBool isDG;

PetscFunctionBegin;
PetscCall(DMGetDimension(dm, &dim));

PetscCall(DMPlexGetCellCoordinates(dm, f, &isDG, &numCoords, &array, &coords));

switch (dim) {
case 1: {
*hasIntersection = PETSC_FALSE;
} break;
case 2: {
PetscCall(DMPlexGetLineIntersection_2D_Internal(coords, segment, intersection, hasIntersection));
} break;
case 3: {
PetscReal segmentA[6], segmentB[6];

// Triangle 0 - 1 - 2
segmentA[0] = segmentB[0] = coords[3];
segmentA[1] = segmentB[1] = coords[4];
segmentA[2] = segmentB[2] = coords[5];

segmentA[3] = coords[0];
segmentA[4] = coords[1];
segmentA[5] = coords[2];

segmentB[3] = coords[6];
segmentB[4] = coords[7];
segmentB[5] = coords[8];

PetscCall(DMPlexGetLineTriangleIntersection_3D_Internal(segment, segmentA, segmentB, intersection, hasIntersection));

if (numCoords == 12 && !(*hasIntersection)) { // The face is a quad, so try triangle 0 - 3 - 2
segmentA[0] = segmentB[0] = coords[9];
segmentA[1] = segmentB[1] = coords[10];
segmentA[2] = segmentB[2] = coords[11];
PetscCall(DMPlexGetLineTriangleIntersection_3D_Internal(segment, segmentA, segmentB, intersection, hasIntersection));
}
} break;
default:
PetscFunctionReturn(PETSC_ERR_SUP);
}

PetscCall(DMPlexRestoreCellCoordinates(dm, f, &isDG, &numCoords, &array, &coords));

PetscFunctionReturn(PETSC_SUCCESS);
}

#include <petsc/private/hashmapi.h>
// This isn't the most accurate as the face center may not lie along the vector connecting two adjacent cells
PetscErrorCode DMPlexCellGradFromCell(DM dm, const PetscInt c, Vec data, PetscInt fID, PetscInt offset, PetscScalar g[]) {
PetscFunctionBegin;

Expand All @@ -1324,73 +1145,47 @@ PetscErrorCode DMPlexCellGradFromCell(DM dm, const PetscInt c, Vec data, PetscIn
const PetscScalar *dataArray;
PetscCall(VecGetArrayRead(data, &dataArray));

// Get all vertices of the cell
PetscInt nVert, *verts;
PetscCall(DMPlexCellGetVertices(dm, c, &nVert, &verts));

PetscHMapI hash = NULL; // Used to convert from vertex numbering to 0->nVert-1
PetscCall(PetscHMapICreate(&hash));
// Get all faces of the cell
PetscInt nFaces;
const PetscInt *faces;
PetscCall(DMPlexGetConeSize(dm, c, &nFaces));
PetscCall(DMPlexGetCone(dm, c, &faces));

PetscReal *vertVals;
PetscCall(DMGetWorkArray(dm, nVert, MPIU_REAL, &vertVals));
for (PetscInt d = 0; d < dim; ++d) g[d] = 0.0;

// Locations of the vertices
PetscReal *vertCoords;
PetscCall(DMPlexVertexGetCoordinates(dm, nVert, verts, &vertCoords));
for (PetscInt v = 0; v < nVert; ++v) {
PetscCall(PetscHMapISet(hash, verts[v], v));
for (PetscInt f = 0; f < nFaces; ++f) {
// Compute the face center location and the outward surface area normal
PetscReal S[dim], centroid[dim];
PetscCall(DMPlexFaceCentroidOutwardAreaNormal(dm, c, faces[f], centroid, S));

vertVals[v] = 0.0;
// The cells sharing this face
PetscInt nSharedCells;
const PetscInt *sharedCells;
PetscCall(DMPlexGetSupportSize(dm, faces[f], &nSharedCells));
PetscCall(DMPlexGetSupport(dm, faces[f], &sharedCells));
PetscCheck(nSharedCells < 3, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "More than two cells are sharing a face.");

PetscInt nCells, *cells;
PetscCall(DMPlexVertexGetCells(dm, verts[v], &nCells, &cells));
PetscReal vAve = 0.0, vSum = 0.0;
for (PetscInt j = 0; j < nSharedCells; ++j) {
PetscInt sc = sharedCells[j];

PetscReal totalWt = 0.0;
for (PetscInt i = 0; i < nCells; ++i) {
PetscReal *cellVal, cellCenter[3], wt = 0.0;
PetscCall(xDMPlexPointLocalRead(dm, cells[i], fID, dataArray, &cellVal));
PetscCall(DMPlexComputeCellGeometryFVM(dm, cells[i], NULL, cellCenter, NULL));
for (PetscInt d = 0; d < dim; ++d) wt += PetscSqr(cellCenter[d] - vertCoords[v * dim + d]);
wt = 1.0 / PetscSqrtReal(wt);
totalWt += wt;
vertVals[v] += wt * cellVal[offset];
}
PetscCall(DMPlexVertexRestoreCells(dm, verts[v], &nCells, &cells));
PetscReal x[dim];
PetscCall(DMPlexComputeCellGeometryFVM(dm, sc, NULL, x, NULL)); // Center of the candidate cell.

vertVals[v] /= totalWt;
}
PetscCall(DMPlexVertexRestoreCoordinates(dm, nVert, verts, &vertCoords));
PetscCall(DMPlexCellRestoreVertices(dm, c, &nVert, &verts));
PetscReal dist = 0.0;
for (PetscInt d = 0; d < dim; ++d) dist += PetscSqr(x[d] - centroid[d]);
dist = PetscSqrtReal(dist);

for (PetscInt d = 0; d < dim; ++d) g[d] = 0.0;
PetscReal *val;
PetscCall(xDMPlexPointLocalRead(dm, sc, fID, dataArray, &val));

PetscInt nFaces;
const PetscInt *faces;
PetscCall(DMPlexGetConeSize(dm, c, &nFaces));
PetscCall(DMPlexGetCone(dm, c, &faces));
for (PetscInt f = 0; f < nFaces; ++f) {
PetscInt nVert, *verts;
PetscCall(DMPlexCellGetVertices(dm, faces[f], &nVert, &verts));

PetscReal faceValue = 0.0;
for (PetscInt v = 0; v < nVert; ++v) {
PetscInt id;
PetscHMapIGet(hash, verts[v], &id);
faceValue += vertVals[id];
vAve += val[offset] * dist;
vSum += dist;
}
faceValue /= nVert;

PetscCall(DMPlexCellRestoreVertices(dm, faces[f], &nVert, &verts));

PetscReal N[dim];
PetscCall(DMPlexFaceCentroidOutwardAreaNormal(dm, c, faces[f], NULL, N));

for (PetscInt d = 0; d < dim; ++d) g[d] += faceValue * N[d];
for (PetscInt d = 0; d < dim; ++d) g[d] += vAve * S[d] / vSum;
}

PetscHMapIDestroy(&hash);
PetscCall(DMRestoreWorkArray(dm, nVert, MPIU_REAL, &vertVals));

// Center of the cell
PetscReal cellVolume;
PetscCall(DMPlexComputeCellGeometryFVM(dm, c, &cellVolume, NULL, NULL));
Expand Down

0 comments on commit 33e8421

Please sign in to comment.