xref: /petsc/src/dm/impls/plex/plexreorder.c (revision d02c7345ebed081476ce6bcdfd784b34c024d6ef)
1af0996ceSBarry Smith #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2af0996ceSBarry Smith #include <petsc/private/matorderimpl.h> /*I      "petscmat.h"      I*/
3695799ffSMatthew G. Knepley #include <petsc/private/dmlabelimpl.h>
48ed5f475SMatthew G. Knepley 
5d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMPlexCreateOrderingClosure_Static(DM dm, PetscInt numPoints, const PetscInt pperm[], PetscInt **clperm, PetscInt **invclperm)
6d71ae5a4SJacob Faibussowitsch {
78ed5f475SMatthew G. Knepley   PetscInt *perm, *iperm;
88ed5f475SMatthew G. Knepley   PetscInt  depth, d, pStart, pEnd, fStart, fMax, fEnd, p;
98ed5f475SMatthew G. Knepley 
108ed5f475SMatthew G. Knepley   PetscFunctionBegin;
119566063dSJacob Faibussowitsch   PetscCall(DMPlexGetDepth(dm, &depth));
129566063dSJacob Faibussowitsch   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
139566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(pEnd - pStart, &perm));
149566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(pEnd - pStart, &iperm));
158ed5f475SMatthew G. Knepley   for (p = pStart; p < pEnd; ++p) iperm[p] = -1;
168ed5f475SMatthew G. Knepley   for (d = depth; d > 0; --d) {
179566063dSJacob Faibussowitsch     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
189566063dSJacob Faibussowitsch     PetscCall(DMPlexGetDepthStratum(dm, d - 1, &fStart, &fEnd));
198ed5f475SMatthew G. Knepley     fMax = fStart;
208ed5f475SMatthew G. Knepley     for (p = pStart; p < pEnd; ++p) {
218ed5f475SMatthew G. Knepley       const PetscInt *cone;
228ed5f475SMatthew G. Knepley       PetscInt        point, coneSize, c;
238ed5f475SMatthew G. Knepley 
248ed5f475SMatthew G. Knepley       if (d == depth) {
258ed5f475SMatthew G. Knepley         perm[p]         = pperm[p];
268ed5f475SMatthew G. Knepley         iperm[pperm[p]] = p;
278ed5f475SMatthew G. Knepley       }
288ed5f475SMatthew G. Knepley       point = perm[p];
299566063dSJacob Faibussowitsch       PetscCall(DMPlexGetConeSize(dm, point, &coneSize));
309566063dSJacob Faibussowitsch       PetscCall(DMPlexGetCone(dm, point, &cone));
318ed5f475SMatthew G. Knepley       for (c = 0; c < coneSize; ++c) {
328ed5f475SMatthew G. Knepley         const PetscInt oldc = cone[c];
338ed5f475SMatthew G. Knepley         const PetscInt newc = iperm[oldc];
348ed5f475SMatthew G. Knepley 
358ed5f475SMatthew G. Knepley         if (newc < 0) {
368ed5f475SMatthew G. Knepley           perm[fMax]  = oldc;
378ed5f475SMatthew G. Knepley           iperm[oldc] = fMax++;
388ed5f475SMatthew G. Knepley         }
398ed5f475SMatthew G. Knepley       }
408ed5f475SMatthew G. Knepley     }
4163a3b9bcSJacob Faibussowitsch     PetscCheck(fMax == fEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of depth %" PetscInt_FMT " faces %" PetscInt_FMT " does not match permuted number %" PetscInt_FMT, d, fEnd - fStart, fMax - fStart);
428ed5f475SMatthew G. Knepley   }
438ed5f475SMatthew G. Knepley   *clperm    = perm;
448ed5f475SMatthew G. Knepley   *invclperm = iperm;
453ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
468ed5f475SMatthew G. Knepley }
478ed5f475SMatthew G. Knepley 
4860225df5SJacob Faibussowitsch /*@C
498ed5f475SMatthew G. Knepley   DMPlexGetOrdering - Calculate a reordering of the mesh
508ed5f475SMatthew G. Knepley 
5120f4b53cSBarry Smith   Collective
528ed5f475SMatthew G. Knepley 
53d8d19677SJose E. Roman   Input Parameters:
548ed5f475SMatthew G. Knepley + dm    - The DMPlex object
5520f4b53cSBarry Smith . otype - type of reordering, see `MatOrderingType`
5620f4b53cSBarry Smith - label - [Optional] Label used to segregate ordering into sets, or `NULL`
578ed5f475SMatthew G. Knepley 
588ed5f475SMatthew G. Knepley   Output Parameter:
5920f4b53cSBarry Smith . perm - The point permutation as an `IS`, `perm`[old point number] = new point number
608ed5f475SMatthew G. Knepley 
618ed5f475SMatthew G. Knepley   Level: intermediate
628ed5f475SMatthew G. Knepley 
6320f4b53cSBarry Smith   Note:
6420f4b53cSBarry Smith   The label is used to group sets of points together by label value. This makes it easy to reorder a mesh which
6520f4b53cSBarry Smith   has different types of cells, and then loop over each set of reordered cells for assembly.
6620f4b53cSBarry Smith 
6720f4b53cSBarry Smith .seealso: `DMPLEX`, `DMPlexPermute()`, `MatOrderingType`, `MatGetOrdering()`
688ed5f475SMatthew G. Knepley @*/
69d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm)
70d71ae5a4SJacob Faibussowitsch {
718ed5f475SMatthew G. Knepley   PetscInt  numCells = 0;
723c67d5adSSatish Balay   PetscInt *start = NULL, *adjacency = NULL, *cperm, *clperm = NULL, *invclperm = NULL, *mask, *xls, pStart, pEnd, c, i;
738ed5f475SMatthew G. Knepley 
748ed5f475SMatthew G. Knepley   PetscFunctionBegin;
758ed5f475SMatthew G. Knepley   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
764f572ea9SToby Isaac   PetscAssertPointer(perm, 4);
779566063dSJacob Faibussowitsch   PetscCall(DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency));
789566063dSJacob Faibussowitsch   PetscCall(PetscMalloc3(numCells, &cperm, numCells, &mask, numCells * 2, &xls));
79d80c2872SMatthew G. Knepley   if (numCells) {
808ed5f475SMatthew G. Knepley     /* Shift for Fortran numbering */
818ed5f475SMatthew G. Knepley     for (i = 0; i < start[numCells]; ++i) ++adjacency[i];
828ed5f475SMatthew G. Knepley     for (i = 0; i <= numCells; ++i) ++start[i];
839566063dSJacob Faibussowitsch     PetscCall(SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls));
84d80c2872SMatthew G. Knepley   }
859566063dSJacob Faibussowitsch   PetscCall(PetscFree(start));
869566063dSJacob Faibussowitsch   PetscCall(PetscFree(adjacency));
878ed5f475SMatthew G. Knepley   /* Shift for Fortran numbering */
888ed5f475SMatthew G. Knepley   for (c = 0; c < numCells; ++c) --cperm[c];
89c99f2573SMatthew G. Knepley   /* Segregate */
90c99f2573SMatthew G. Knepley   if (label) {
91c99f2573SMatthew G. Knepley     IS              valueIS;
9250be1cf5SMatthew G. Knepley     const PetscInt *valuesTmp;
9350be1cf5SMatthew G. Knepley     PetscInt       *values;
94c99f2573SMatthew G. Knepley     PetscInt        numValues, numPoints = 0;
95c99f2573SMatthew G. Knepley     PetscInt       *sperm, *vsize, *voff, v;
96c99f2573SMatthew G. Knepley 
9750be1cf5SMatthew G. Knepley     // Can't directly sort the valueIS, since it is a view into the DMLabel
989566063dSJacob Faibussowitsch     PetscCall(DMLabelGetValueIS(label, &valueIS));
999566063dSJacob Faibussowitsch     PetscCall(ISGetLocalSize(valueIS, &numValues));
10050be1cf5SMatthew G. Knepley     PetscCall(ISGetIndices(valueIS, &valuesTmp));
10150be1cf5SMatthew G. Knepley     PetscCall(PetscCalloc4(numCells, &sperm, numValues, &values, numValues, &vsize, numValues + 1, &voff));
10250be1cf5SMatthew G. Knepley     PetscCall(PetscArraycpy(values, valuesTmp, numValues));
10350be1cf5SMatthew G. Knepley     PetscCall(PetscSortInt(numValues, values));
10450be1cf5SMatthew G. Knepley     PetscCall(ISRestoreIndices(valueIS, &valuesTmp));
10550be1cf5SMatthew G. Knepley     PetscCall(ISDestroy(&valueIS));
106c99f2573SMatthew G. Knepley     for (v = 0; v < numValues; ++v) {
1079566063dSJacob Faibussowitsch       PetscCall(DMLabelGetStratumSize(label, values[v], &vsize[v]));
108c99f2573SMatthew G. Knepley       if (v < numValues - 1) voff[v + 2] += vsize[v] + voff[v + 1];
109c99f2573SMatthew G. Knepley       numPoints += vsize[v];
110c99f2573SMatthew G. Knepley     }
1116469bafaSMatthew G. Knepley     PetscCheck(numPoints == numCells, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label only covers %" PetscInt_FMT " cells != %" PetscInt_FMT " total cells", numPoints, numCells);
112c99f2573SMatthew G. Knepley     for (c = 0; c < numCells; ++c) {
113c99f2573SMatthew G. Knepley       const PetscInt oldc = cperm[c];
114c99f2573SMatthew G. Knepley       PetscInt       val, vloc;
115c99f2573SMatthew G. Knepley 
1169566063dSJacob Faibussowitsch       PetscCall(DMLabelGetValue(label, oldc, &val));
11763a3b9bcSJacob Faibussowitsch       PetscCheck(val != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %" PetscInt_FMT " not present in label", oldc);
1189566063dSJacob Faibussowitsch       PetscCall(PetscFindInt(val, numValues, values, &vloc));
11963a3b9bcSJacob Faibussowitsch       PetscCheck(vloc >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %" PetscInt_FMT " not present label", val);
120c99f2573SMatthew G. Knepley       sperm[voff[vloc + 1]++] = oldc;
121c99f2573SMatthew G. Knepley     }
122ad540459SPierre Jolivet     for (v = 0; v < numValues; ++v) PetscCheck(voff[v + 1] - voff[v] == vsize[v], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of %" PetscInt_FMT " values found is %" PetscInt_FMT " != %" PetscInt_FMT, values[v], voff[v + 1] - voff[v], vsize[v]);
1239566063dSJacob Faibussowitsch     PetscCall(PetscArraycpy(cperm, sperm, numCells));
12450be1cf5SMatthew G. Knepley     PetscCall(PetscFree4(sperm, values, vsize, voff));
125c99f2573SMatthew G. Knepley   }
1268ed5f475SMatthew G. Knepley   /* Construct closure */
1279566063dSJacob Faibussowitsch   PetscCall(DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm));
1289566063dSJacob Faibussowitsch   PetscCall(PetscFree3(cperm, mask, xls));
1299566063dSJacob Faibussowitsch   PetscCall(PetscFree(clperm));
1308ed5f475SMatthew G. Knepley   /* Invert permutation */
1319566063dSJacob Faibussowitsch   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1329566063dSJacob Faibussowitsch   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, invclperm, PETSC_OWN_POINTER, perm));
1333ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1348ed5f475SMatthew G. Knepley }
1358ed5f475SMatthew G. Knepley 
1368ed5f475SMatthew G. Knepley /*@
1376913077dSMatthew G. Knepley   DMPlexGetOrdering1D - Reorder the vertices so that the mesh is in a line
1386913077dSMatthew G. Knepley 
13920f4b53cSBarry Smith   Collective
1406913077dSMatthew G. Knepley 
1416913077dSMatthew G. Knepley   Input Parameter:
14220f4b53cSBarry Smith . dm - The `DMPLEX` object
1436913077dSMatthew G. Knepley 
1446913077dSMatthew G. Knepley   Output Parameter:
14520f4b53cSBarry Smith . perm - The point permutation as an `IS`, `perm`[old point number] = new point number
1466913077dSMatthew G. Knepley 
1476913077dSMatthew G. Knepley   Level: intermediate
1486913077dSMatthew G. Knepley 
14920f4b53cSBarry Smith .seealso: `DMPLEX`, `DMPlexGetOrdering()`, `DMPlexPermute()`, `MatGetOrdering()`
1506913077dSMatthew G. Knepley @*/
151d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexGetOrdering1D(DM dm, IS *perm)
152d71ae5a4SJacob Faibussowitsch {
1536913077dSMatthew G. Knepley   PetscInt       *points;
1546913077dSMatthew G. Knepley   const PetscInt *support, *cone;
1556913077dSMatthew G. Knepley   PetscInt        dim, pStart, pEnd, cStart, cEnd, c, vStart, vEnd, v, suppSize, lastCell = 0;
1566913077dSMatthew G. Knepley 
1576913077dSMatthew G. Knepley   PetscFunctionBegin;
1589566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(dm, &dim));
1596913077dSMatthew G. Knepley   PetscCheck(dim == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Input mesh must be one dimensional, not %" PetscInt_FMT, dim);
1609566063dSJacob Faibussowitsch   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1619566063dSJacob Faibussowitsch   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1629566063dSJacob Faibussowitsch   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1639566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(pEnd - pStart, &points));
1646913077dSMatthew G. Knepley   for (c = cStart; c < cEnd; ++c) points[c] = c;
1656913077dSMatthew G. Knepley   for (v = vStart; v < vEnd; ++v) points[v] = v;
1666913077dSMatthew G. Knepley   for (v = vStart; v < vEnd; ++v) {
1679566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupportSize(dm, v, &suppSize));
1689566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupport(dm, v, &support));
1699371c9d4SSatish Balay     if (suppSize == 1) {
1709371c9d4SSatish Balay       lastCell = support[0];
1719371c9d4SSatish Balay       break;
1729371c9d4SSatish Balay     }
1736913077dSMatthew G. Knepley   }
1746913077dSMatthew G. Knepley   if (v < vEnd) {
1756913077dSMatthew G. Knepley     PetscInt pos = cEnd;
1766913077dSMatthew G. Knepley 
1776913077dSMatthew G. Knepley     points[v] = pos++;
1786913077dSMatthew G. Knepley     while (lastCell >= cStart) {
1799566063dSJacob Faibussowitsch       PetscCall(DMPlexGetCone(dm, lastCell, &cone));
1806913077dSMatthew G. Knepley       if (cone[0] == v) v = cone[1];
1816913077dSMatthew G. Knepley       else v = cone[0];
1829566063dSJacob Faibussowitsch       PetscCall(DMPlexGetSupport(dm, v, &support));
1839566063dSJacob Faibussowitsch       PetscCall(DMPlexGetSupportSize(dm, v, &suppSize));
1849371c9d4SSatish Balay       if (suppSize == 1) {
1859371c9d4SSatish Balay         lastCell = -1;
1869371c9d4SSatish Balay       } else {
1876913077dSMatthew G. Knepley         if (support[0] == lastCell) lastCell = support[1];
1886913077dSMatthew G. Knepley         else lastCell = support[0];
1896913077dSMatthew G. Knepley       }
1906913077dSMatthew G. Knepley       points[v] = pos++;
1916913077dSMatthew G. Knepley     }
1926913077dSMatthew G. Knepley     PetscCheck(pos == pEnd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Last vertex was %" PetscInt_FMT ", not %" PetscInt_FMT, pos, pEnd);
1936913077dSMatthew G. Knepley   }
1949566063dSJacob Faibussowitsch   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, points, PETSC_OWN_POINTER, perm));
1953ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1966913077dSMatthew G. Knepley }
1976913077dSMatthew G. Knepley 
198d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMPlexRemapCoordinates_Private(IS perm, PetscSection cs, Vec coordinates, PetscSection *csNew, Vec *coordinatesNew)
199d71ae5a4SJacob Faibussowitsch {
2006858538eSMatthew G. Knepley   PetscScalar    *coords, *coordsNew;
2016858538eSMatthew G. Knepley   const PetscInt *pperm;
2026858538eSMatthew G. Knepley   PetscInt        pStart, pEnd, p;
2036858538eSMatthew G. Knepley   const char     *name;
2046858538eSMatthew G. Knepley 
2056858538eSMatthew G. Knepley   PetscFunctionBegin;
2066858538eSMatthew G. Knepley   PetscCall(PetscSectionPermute(cs, perm, csNew));
2076858538eSMatthew G. Knepley   PetscCall(VecDuplicate(coordinates, coordinatesNew));
2086858538eSMatthew G. Knepley   PetscCall(PetscObjectGetName((PetscObject)coordinates, &name));
2096858538eSMatthew G. Knepley   PetscCall(PetscObjectSetName((PetscObject)*coordinatesNew, name));
2106858538eSMatthew G. Knepley   PetscCall(VecGetArray(coordinates, &coords));
2116858538eSMatthew G. Knepley   PetscCall(VecGetArray(*coordinatesNew, &coordsNew));
2126858538eSMatthew G. Knepley   PetscCall(PetscSectionGetChart(*csNew, &pStart, &pEnd));
2136858538eSMatthew G. Knepley   PetscCall(ISGetIndices(perm, &pperm));
2146858538eSMatthew G. Knepley   for (p = pStart; p < pEnd; ++p) {
2156858538eSMatthew G. Knepley     PetscInt dof, off, offNew, d;
2166858538eSMatthew G. Knepley 
2176858538eSMatthew G. Knepley     PetscCall(PetscSectionGetDof(*csNew, p, &dof));
2186858538eSMatthew G. Knepley     PetscCall(PetscSectionGetOffset(cs, p, &off));
2196858538eSMatthew G. Knepley     PetscCall(PetscSectionGetOffset(*csNew, pperm[p], &offNew));
2206858538eSMatthew G. Knepley     for (d = 0; d < dof; ++d) coordsNew[offNew + d] = coords[off + d];
2216858538eSMatthew G. Knepley   }
2226858538eSMatthew G. Knepley   PetscCall(ISRestoreIndices(perm, &pperm));
2236858538eSMatthew G. Knepley   PetscCall(VecRestoreArray(coordinates, &coords));
2246858538eSMatthew G. Knepley   PetscCall(VecRestoreArray(*coordinatesNew, &coordsNew));
2253ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
2266858538eSMatthew G. Knepley }
2276858538eSMatthew G. Knepley 
2286913077dSMatthew G. Knepley /*@
2298ed5f475SMatthew G. Knepley   DMPlexPermute - Reorder the mesh according to the input permutation
2308ed5f475SMatthew G. Knepley 
23120f4b53cSBarry Smith   Collective
2328ed5f475SMatthew G. Knepley 
233d8d19677SJose E. Roman   Input Parameters:
23420f4b53cSBarry Smith + dm   - The `DMPLEX` object
23520f4b53cSBarry Smith - perm - The point permutation, `perm`[old point number] = new point number
2368ed5f475SMatthew G. Knepley 
2378ed5f475SMatthew G. Knepley   Output Parameter:
23820f4b53cSBarry Smith . pdm - The permuted `DM`
2398ed5f475SMatthew G. Knepley 
2408ed5f475SMatthew G. Knepley   Level: intermediate
2418ed5f475SMatthew G. Knepley 
24220f4b53cSBarry Smith .seealso: `DMPLEX`, `MatPermute()`
2438ed5f475SMatthew G. Knepley @*/
244d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
245d71ae5a4SJacob Faibussowitsch {
2468ed5f475SMatthew G. Knepley   DM_Plex    *plex = (DM_Plex *)dm->data, *plexNew;
2472f300c92SMatthew Knepley   PetscInt    dim, cdim;
2486b74b800SMatthew G. Knepley   const char *name;
2498ed5f475SMatthew G. Knepley 
2508ed5f475SMatthew G. Knepley   PetscFunctionBegin;
2518ed5f475SMatthew G. Knepley   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
2528ed5f475SMatthew G. Knepley   PetscValidHeaderSpecific(perm, IS_CLASSID, 2);
2534f572ea9SToby Isaac   PetscAssertPointer(pdm, 3);
2549566063dSJacob Faibussowitsch   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), pdm));
2559566063dSJacob Faibussowitsch   PetscCall(DMSetType(*pdm, DMPLEX));
2569566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetName((PetscObject)dm, &name));
2579566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)*pdm, name));
2589566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(dm, &dim));
2599566063dSJacob Faibussowitsch   PetscCall(DMSetDimension(*pdm, dim));
2609566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateDim(dm, &cdim));
2619566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinateDim(*pdm, cdim));
2629566063dSJacob Faibussowitsch   PetscCall(DMCopyDisc(dm, *pdm));
2636b74b800SMatthew G. Knepley   if (dm->localSection) {
2646b74b800SMatthew G. Knepley     PetscSection section, sectionNew;
2656b74b800SMatthew G. Knepley 
2669566063dSJacob Faibussowitsch     PetscCall(DMGetLocalSection(dm, &section));
2679566063dSJacob Faibussowitsch     PetscCall(PetscSectionPermute(section, perm, &sectionNew));
2689566063dSJacob Faibussowitsch     PetscCall(DMSetLocalSection(*pdm, sectionNew));
2699566063dSJacob Faibussowitsch     PetscCall(PetscSectionDestroy(&sectionNew));
270dd742cf1SMatthew G. Knepley   }
2718ed5f475SMatthew G. Knepley   plexNew = (DM_Plex *)(*pdm)->data;
2728ed5f475SMatthew G. Knepley   /* Ignore ltogmap, ltogmapb */
2731bb6d2a8SBarry Smith   /* Ignore sf, sectionSF */
2748ed5f475SMatthew G. Knepley   /* Ignore globalVertexNumbers, globalCellNumbers */
2758ed5f475SMatthew G. Knepley   /* Reorder labels */
2768ed5f475SMatthew G. Knepley   {
277a85475f2SMatthew G. Knepley     PetscInt numLabels, l;
278a85475f2SMatthew G. Knepley     DMLabel  label, labelNew;
2798ed5f475SMatthew G. Knepley 
2809566063dSJacob Faibussowitsch     PetscCall(DMGetNumLabels(dm, &numLabels));
2815d80c0bfSVaclav Hapla     for (l = 0; l < numLabels; ++l) {
2829566063dSJacob Faibussowitsch       PetscCall(DMGetLabelByNum(dm, l, &label));
2839566063dSJacob Faibussowitsch       PetscCall(DMLabelPermute(label, perm, &labelNew));
2849566063dSJacob Faibussowitsch       PetscCall(DMAddLabel(*pdm, labelNew));
2859566063dSJacob Faibussowitsch       PetscCall(DMLabelDestroy(&labelNew));
2868ed5f475SMatthew G. Knepley     }
2879566063dSJacob Faibussowitsch     PetscCall(DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel));
2889566063dSJacob Faibussowitsch     if (plex->subpointMap) PetscCall(DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap));
2898ed5f475SMatthew G. Knepley   }
290695799ffSMatthew G. Knepley   if ((*pdm)->celltypeLabel) {
291695799ffSMatthew G. Knepley     DMLabel ctLabel;
292695799ffSMatthew G. Knepley 
293695799ffSMatthew G. Knepley     // Reset label for fast lookup
294695799ffSMatthew G. Knepley     PetscCall(DMPlexGetCellTypeLabel(*pdm, &ctLabel));
295695799ffSMatthew G. Knepley     PetscCall(DMLabelMakeAllInvalid_Internal(ctLabel));
296695799ffSMatthew G. Knepley   }
2978ed5f475SMatthew G. Knepley   /* Reorder topology */
2988ed5f475SMatthew G. Knepley   {
2998ed5f475SMatthew G. Knepley     const PetscInt *pperm;
3006302a7fbSVaclav Hapla     PetscInt        n, pStart, pEnd, p;
3018ed5f475SMatthew G. Knepley 
3029566063dSJacob Faibussowitsch     PetscCall(PetscSectionDestroy(&plexNew->coneSection));
3039566063dSJacob Faibussowitsch     PetscCall(PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection));
3049566063dSJacob Faibussowitsch     PetscCall(PetscSectionGetStorageSize(plexNew->coneSection, &n));
3059566063dSJacob Faibussowitsch     PetscCall(PetscMalloc1(n, &plexNew->cones));
3069566063dSJacob Faibussowitsch     PetscCall(PetscMalloc1(n, &plexNew->coneOrientations));
3079566063dSJacob Faibussowitsch     PetscCall(ISGetIndices(perm, &pperm));
3089566063dSJacob Faibussowitsch     PetscCall(PetscSectionGetChart(plex->coneSection, &pStart, &pEnd));
3098ed5f475SMatthew G. Knepley     for (p = pStart; p < pEnd; ++p) {
3108ed5f475SMatthew G. Knepley       PetscInt dof, off, offNew, d;
3118ed5f475SMatthew G. Knepley 
3129566063dSJacob Faibussowitsch       PetscCall(PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof));
3139566063dSJacob Faibussowitsch       PetscCall(PetscSectionGetOffset(plex->coneSection, p, &off));
3149566063dSJacob Faibussowitsch       PetscCall(PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew));
3158ed5f475SMatthew G. Knepley       for (d = 0; d < dof; ++d) {
3168ed5f475SMatthew G. Knepley         plexNew->cones[offNew + d]            = pperm[plex->cones[off + d]];
3178ed5f475SMatthew G. Knepley         plexNew->coneOrientations[offNew + d] = plex->coneOrientations[off + d];
3188ed5f475SMatthew G. Knepley       }
3198ed5f475SMatthew G. Knepley     }
3209566063dSJacob Faibussowitsch     PetscCall(PetscSectionDestroy(&plexNew->supportSection));
3219566063dSJacob Faibussowitsch     PetscCall(PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection));
3229566063dSJacob Faibussowitsch     PetscCall(PetscSectionGetStorageSize(plexNew->supportSection, &n));
3239566063dSJacob Faibussowitsch     PetscCall(PetscMalloc1(n, &plexNew->supports));
3249566063dSJacob Faibussowitsch     PetscCall(PetscSectionGetChart(plex->supportSection, &pStart, &pEnd));
3258ed5f475SMatthew G. Knepley     for (p = pStart; p < pEnd; ++p) {
3268ed5f475SMatthew G. Knepley       PetscInt dof, off, offNew, d;
3278ed5f475SMatthew G. Knepley 
3289566063dSJacob Faibussowitsch       PetscCall(PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof));
3299566063dSJacob Faibussowitsch       PetscCall(PetscSectionGetOffset(plex->supportSection, p, &off));
3309566063dSJacob Faibussowitsch       PetscCall(PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew));
331ad540459SPierre Jolivet       for (d = 0; d < dof; ++d) plexNew->supports[offNew + d] = pperm[plex->supports[off + d]];
3328ed5f475SMatthew G. Knepley     }
3339566063dSJacob Faibussowitsch     PetscCall(ISRestoreIndices(perm, &pperm));
3348ed5f475SMatthew G. Knepley   }
335a6e0b375SMatthew G. Knepley   /* Remap coordinates */
336a6e0b375SMatthew G. Knepley   {
337a6e0b375SMatthew G. Knepley     DM           cdm, cdmNew;
3386858538eSMatthew G. Knepley     PetscSection cs, csNew;
339a6e0b375SMatthew G. Knepley     Vec          coordinates, coordinatesNew;
340a6e0b375SMatthew G. Knepley 
3416858538eSMatthew G. Knepley     PetscCall(DMGetCoordinateSection(dm, &cs));
3429566063dSJacob Faibussowitsch     PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
3436858538eSMatthew G. Knepley     PetscCall(DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew));
3446858538eSMatthew G. Knepley     PetscCall(DMSetCoordinateSection(*pdm, PETSC_DETERMINE, csNew));
3459566063dSJacob Faibussowitsch     PetscCall(DMSetCoordinatesLocal(*pdm, coordinatesNew));
3466858538eSMatthew G. Knepley     PetscCall(PetscSectionDestroy(&csNew));
3479566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&coordinatesNew));
3486858538eSMatthew G. Knepley 
3496858538eSMatthew G. Knepley     PetscCall(DMGetCellCoordinateDM(dm, &cdm));
3506858538eSMatthew G. Knepley     if (cdm) {
3516858538eSMatthew G. Knepley       PetscCall(DMGetCoordinateDM(*pdm, &cdm));
3526858538eSMatthew G. Knepley       PetscCall(DMClone(cdm, &cdmNew));
3536858538eSMatthew G. Knepley       PetscCall(DMSetCellCoordinateDM(*pdm, cdmNew));
3546858538eSMatthew G. Knepley       PetscCall(DMDestroy(&cdmNew));
3556858538eSMatthew G. Knepley       PetscCall(DMGetCellCoordinateSection(dm, &cs));
3566858538eSMatthew G. Knepley       PetscCall(DMGetCellCoordinatesLocal(dm, &coordinates));
3576858538eSMatthew G. Knepley       PetscCall(DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew));
3586858538eSMatthew G. Knepley       PetscCall(DMSetCellCoordinateSection(*pdm, PETSC_DETERMINE, csNew));
3596858538eSMatthew G. Knepley       PetscCall(DMSetCellCoordinatesLocal(*pdm, coordinatesNew));
3606858538eSMatthew G. Knepley       PetscCall(PetscSectionDestroy(&csNew));
3616858538eSMatthew G. Knepley       PetscCall(VecDestroy(&coordinatesNew));
3626858538eSMatthew G. Knepley     }
363a6e0b375SMatthew G. Knepley   }
3645de52c6dSVaclav Hapla   PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *pdm));
365e66f2fa0SMatthew G. Knepley   (*pdm)->setupcalled = PETSC_TRUE;
3663ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3678ed5f475SMatthew G. Knepley }
3686bc1bd01Sksagiyam 
369d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexReorderSetDefault_Plex(DM dm, DMPlexReorderDefaultFlag reorder)
370d71ae5a4SJacob Faibussowitsch {
3716bc1bd01Sksagiyam   DM_Plex *mesh = (DM_Plex *)dm->data;
3726bc1bd01Sksagiyam 
3736bc1bd01Sksagiyam   PetscFunctionBegin;
3746bc1bd01Sksagiyam   mesh->reorderDefault = reorder;
3753ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3766bc1bd01Sksagiyam }
3776bc1bd01Sksagiyam 
3786bc1bd01Sksagiyam /*@
3796bc1bd01Sksagiyam   DMPlexReorderSetDefault - Set flag indicating whether the DM should be reordered by default
3806bc1bd01Sksagiyam 
38120f4b53cSBarry Smith   Logically Collective
3826bc1bd01Sksagiyam 
3836bc1bd01Sksagiyam   Input Parameters:
38420f4b53cSBarry Smith + dm      - The `DM`
3856bc1bd01Sksagiyam - reorder - Flag for reordering
3866bc1bd01Sksagiyam 
3876bc1bd01Sksagiyam   Level: intermediate
3886bc1bd01Sksagiyam 
3896bc1bd01Sksagiyam .seealso: `DMPlexReorderGetDefault()`
3906bc1bd01Sksagiyam @*/
391d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexReorderSetDefault(DM dm, DMPlexReorderDefaultFlag reorder)
392d71ae5a4SJacob Faibussowitsch {
3936bc1bd01Sksagiyam   PetscFunctionBegin;
3946bc1bd01Sksagiyam   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
3956bc1bd01Sksagiyam   PetscTryMethod(dm, "DMPlexReorderSetDefault_C", (DM, DMPlexReorderDefaultFlag), (dm, reorder));
3963ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3976bc1bd01Sksagiyam }
3986bc1bd01Sksagiyam 
399d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexReorderGetDefault_Plex(DM dm, DMPlexReorderDefaultFlag *reorder)
400d71ae5a4SJacob Faibussowitsch {
4016bc1bd01Sksagiyam   DM_Plex *mesh = (DM_Plex *)dm->data;
4026bc1bd01Sksagiyam 
4036bc1bd01Sksagiyam   PetscFunctionBegin;
4046bc1bd01Sksagiyam   *reorder = mesh->reorderDefault;
4053ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
4066bc1bd01Sksagiyam }
4076bc1bd01Sksagiyam 
4086bc1bd01Sksagiyam /*@
4096bc1bd01Sksagiyam   DMPlexReorderGetDefault - Get flag indicating whether the DM should be reordered by default
4106bc1bd01Sksagiyam 
41120f4b53cSBarry Smith   Not Collective
4126bc1bd01Sksagiyam 
4136bc1bd01Sksagiyam   Input Parameter:
41420f4b53cSBarry Smith . dm - The `DM`
4156bc1bd01Sksagiyam 
4166bc1bd01Sksagiyam   Output Parameter:
4176bc1bd01Sksagiyam . reorder - Flag for reordering
4186bc1bd01Sksagiyam 
4196bc1bd01Sksagiyam   Level: intermediate
4206bc1bd01Sksagiyam 
4216bc1bd01Sksagiyam .seealso: `DMPlexReorderSetDefault()`
4226bc1bd01Sksagiyam @*/
423d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexReorderGetDefault(DM dm, DMPlexReorderDefaultFlag *reorder)
424d71ae5a4SJacob Faibussowitsch {
4256bc1bd01Sksagiyam   PetscFunctionBegin;
4266bc1bd01Sksagiyam   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
4274f572ea9SToby Isaac   PetscAssertPointer(reorder, 2);
4286bc1bd01Sksagiyam   PetscUseMethod(dm, "DMPlexReorderGetDefault_C", (DM, DMPlexReorderDefaultFlag *), (dm, reorder));
4293ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
4306bc1bd01Sksagiyam }
431*d02c7345SMatthew G. Knepley 
432*d02c7345SMatthew G. Knepley // Reorder to group split nodes
433*d02c7345SMatthew G. Knepley PetscErrorCode DMPlexCreateSectionPermutation_Internal(DM dm, IS *permutation, PetscBT *blockStarts)
434*d02c7345SMatthew G. Knepley {
435*d02c7345SMatthew G. Knepley   IS        permIS;
436*d02c7345SMatthew G. Knepley   PetscBT   bt, blst;
437*d02c7345SMatthew G. Knepley   PetscInt *perm;
438*d02c7345SMatthew G. Knepley   PetscInt  pStart, pEnd, i = 0;
439*d02c7345SMatthew G. Knepley 
440*d02c7345SMatthew G. Knepley   PetscFunctionBegin;
441*d02c7345SMatthew G. Knepley   *permutation = NULL;
442*d02c7345SMatthew G. Knepley   *blockStarts = NULL;
443*d02c7345SMatthew G. Knepley   {
444*d02c7345SMatthew G. Knepley     DMPlexReorderDefaultFlag reorder;
445*d02c7345SMatthew G. Knepley     PetscCall(DMPlexReorderSectionGetDefault(dm, &reorder));
446*d02c7345SMatthew G. Knepley     if (reorder != DMPLEX_REORDER_DEFAULT_TRUE) PetscFunctionReturn(PETSC_SUCCESS);
447*d02c7345SMatthew G. Knepley   }
448*d02c7345SMatthew G. Knepley   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
449*d02c7345SMatthew G. Knepley   PetscCall(PetscMalloc1(pEnd - pStart, &perm));
450*d02c7345SMatthew G. Knepley   PetscCall(PetscBTCreate(pEnd - pStart, &bt));
451*d02c7345SMatthew G. Knepley   PetscCall(PetscBTCreate(pEnd - pStart, &blst));
452*d02c7345SMatthew G. Knepley   for (PetscInt p = pStart; p < pEnd; ++p) {
453*d02c7345SMatthew G. Knepley     const PetscInt *supp, *cone;
454*d02c7345SMatthew G. Knepley     PetscInt        suppSize;
455*d02c7345SMatthew G. Knepley 
456*d02c7345SMatthew G. Knepley     if (PetscBTLookupSet(bt, p)) continue;
457*d02c7345SMatthew G. Knepley     PetscCall(PetscBTSet(blst, p));
458*d02c7345SMatthew G. Knepley     perm[i++] = p;
459*d02c7345SMatthew G. Knepley     // Check for tensor cells in the support
460*d02c7345SMatthew G. Knepley     PetscCall(DMPlexGetSupport(dm, p, &supp));
461*d02c7345SMatthew G. Knepley     PetscCall(DMPlexGetSupportSize(dm, p, &suppSize));
462*d02c7345SMatthew G. Knepley     for (PetscInt s = 0; s < suppSize; ++s) {
463*d02c7345SMatthew G. Knepley       DMPolytopeType ct;
464*d02c7345SMatthew G. Knepley       PetscInt       q, qq;
465*d02c7345SMatthew G. Knepley 
466*d02c7345SMatthew G. Knepley       PetscCall(DMPlexGetCellType(dm, supp[s], &ct));
467*d02c7345SMatthew G. Knepley       switch (ct) {
468*d02c7345SMatthew G. Knepley       case DM_POLYTOPE_POINT_PRISM_TENSOR:
469*d02c7345SMatthew G. Knepley       case DM_POLYTOPE_SEG_PRISM_TENSOR:
470*d02c7345SMatthew G. Knepley       case DM_POLYTOPE_TRI_PRISM_TENSOR:
471*d02c7345SMatthew G. Knepley       case DM_POLYTOPE_QUAD_PRISM_TENSOR:
472*d02c7345SMatthew G. Knepley         // If found, move up the split partner of the tensor cell, and the cell itself
473*d02c7345SMatthew G. Knepley         PetscCall(DMPlexGetCone(dm, supp[s], &cone));
474*d02c7345SMatthew G. Knepley         qq = supp[s];
475*d02c7345SMatthew G. Knepley         q  = (cone[0] == p) ? cone[1] : cone[0];
476*d02c7345SMatthew G. Knepley         if (!PetscBTLookupSet(bt, q)) {
477*d02c7345SMatthew G. Knepley           perm[i++] = q;
478*d02c7345SMatthew G. Knepley           s         = suppSize;
479*d02c7345SMatthew G. Knepley         }
480*d02c7345SMatthew G. Knepley         if (!PetscBTLookupSet(bt, qq)) {
481*d02c7345SMatthew G. Knepley           perm[i++] = qq;
482*d02c7345SMatthew G. Knepley           s         = suppSize;
483*d02c7345SMatthew G. Knepley         }
484*d02c7345SMatthew G. Knepley         break;
485*d02c7345SMatthew G. Knepley       default:
486*d02c7345SMatthew G. Knepley         break;
487*d02c7345SMatthew G. Knepley       }
488*d02c7345SMatthew G. Knepley     }
489*d02c7345SMatthew G. Knepley   }
490*d02c7345SMatthew G. Knepley   PetscCall(PetscBTDestroy(&bt));
491*d02c7345SMatthew G. Knepley   PetscCheck(i == pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of points in permutation %" PetscInt_FMT " does not match chart size %" PetscInt_FMT, i, pEnd - pStart);
492*d02c7345SMatthew G. Knepley   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS));
493*d02c7345SMatthew G. Knepley   PetscCall(ISSetPermutation(permIS));
494*d02c7345SMatthew G. Knepley   *permutation = permIS;
495*d02c7345SMatthew G. Knepley   *blockStarts = blst;
496*d02c7345SMatthew G. Knepley   PetscFunctionReturn(PETSC_SUCCESS);
497*d02c7345SMatthew G. Knepley }
498*d02c7345SMatthew G. Knepley 
499*d02c7345SMatthew G. Knepley PetscErrorCode DMPlexReorderSectionSetDefault_Plex(DM dm, DMPlexReorderDefaultFlag reorder)
500*d02c7345SMatthew G. Knepley {
501*d02c7345SMatthew G. Knepley   DM_Plex *mesh = (DM_Plex *)dm->data;
502*d02c7345SMatthew G. Knepley 
503*d02c7345SMatthew G. Knepley   PetscFunctionBegin;
504*d02c7345SMatthew G. Knepley   mesh->reorderSection = reorder;
505*d02c7345SMatthew G. Knepley   PetscFunctionReturn(PETSC_SUCCESS);
506*d02c7345SMatthew G. Knepley }
507*d02c7345SMatthew G. Knepley 
508*d02c7345SMatthew G. Knepley /*@
509*d02c7345SMatthew G. Knepley   DMPlexReorderSectionSetDefault - Set flag indicating whether the local section should be reordered by default
510*d02c7345SMatthew G. Knepley 
511*d02c7345SMatthew G. Knepley   Logically collective
512*d02c7345SMatthew G. Knepley 
513*d02c7345SMatthew G. Knepley   Input Parameters:
514*d02c7345SMatthew G. Knepley + dm      - The DM
515*d02c7345SMatthew G. Knepley - reorder - Flag for reordering
516*d02c7345SMatthew G. Knepley 
517*d02c7345SMatthew G. Knepley   Level: intermediate
518*d02c7345SMatthew G. Knepley 
519*d02c7345SMatthew G. Knepley .seealso: `DMPlexReorderSectionGetDefault()`
520*d02c7345SMatthew G. Knepley @*/
521*d02c7345SMatthew G. Knepley PetscErrorCode DMPlexReorderSectionSetDefault(DM dm, DMPlexReorderDefaultFlag reorder)
522*d02c7345SMatthew G. Knepley {
523*d02c7345SMatthew G. Knepley   PetscFunctionBegin;
524*d02c7345SMatthew G. Knepley   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
525*d02c7345SMatthew G. Knepley   PetscTryMethod(dm, "DMPlexReorderSectionSetDefault_C", (DM, DMPlexReorderDefaultFlag), (dm, reorder));
526*d02c7345SMatthew G. Knepley   PetscFunctionReturn(PETSC_SUCCESS);
527*d02c7345SMatthew G. Knepley }
528*d02c7345SMatthew G. Knepley 
529*d02c7345SMatthew G. Knepley PetscErrorCode DMPlexReorderSectionGetDefault_Plex(DM dm, DMPlexReorderDefaultFlag *reorder)
530*d02c7345SMatthew G. Knepley {
531*d02c7345SMatthew G. Knepley   DM_Plex *mesh = (DM_Plex *)dm->data;
532*d02c7345SMatthew G. Knepley 
533*d02c7345SMatthew G. Knepley   PetscFunctionBegin;
534*d02c7345SMatthew G. Knepley   *reorder = mesh->reorderSection;
535*d02c7345SMatthew G. Knepley   PetscFunctionReturn(PETSC_SUCCESS);
536*d02c7345SMatthew G. Knepley }
537*d02c7345SMatthew G. Knepley 
538*d02c7345SMatthew G. Knepley /*@
539*d02c7345SMatthew G. Knepley   DMPlexReorderSectionGetDefault - Get flag indicating whether the local section should be reordered by default
540*d02c7345SMatthew G. Knepley 
541*d02c7345SMatthew G. Knepley   Not collective
542*d02c7345SMatthew G. Knepley 
543*d02c7345SMatthew G. Knepley   Input Parameter:
544*d02c7345SMatthew G. Knepley . dm - The DM
545*d02c7345SMatthew G. Knepley 
546*d02c7345SMatthew G. Knepley   Output Parameter:
547*d02c7345SMatthew G. Knepley . reorder - Flag for reordering
548*d02c7345SMatthew G. Knepley 
549*d02c7345SMatthew G. Knepley   Level: intermediate
550*d02c7345SMatthew G. Knepley 
551*d02c7345SMatthew G. Knepley .seealso: `DMPlexReorderSetDefault()`
552*d02c7345SMatthew G. Knepley @*/
553*d02c7345SMatthew G. Knepley PetscErrorCode DMPlexReorderSectionGetDefault(DM dm, DMPlexReorderDefaultFlag *reorder)
554*d02c7345SMatthew G. Knepley {
555*d02c7345SMatthew G. Knepley   PetscFunctionBegin;
556*d02c7345SMatthew G. Knepley   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
557*d02c7345SMatthew G. Knepley   PetscAssertPointer(reorder, 2);
558*d02c7345SMatthew G. Knepley   PetscUseMethod(dm, "DMPlexReorderSectionGetDefault_C", (DM, DMPlexReorderDefaultFlag *), (dm, reorder));
559*d02c7345SMatthew G. Knepley   PetscFunctionReturn(PETSC_SUCCESS);
560*d02c7345SMatthew G. Knepley }
561