plexsubmesh.c (830e53ef5b16ac2200f70346ac25d1dd676ccb95) plexsubmesh.c (87feddfd2bc7dbf59fdfd7f81a9cb7f0c426fc5c)
1#include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2
3#undef __FUNCT__
4#define __FUNCT__ "DMPlexMarkBoundaryFaces"
5/*@
6 DMPlexMarkBoundaryFaces - Mark all faces on the boundary
7
8 Not Collective

--- 2077 unchanged lines hidden (view full) ---

2086static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, DM subdm)
2087{
2088 MPI_Comm comm;
2089 DMLabel subpointMap;
2090 IS subvertexIS;
2091 const PetscInt *subVertices;
2092 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells;
2093 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
1#include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2
3#undef __FUNCT__
4#define __FUNCT__ "DMPlexMarkBoundaryFaces"
5/*@
6 DMPlexMarkBoundaryFaces - Mark all faces on the boundary
7
8 Not Collective

--- 2077 unchanged lines hidden (view full) ---

2086static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, DM subdm)
2087{
2088 MPI_Comm comm;
2089 DMLabel subpointMap;
2090 IS subvertexIS;
2091 const PetscInt *subVertices;
2092 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells;
2093 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
2094 PetscInt vStart, vEnd, c, f;
2094 PetscInt cMax, c, f;
2095 PetscErrorCode ierr;
2096
2097 PetscFunctionBegin;
2098 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr);
2099 /* Create subpointMap which marks the submesh */
2100 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr);
2101 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr);
2102 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr);

--- 4 unchanged lines hidden (view full) ---

2107 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr);
2108 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr);
2109 /* Set cone sizes */
2110 firstSubVertex = numSubCells;
2111 firstSubFace = numSubCells+numSubVertices;
2112 newFacePoint = firstSubFace;
2113 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr);
2114 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);}
2095 PetscErrorCode ierr;
2096
2097 PetscFunctionBegin;
2098 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr);
2099 /* Create subpointMap which marks the submesh */
2100 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr);
2101 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr);
2102 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr);

--- 4 unchanged lines hidden (view full) ---

2107 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr);
2108 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr);
2109 /* Set cone sizes */
2110 firstSubVertex = numSubCells;
2111 firstSubFace = numSubCells+numSubVertices;
2112 newFacePoint = firstSubFace;
2113 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr);
2114 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);}
2115#if 0
2116 ierr = DMLabelGetStratumIS(subpointMap, 2, &subcellIS);CHKERRQ(ierr);
2117 if (subcellIS) {ierr = ISGetIndices(subcellIS, &subCells);CHKERRQ(ierr);}
2118#endif
2119 for (c = 0; c < numSubCells; ++c) {
2120 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr);
2121 }
2122 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2123 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr);
2124 }
2125 ierr = DMSetUp(subdm);CHKERRQ(ierr);
2126 /* Create face cones */
2115 for (c = 0; c < numSubCells; ++c) {
2116 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr);
2117 }
2118 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2119 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr);
2120 }
2121 ierr = DMSetUp(subdm);CHKERRQ(ierr);
2122 /* Create face cones */
2127 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
2128 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr);
2123 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr);
2124 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr);
2129 ierr = DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr);
2130 for (c = 0; c < numSubCells; ++c) {
2125 ierr = DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr);
2126 for (c = 0; c < numSubCells; ++c) {
2131 const PetscInt cell = subCells[c];
2132 const PetscInt subcell = c;
2133 PetscInt *closure = NULL;
2134 PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
2127 const PetscInt cell = subCells[c];
2128 const PetscInt subcell = c;
2129 const PetscInt *cone, *cells;
2130 PetscInt numCells, subVertex, p, v;
2135
2131
2136 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
2137 for (cl = 0; cl < closureSize*2; cl += 2) {
2138 const PetscInt point = closure[cl];
2139 PetscInt subVertex;
2132 if (cell < cMax) continue;
2133 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr);
2134 for (v = 0; v < nFV; ++v) {
2135 ierr = PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);CHKERRQ(ierr);
2136 subface[v] = firstSubVertex+subVertex;
2137 }
2138 ierr = DMPlexSetCone(subdm, newFacePoint, subface);CHKERRQ(ierr);
2139 ierr = DMPlexSetCone(subdm, subcell, &newFacePoint);CHKERRQ(ierr);
2140 ierr = DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr);
2141 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2142 for (p = 0; p < numCells; ++p) {
2143 PetscInt negsubcell;
2140
2144
2141 if ((point >= vStart) && (point < vEnd)) {
2142 ++numCorners;
2143 ierr = PetscFindInt(point, numSubVertices, subVertices, &subVertex);CHKERRQ(ierr);
2144 if (subVertex >= 0) {
2145 closure[faceSize] = point;
2146 subface[faceSize] = firstSubVertex+subVertex;
2147 ++faceSize;
2148 }
2145 if (cells[p] >= cMax) continue;
2146 /* I know this is a crap search */
2147 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
2148 if (subCells[negsubcell] == cells[p]) break;
2149 }
2149 }
2150 if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
2151 ierr = DMPlexSetCone(subdm, negsubcell, &newFacePoint);CHKERRQ(ierr);
2150 }
2152 }
2151 if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2152 if (faceSize == nFV) {
2153 ierr = DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);CHKERRQ(ierr);
2154 }
2155 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
2153 ierr = DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr);
2154 ++newFacePoint;
2156 }
2157 ierr = DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr);
2158 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr);
2159 ierr = DMPlexStratify(subdm);CHKERRQ(ierr);
2160 /* Build coordinates */
2161 {
2162 PetscSection coordSection, subCoordSection;
2163 Vec coordinates, subCoordinates;

--- 36 unchanged lines hidden (view full) ---

2200 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr);
2201 }
2202 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr);
2203 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr);
2204 }
2205 /* Cleanup */
2206 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);}
2207 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr);
2155 }
2156 ierr = DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr);
2157 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr);
2158 ierr = DMPlexStratify(subdm);CHKERRQ(ierr);
2159 /* Build coordinates */
2160 {
2161 PetscSection coordSection, subCoordSection;
2162 Vec coordinates, subCoordinates;

--- 36 unchanged lines hidden (view full) ---

2199 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr);
2200 }
2201 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr);
2202 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr);
2203 }
2204 /* Cleanup */
2205 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);}
2206 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr);
2208#if 0
2209 if (subcellIS) {ierr = ISRestoreIndices(subcellIS, &subCells);CHKERRQ(ierr);}
2210 ierr = ISDestroy(&subcellIS);CHKERRQ(ierr);
2211#else
2212 ierr = PetscFree(subCells);CHKERRQ(ierr);
2207 ierr = PetscFree(subCells);CHKERRQ(ierr);
2213#endif
2214 PetscFunctionReturn(0);
2215}
2216
2217#undef __FUNCT__
2218#define __FUNCT__ "DMPlexCreateCohesiveSubmesh_Interpolated"
2219static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, PetscBool hasLagrange, DM subdm)
2220{
2221 MPI_Comm comm;

--- 266 unchanged lines hidden --- 		2208 PetscFunctionReturn(0);
2209}
2210
2211#undef __FUNCT__
2212#define __FUNCT__ "DMPlexCreateCohesiveSubmesh_Interpolated"
2213static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, PetscBool hasLagrange, DM subdm)
2214{
2215 MPI_Comm comm;

--- 266 unchanged lines hidden ---