|
| 1 | +diff --git a/src/dm/impls/plex/plexfem.c b/src/dm/impls/plex/plexfem.c |
| 2 | +index 0299f569070..1f07f49f60f 100644 |
| 3 | +--- a/src/dm/impls/plex/plexfem.c |
| 4 | ++++ b/src/dm/impls/plex/plexfem.c |
| 5 | +@@ -2874,6 +2874,24 @@ PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt num |
| 6 | + PetscCall(DMPlexGetDepthLabel(dm, &depthLabel)); |
| 7 | + PetscCall(DMGetDimension(dm, &dim)); |
| 8 | + PetscCall(DMLabelGetStratumIS(depthLabel, dim - 1, &facetIS)); |
| 9 | ++ /* Filter out ghost facets (SF leaves) so each facet contributes once */ |
| 10 | ++ if (facetIS) { |
| 11 | ++ PetscSF sf; |
| 12 | ++ PetscInt nleaves; |
| 13 | ++ const PetscInt *leaves; |
| 14 | ++ |
| 15 | ++ PetscCall(DMGetPointSF(dm, &sf)); |
| 16 | ++ PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL)); |
| 17 | ++ if (nleaves > 0 && leaves) { |
| 18 | ++ IS leafIS, ownedFacetIS; |
| 19 | ++ |
| 20 | ++ PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nleaves, leaves, PETSC_USE_POINTER, &leafIS)); |
| 21 | ++ PetscCall(ISDifference(facetIS, leafIS, &ownedFacetIS)); |
| 22 | ++ PetscCall(ISDestroy(&leafIS)); |
| 23 | ++ PetscCall(ISDestroy(&facetIS)); |
| 24 | ++ facetIS = ownedFacetIS; |
| 25 | ++ } |
| 26 | ++ } |
| 27 | + PetscCall(DMGetLocalSection(dm, §ion)); |
| 28 | + PetscCall(PetscSectionGetNumFields(section, &Nf)); |
| 29 | + /* Get local solution with boundary values */ |
| 30 | +@@ -4961,6 +4979,28 @@ PetscErrorCode DMPlexComputeBdResidualSingleByKey(DM dm, PetscWeakForm wf, Petsc |
| 31 | + PetscCall(ISDestroy(&pointIS)); |
| 32 | + pointIS = isectIS; |
| 33 | + } |
| 34 | ++ if (key.part > 0) { |
| 35 | ++ IS filteredIS = NULL; |
| 36 | ++ const PetscInt *allPoints; |
| 37 | ++ PetscInt *keptPoints = NULL; |
| 38 | ++ PetscInt nAllFaces, nKept = 0, fidx; |
| 39 | ++ |
| 40 | ++ PetscCall(ISGetLocalSize(pointIS, &nAllFaces)); |
| 41 | ++ PetscCall(ISGetIndices(pointIS, &allPoints)); |
| 42 | ++ PetscCall(PetscMalloc1(nAllFaces, &keptPoints)); |
| 43 | ++ for (fidx = 0; fidx < nAllFaces; ++fidx) { |
| 44 | ++ const PetscInt point = allPoints[fidx]; |
| 45 | ++ PetscInt supportSize; |
| 46 | ++ |
| 47 | ++ PetscCall(DMPlexGetSupportSize(dm, point, &supportSize)); |
| 48 | ++ if (supportSize <= key.part) continue; |
| 49 | ++ keptPoints[nKept++] = point; |
| 50 | ++ } |
| 51 | ++ PetscCall(ISRestoreIndices(pointIS, &allPoints)); |
| 52 | ++ PetscCall(ISDestroy(&pointIS)); |
| 53 | ++ PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nKept, keptPoints, PETSC_OWN_POINTER, &filteredIS)); |
| 54 | ++ pointIS = filteredIS; |
| 55 | ++ } |
| 56 | + PetscCall(ISGetLocalSize(pointIS, &numFaces)); |
| 57 | + PetscCall(ISGetIndices(pointIS, &points)); |
| 58 | + PetscCall(PetscMalloc4(numFaces * totDim, &u, (locX_t ? (size_t)numFaces * totDim : 0), &u_t, numFaces * totDim, &elemVec, (locA ? (size_t)numFaces * totDimAux : 0), &a)); |
| 59 | +@@ -4978,21 +5018,24 @@ PetscErrorCode DMPlexComputeBdResidualSingleByKey(DM dm, PetscWeakForm wf, Petsc |
| 60 | + for (face = 0; face < numFaces; ++face) { |
| 61 | + const PetscInt point = points[face], *support; |
| 62 | + PetscScalar *x = NULL; |
| 63 | +- PetscInt i; |
| 64 | ++ PetscInt i, supportSize, cell; |
| 65 | + |
| 66 | + PetscCall(DMPlexGetSupport(dm, point, &support)); |
| 67 | +- PetscCall(DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x)); |
| 68 | ++ PetscCall(DMPlexGetSupportSize(dm, point, &supportSize)); |
| 69 | ++ PetscCheck(supportSize > key.part, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Boundary point %" PetscInt_FMT " has support size %" PetscInt_FMT " but requested part %" PetscInt_FMT, point, supportSize, key.part); |
| 70 | ++ cell = support[key.part]; |
| 71 | ++ PetscCall(DMPlexVecGetClosure(plex, section, locX, cell, NULL, &x)); |
| 72 | + for (i = 0; i < totDim; ++i) u[face * totDim + i] = x[i]; |
| 73 | +- PetscCall(DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x)); |
| 74 | ++ PetscCall(DMPlexVecRestoreClosure(plex, section, locX, cell, NULL, &x)); |
| 75 | + if (locX_t) { |
| 76 | +- PetscCall(DMPlexVecGetClosure(plex, section, locX_t, support[0], NULL, &x)); |
| 77 | ++ PetscCall(DMPlexVecGetClosure(plex, section, locX_t, cell, NULL, &x)); |
| 78 | + for (i = 0; i < totDim; ++i) u_t[face * totDim + i] = x[i]; |
| 79 | +- PetscCall(DMPlexVecRestoreClosure(plex, section, locX_t, support[0], NULL, &x)); |
| 80 | ++ PetscCall(DMPlexVecRestoreClosure(plex, section, locX_t, cell, NULL, &x)); |
| 81 | + } |
| 82 | + if (locA) { |
| 83 | + PetscInt subp; |
| 84 | + |
| 85 | +- PetscCall(DMGetEnclosurePoint(plexA, dm, encAux, support[0], &subp)); |
| 86 | ++ PetscCall(DMGetEnclosurePoint(plexA, dm, encAux, cell, &subp)); |
| 87 | + PetscCall(DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x)); |
| 88 | + for (i = 0; i < totDimAux; ++i) a[face * totDimAux + i] = x[i]; |
| 89 | + PetscCall(DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x)); |
| 90 | +@@ -5028,10 +5071,13 @@ PetscErrorCode DMPlexComputeBdResidualSingleByKey(DM dm, PetscWeakForm wf, Petsc |
| 91 | + } |
| 92 | + for (face = 0; face < numFaces; ++face) { |
| 93 | + const PetscInt point = points[face], *support; |
| 94 | ++ PetscInt supportSize; |
| 95 | + |
| 96 | + if (mesh->printFEM > 1) PetscCall(DMPrintCellVector(point, name, totDim, &elemVec[face * totDim])); |
| 97 | + PetscCall(DMPlexGetSupport(plex, point, &support)); |
| 98 | +- PetscCall(DMPlexVecSetClosure(plex, NULL, locF, support[0], &elemVec[face * totDim], ADD_ALL_VALUES)); |
| 99 | ++ PetscCall(DMPlexGetSupportSize(dm, point, &supportSize)); |
| 100 | ++ PetscCheck(supportSize > key.part, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Boundary point %" PetscInt_FMT " has support size %" PetscInt_FMT " but requested part %" PetscInt_FMT, point, supportSize, key.part); |
| 101 | ++ PetscCall(DMPlexVecSetClosure(plex, NULL, locF, support[key.part], &elemVec[face * totDim], ADD_ALL_VALUES)); |
| 102 | + } |
| 103 | + PetscCall(DMSNESRestoreFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom)); |
| 104 | + PetscCall(PetscQuadratureDestroy(&qGeom)); |
| 105 | +@@ -5113,6 +5159,24 @@ static PetscErrorCode DMPlexComputeBdResidual_Internal(DM dm, Vec locX, Vec locX |
| 106 | + PetscCall(DMPlexGetDepthLabel(dm, &depthLabel)); |
| 107 | + PetscCall(DMGetDimension(dm, &dim)); |
| 108 | + PetscCall(DMLabelGetStratumIS(depthLabel, dim - 1, &facetIS)); |
| 109 | ++ /* Filter out ghost facets (SF leaves) so each facet contributes once */ |
| 110 | ++ if (facetIS) { |
| 111 | ++ PetscSF sf; |
| 112 | ++ PetscInt nleaves; |
| 113 | ++ const PetscInt *leaves; |
| 114 | ++ |
| 115 | ++ PetscCall(DMGetPointSF(dm, &sf)); |
| 116 | ++ PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL)); |
| 117 | ++ if (nleaves > 0 && leaves) { |
| 118 | ++ IS leafIS, ownedFacetIS; |
| 119 | ++ |
| 120 | ++ PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nleaves, leaves, PETSC_USE_POINTER, &leafIS)); |
| 121 | ++ PetscCall(ISDifference(facetIS, leafIS, &ownedFacetIS)); |
| 122 | ++ PetscCall(ISDestroy(&leafIS)); |
| 123 | ++ PetscCall(ISDestroy(&facetIS)); |
| 124 | ++ facetIS = ownedFacetIS; |
| 125 | ++ } |
| 126 | ++ } |
| 127 | + PetscCall(PetscDSGetNumBoundary(prob, &numBd)); |
| 128 | + for (bd = 0; bd < numBd; ++bd) { |
| 129 | + PetscWeakForm wf; |
| 130 | +@@ -6117,13 +6181,31 @@ static PetscErrorCode DMPlexComputeBdJacobian_Internal(DM dm, Vec locX, Vec locX |
| 131 | + PetscInt dim, numBd, bd; |
| 132 | + DMLabel depthLabel; |
| 133 | + DMField coordField = NULL; |
| 134 | +- IS facetIS; |
| 135 | ++ IS facetIS = NULL; |
| 136 | + |
| 137 | + PetscFunctionBegin; |
| 138 | + PetscCall(DMGetDS(dm, &prob)); |
| 139 | + PetscCall(DMPlexGetDepthLabel(dm, &depthLabel)); |
| 140 | + PetscCall(DMGetDimension(dm, &dim)); |
| 141 | + PetscCall(DMLabelGetStratumIS(depthLabel, dim - 1, &facetIS)); |
| 142 | ++ /* Filter out ghost facets (SF leaves) so each facet contributes once */ |
| 143 | ++ if (facetIS) { |
| 144 | ++ PetscSF sf; |
| 145 | ++ PetscInt nleaves; |
| 146 | ++ const PetscInt *leaves; |
| 147 | ++ |
| 148 | ++ PetscCall(DMGetPointSF(dm, &sf)); |
| 149 | ++ PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL)); |
| 150 | ++ if (nleaves > 0 && leaves) { |
| 151 | ++ IS leafIS, ownedFacetIS; |
| 152 | ++ |
| 153 | ++ PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nleaves, leaves, PETSC_USE_POINTER, &leafIS)); |
| 154 | ++ PetscCall(ISDifference(facetIS, leafIS, &ownedFacetIS)); |
| 155 | ++ PetscCall(ISDestroy(&leafIS)); |
| 156 | ++ PetscCall(ISDestroy(&facetIS)); |
| 157 | ++ facetIS = ownedFacetIS; |
| 158 | ++ } |
| 159 | ++ } |
| 160 | + PetscCall(PetscDSGetNumBoundary(prob, &numBd)); |
| 161 | + PetscCall(DMGetCoordinateField(dm, &coordField)); |
| 162 | + for (bd = 0; bd < numBd; ++bd) { |
0 commit comments