xref: /petsc/include/petscdm.h (revision d410b0cf18e1798d3d4c14858e0c2ffdbe2fea69)
1 /*
2       Objects to manage the interactions between the mesh data structures and the algebraic objects
3 */
4 #if !defined(PETSCDM_H)
5 #define PETSCDM_H
6 #include <petscmat.h>
7 #include <petscdmtypes.h>
8 #include <petscfetypes.h>
9 #include <petscdstypes.h>
10 #include <petscdmlabel.h>
11 
12 PETSC_EXTERN PetscErrorCode DMInitializePackage(void);
13 
14 PETSC_EXTERN PetscClassId DM_CLASSID;
15 PETSC_EXTERN PetscClassId DMLABEL_CLASSID;
16 
17 #define DMLOCATEPOINT_POINT_NOT_FOUND -367
18 
19 /*J
20     DMType - String with the name of a PETSc DM
21 
22    Level: beginner
23 
24 .seealso: DMSetType(), DM
25 J*/
26 typedef const char* DMType;
27 #define DMDA        "da"
28 #define DMCOMPOSITE "composite"
29 #define DMSLICED    "sliced"
30 #define DMSHELL     "shell"
31 #define DMPLEX      "plex"
32 #define DMREDUNDANT "redundant"
33 #define DMPATCH     "patch"
34 #define DMMOAB      "moab"
35 #define DMNETWORK   "network"
36 #define DMFOREST    "forest"
37 #define DMP4EST     "p4est"
38 #define DMP8EST     "p8est"
39 #define DMSWARM     "swarm"
40 #define DMPRODUCT   "product"
41 #define DMSTAG      "stag"
42 
43 PETSC_EXTERN const char *const DMBoundaryTypes[];
44 PETSC_EXTERN const char *const DMBoundaryConditionTypes[];
45 PETSC_EXTERN PetscFunctionList DMList;
46 PETSC_EXTERN PetscErrorCode DMCreate(MPI_Comm,DM*);
47 PETSC_EXTERN PetscErrorCode DMClone(DM,DM*);
48 PETSC_EXTERN PetscErrorCode DMSetType(DM, DMType);
49 PETSC_EXTERN PetscErrorCode DMGetType(DM, DMType *);
50 PETSC_EXTERN PetscErrorCode DMRegister(const char[],PetscErrorCode (*)(DM));
51 PETSC_EXTERN PetscErrorCode DMRegisterDestroy(void);
52 
53 PETSC_EXTERN PetscErrorCode DMView(DM,PetscViewer);
54 PETSC_EXTERN PetscErrorCode DMLoad(DM,PetscViewer);
55 PETSC_EXTERN PetscErrorCode DMDestroy(DM*);
56 PETSC_EXTERN PetscErrorCode DMCreateGlobalVector(DM,Vec*);
57 PETSC_EXTERN PetscErrorCode DMCreateLocalVector(DM,Vec*);
58 PETSC_EXTERN PetscErrorCode DMGetLocalVector(DM,Vec *);
59 PETSC_EXTERN PetscErrorCode DMRestoreLocalVector(DM,Vec *);
60 PETSC_EXTERN PetscErrorCode DMGetGlobalVector(DM,Vec *);
61 PETSC_EXTERN PetscErrorCode DMRestoreGlobalVector(DM,Vec *);
62 PETSC_EXTERN PetscErrorCode DMClearGlobalVectors(DM);
63 PETSC_EXTERN PetscErrorCode DMClearLocalVectors(DM);
64 PETSC_EXTERN PetscErrorCode DMHasNamedGlobalVector(DM,const char*,PetscBool*);
65 PETSC_EXTERN PetscErrorCode DMGetNamedGlobalVector(DM,const char*,Vec*);
66 PETSC_EXTERN PetscErrorCode DMRestoreNamedGlobalVector(DM,const char*,Vec*);
67 PETSC_EXTERN PetscErrorCode DMHasNamedLocalVector(DM,const char*,PetscBool*);
68 PETSC_EXTERN PetscErrorCode DMGetNamedLocalVector(DM,const char*,Vec*);
69 PETSC_EXTERN PetscErrorCode DMRestoreNamedLocalVector(DM,const char*,Vec*);
70 PETSC_EXTERN PetscErrorCode DMGetLocalToGlobalMapping(DM,ISLocalToGlobalMapping*);
71 PETSC_EXTERN PetscErrorCode DMCreateFieldIS(DM,PetscInt*,char***,IS**);
72 PETSC_EXTERN PetscErrorCode DMGetBlockSize(DM,PetscInt*);
73 PETSC_EXTERN PetscErrorCode DMCreateColoring(DM,ISColoringType,ISColoring*);
74 PETSC_EXTERN PetscErrorCode DMCreateMatrix(DM,Mat*);
75 PETSC_EXTERN PetscErrorCode DMSetMatrixPreallocateOnly(DM,PetscBool);
76 PETSC_EXTERN PetscErrorCode DMSetMatrixStructureOnly(DM,PetscBool);
77 PETSC_EXTERN PetscErrorCode DMCreateInterpolation(DM,DM,Mat*,Vec*);
78 PETSC_EXTERN PetscErrorCode DMCreateRestriction(DM,DM,Mat*);
79 PETSC_EXTERN PetscErrorCode DMCreateInjection(DM,DM,Mat*);
80 PETSC_EXTERN PetscErrorCode DMCreateMassMatrix(DM,DM,Mat*);
81 PETSC_EXTERN PetscErrorCode DMGetWorkArray(DM,PetscInt,MPI_Datatype,void*);
82 PETSC_EXTERN PetscErrorCode DMRestoreWorkArray(DM,PetscInt,MPI_Datatype,void*);
83 PETSC_EXTERN PetscErrorCode DMRefine(DM,MPI_Comm,DM*);
84 PETSC_EXTERN PetscErrorCode DMCoarsen(DM,MPI_Comm,DM*);
85 PETSC_EXTERN PetscErrorCode DMGetCoarseDM(DM,DM*);
86 PETSC_EXTERN PetscErrorCode DMSetCoarseDM(DM,DM);
87 PETSC_EXTERN PetscErrorCode DMGetFineDM(DM,DM*);
88 PETSC_EXTERN PetscErrorCode DMSetFineDM(DM,DM);
89 PETSC_EXTERN PetscErrorCode DMRefineHierarchy(DM,PetscInt,DM[]);
90 PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy(DM,PetscInt,DM[]);
91 PETSC_EXTERN PetscErrorCode DMCoarsenHookAdd(DM,PetscErrorCode (*)(DM,DM,void*),PetscErrorCode (*)(DM,Mat,Vec,Mat,DM,void*),void*);
92 PETSC_EXTERN PetscErrorCode DMCoarsenHookRemove(DM,PetscErrorCode (*)(DM,DM,void*),PetscErrorCode (*)(DM,Mat,Vec,Mat,DM,void*),void*);
93 PETSC_EXTERN PetscErrorCode DMRefineHookAdd(DM,PetscErrorCode (*)(DM,DM,void*),PetscErrorCode (*)(DM,Mat,DM,void*),void*);
94 PETSC_EXTERN PetscErrorCode DMRefineHookRemove(DM,PetscErrorCode (*)(DM,DM,void*),PetscErrorCode (*)(DM,Mat,DM,void*),void*);
95 PETSC_EXTERN PetscErrorCode DMRestrict(DM,Mat,Vec,Mat,DM);
96 PETSC_EXTERN PetscErrorCode DMInterpolate(DM,Mat,DM);
97 PETSC_EXTERN PetscErrorCode DMInterpolateSolution(DM,DM,Mat,Vec,Vec);
98 PETSC_EXTERN PetscErrorCode DMExtrude(DM,PetscInt,DM*);
99 PETSC_EXTERN PetscErrorCode DMSetFromOptions(DM);
100 PETSC_EXTERN PetscErrorCode DMViewFromOptions(DM,PetscObject,const char[]);
101 
102 PETSC_EXTERN PetscErrorCode DMAdaptLabel(DM,DMLabel,DM*);
103 PETSC_EXTERN PetscErrorCode DMAdaptMetric(DM, Vec, DMLabel, DM *);
104 
105 PETSC_EXTERN PetscErrorCode DMSetUp(DM);
106 PETSC_EXTERN PetscErrorCode DMCreateInterpolationScale(DM,DM,Mat,Vec*);
107 PETSC_EXTERN PETSC_DEPRECATED_FUNCTION("Use DMDACreateAggregates() or DMCreateRestriction() (since version 3.12)") PetscErrorCode DMCreateAggregates(DM,DM,Mat*);
108 PETSC_EXTERN PetscErrorCode DMGlobalToLocalHookAdd(DM,PetscErrorCode (*)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*)(DM,Vec,InsertMode,Vec,void*),void*);
109 PETSC_EXTERN PetscErrorCode DMLocalToGlobalHookAdd(DM,PetscErrorCode (*)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*)(DM,Vec,InsertMode,Vec,void*),void*);
110 PETSC_EXTERN PetscErrorCode DMGlobalToLocal(DM,Vec,InsertMode,Vec);
111 PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin(DM,Vec,InsertMode,Vec);
112 PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd(DM,Vec,InsertMode,Vec);
113 PETSC_EXTERN PetscErrorCode DMLocalToGlobal(DM,Vec,InsertMode,Vec);
114 PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin(DM,Vec,InsertMode,Vec);
115 PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd(DM,Vec,InsertMode,Vec);
116 PETSC_EXTERN PetscErrorCode DMLocalToLocalBegin(DM,Vec,InsertMode,Vec);
117 PETSC_EXTERN PetscErrorCode DMLocalToLocalEnd(DM,Vec,InsertMode,Vec);
118 PETSC_EXTERN PetscErrorCode DMConvert(DM,DMType,DM*);
119 
120 /* Topology support */
121 PETSC_EXTERN PetscErrorCode DMGetDimension(DM,PetscInt*);
122 PETSC_EXTERN PetscErrorCode DMSetDimension(DM,PetscInt);
123 PETSC_EXTERN PetscErrorCode DMGetDimPoints(DM,PetscInt,PetscInt*,PetscInt*);
124 PETSC_EXTERN PetscErrorCode DMGetUseNatural(DM,PetscBool*);
125 PETSC_EXTERN PetscErrorCode DMSetUseNatural(DM,PetscBool);
126 
127 /* Coordinate support */
128 PETSC_EXTERN PetscErrorCode DMGetCoordinateDM(DM,DM*);
129 PETSC_EXTERN PetscErrorCode DMSetCoordinateDM(DM,DM);
130 PETSC_EXTERN PetscErrorCode DMGetCoordinateDim(DM,PetscInt*);
131 PETSC_EXTERN PetscErrorCode DMSetCoordinateDim(DM,PetscInt);
132 PETSC_EXTERN PetscErrorCode DMGetCoordinateSection(DM,PetscSection*);
133 PETSC_EXTERN PetscErrorCode DMSetCoordinateSection(DM,PetscInt,PetscSection);
134 PETSC_EXTERN PetscErrorCode DMGetCoordinates(DM,Vec*);
135 PETSC_EXTERN PetscErrorCode DMSetCoordinates(DM,Vec);
136 PETSC_EXTERN PetscErrorCode DMGetCoordinatesLocal(DM,Vec*);
137 PETSC_EXTERN PetscErrorCode DMGetCoordinatesLocalSetUp(DM);
138 PETSC_EXTERN PetscErrorCode DMGetCoordinatesLocalNoncollective(DM,Vec*);
139 PETSC_EXTERN PetscErrorCode DMGetCoordinatesLocalTuple(DM,IS,PetscSection*,Vec*);
140 PETSC_EXTERN PetscErrorCode DMSetCoordinatesLocal(DM,Vec);
141 PETSC_EXTERN PetscErrorCode DMLocatePoints(DM,Vec,DMPointLocationType,PetscSF*);
142 PETSC_EXTERN PetscErrorCode DMGetPeriodicity(DM,PetscBool*,const PetscReal**,const PetscReal**,const DMBoundaryType**);
143 PETSC_EXTERN PetscErrorCode DMSetPeriodicity(DM,PetscBool,const PetscReal[],const PetscReal[],const DMBoundaryType[]);
144 PETSC_EXTERN PetscErrorCode DMLocalizeCoordinate(DM, const PetscScalar[], PetscBool, PetscScalar[]);
145 PETSC_EXTERN PetscErrorCode DMLocalizeCoordinates(DM);
146 PETSC_EXTERN PetscErrorCode DMGetCoordinatesLocalized(DM,PetscBool*);
147 PETSC_EXTERN PetscErrorCode DMGetCoordinatesLocalizedLocal(DM,PetscBool*);
148 PETSC_EXTERN PetscErrorCode DMGetNeighbors(DM,PetscInt*,const PetscMPIInt**);
149 PETSC_EXTERN PetscErrorCode DMGetCoordinateField(DM,DMField*);
150 PETSC_EXTERN PetscErrorCode DMSetCoordinateField(DM,DMField);
151 PETSC_EXTERN PetscErrorCode DMGetBoundingBox(DM,PetscReal[],PetscReal[]);
152 PETSC_EXTERN PetscErrorCode DMGetLocalBoundingBox(DM,PetscReal[],PetscReal[]);
153 PETSC_EXTERN PetscErrorCode DMProjectCoordinates(DM,PetscFE);
154 
155 /* block hook interface */
156 PETSC_EXTERN PetscErrorCode DMSubDomainHookAdd(DM,PetscErrorCode (*)(DM,DM,void*),PetscErrorCode (*)(DM,VecScatter,VecScatter,DM,void*),void*);
157 PETSC_EXTERN PetscErrorCode DMSubDomainHookRemove(DM,PetscErrorCode (*)(DM,DM,void*),PetscErrorCode (*)(DM,VecScatter,VecScatter,DM,void*),void*);
158 PETSC_EXTERN PetscErrorCode DMSubDomainRestrict(DM,VecScatter,VecScatter,DM);
159 
160 PETSC_EXTERN PetscErrorCode DMSetOptionsPrefix(DM,const char []);
161 PETSC_EXTERN PetscErrorCode DMAppendOptionsPrefix(DM,const char []);
162 PETSC_EXTERN PetscErrorCode DMGetOptionsPrefix(DM,const char*[]);
163 PETSC_EXTERN PetscErrorCode DMSetVecType(DM,VecType);
164 PETSC_EXTERN PetscErrorCode DMGetVecType(DM,VecType*);
165 PETSC_EXTERN PetscErrorCode DMSetMatType(DM,MatType);
166 PETSC_EXTERN PetscErrorCode DMGetMatType(DM,MatType*);
167 PETSC_EXTERN PetscErrorCode DMSetISColoringType(DM,ISColoringType);
168 PETSC_EXTERN PetscErrorCode DMGetISColoringType(DM,ISColoringType*);
169 PETSC_EXTERN PetscErrorCode DMSetApplicationContext(DM,void*);
170 PETSC_EXTERN PetscErrorCode DMSetApplicationContextDestroy(DM,PetscErrorCode (*)(void**));
171 PETSC_EXTERN PetscErrorCode DMGetApplicationContext(DM,void*);
172 PETSC_EXTERN PetscErrorCode DMSetVariableBounds(DM,PetscErrorCode (*)(DM,Vec,Vec));
173 PETSC_EXTERN PetscErrorCode DMHasVariableBounds(DM,PetscBool *);
174 PETSC_EXTERN PetscErrorCode DMHasColoring(DM,PetscBool *);
175 PETSC_EXTERN PetscErrorCode DMHasCreateRestriction(DM,PetscBool *);
176 PETSC_EXTERN PetscErrorCode DMHasCreateInjection(DM,PetscBool *);
177 PETSC_EXTERN PetscErrorCode DMComputeVariableBounds(DM,Vec,Vec);
178 
179 PETSC_EXTERN PetscErrorCode DMCreateSubDM(DM, PetscInt, const PetscInt[], IS *, DM *);
180 PETSC_EXTERN PetscErrorCode DMCreateSuperDM(DM[], PetscInt, IS **, DM *);
181 PETSC_EXTERN PetscErrorCode DMCreateSectionSubDM(DM,PetscInt,const PetscInt[],IS*,DM*);
182 PETSC_EXTERN PetscErrorCode DMCreateSectionSuperDM(DM[],PetscInt,IS**,DM*);
183 PETSC_EXTERN PetscErrorCode DMCreateFieldDecomposition(DM,PetscInt*,char***,IS**,DM**);
184 PETSC_EXTERN PetscErrorCode DMCreateDomainDecomposition(DM,PetscInt*,char***,IS**,IS**,DM**);
185 PETSC_EXTERN PetscErrorCode DMCreateDomainDecompositionScatters(DM,PetscInt,DM*,VecScatter**,VecScatter**,VecScatter**);
186 
187 PETSC_EXTERN PetscErrorCode DMGetRefineLevel(DM,PetscInt*);
188 PETSC_EXTERN PetscErrorCode DMSetRefineLevel(DM,PetscInt);
189 PETSC_EXTERN PetscErrorCode DMGetCoarsenLevel(DM,PetscInt*);
190 PETSC_EXTERN PetscErrorCode DMSetCoarsenLevel(DM,PetscInt);
191 PETSC_EXTERN PetscErrorCode DMFinalizePackage(void);
192 
193 PETSC_EXTERN PetscErrorCode VecGetDM(Vec, DM*);
194 PETSC_EXTERN PetscErrorCode VecSetDM(Vec, DM);
195 PETSC_EXTERN PetscErrorCode MatGetDM(Mat, DM*);
196 PETSC_EXTERN PetscErrorCode MatSetDM(Mat, DM);
197 PETSC_EXTERN PetscErrorCode MatFDColoringUseDM(Mat,MatFDColoring);
198 
199 typedef struct NLF_DAAD* NLF;
200 
201 #define DM_FILE_CLASSID 1211221
202 
203 /* FEM support */
204 PETSC_EXTERN PetscErrorCode DMPrintCellVector(PetscInt, const char [], PetscInt, const PetscScalar []);
205 PETSC_EXTERN PetscErrorCode DMPrintCellMatrix(PetscInt, const char [], PetscInt, PetscInt, const PetscScalar []);
206 PETSC_EXTERN PetscErrorCode DMPrintLocalVec(DM, const char [], PetscReal, Vec);
207 
208 PETSC_EXTERN PetscErrorCode DMSetNullSpaceConstructor(DM, PetscInt, PetscErrorCode (*)(DM, PetscInt, PetscInt, MatNullSpace *));
209 PETSC_EXTERN PetscErrorCode DMGetNullSpaceConstructor(DM, PetscInt, PetscErrorCode (**)(DM, PetscInt, PetscInt, MatNullSpace *));
210 PETSC_EXTERN PetscErrorCode DMSetNearNullSpaceConstructor(DM, PetscInt, PetscErrorCode (*)(DM, PetscInt, PetscInt, MatNullSpace *));
211 PETSC_EXTERN PetscErrorCode DMGetNearNullSpaceConstructor(DM, PetscInt, PetscErrorCode (**)(DM, PetscInt, PetscInt, MatNullSpace *));
212 
213 PETSC_EXTERN PetscErrorCode DMGetSection(DM, PetscSection *); /* Use DMGetLocalSection() in new code (since v3.12) */
214 PETSC_EXTERN PetscErrorCode DMSetSection(DM, PetscSection);   /* Use DMSetLocalSection() in new code (since v3.12) */
215 PETSC_EXTERN PetscErrorCode DMGetLocalSection(DM, PetscSection *);
216 PETSC_EXTERN PetscErrorCode DMSetLocalSection(DM, PetscSection);
217 PETSC_EXTERN PetscErrorCode DMGetGlobalSection(DM, PetscSection *);
218 PETSC_EXTERN PetscErrorCode DMSetGlobalSection(DM, PetscSection);
219 PETSC_STATIC_INLINE PETSC_DEPRECATED_FUNCTION("Use DMGetSection() (since v3.9)") PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *s) {return DMGetSection(dm,s);}
220 PETSC_STATIC_INLINE PETSC_DEPRECATED_FUNCTION("Use DMSetSection() (since v3.9)") PetscErrorCode DMSetDefaultSection(DM dm, PetscSection s) {return DMSetSection(dm,s);}
221 PETSC_STATIC_INLINE PETSC_DEPRECATED_FUNCTION("Use DMGetGlobalSection() (since v3.9)") PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *s) {return DMGetGlobalSection(dm,s);}
222 PETSC_STATIC_INLINE PETSC_DEPRECATED_FUNCTION("Use DMSetGlobalSection() (since v3.9)") PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection s) {return DMSetGlobalSection(dm,s);}
223 
224 PETSC_EXTERN PetscErrorCode DMGetSectionSF(DM, PetscSF*);
225 PETSC_EXTERN PetscErrorCode DMSetSectionSF(DM, PetscSF);
226 PETSC_EXTERN PetscErrorCode DMCreateSectionSF(DM, PetscSection, PetscSection);
227 PETSC_STATIC_INLINE PETSC_DEPRECATED_FUNCTION("Use DMGetSectionSF() (since v3.12)") PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *s) {return DMGetSectionSF(dm,s);}
228 PETSC_STATIC_INLINE PETSC_DEPRECATED_FUNCTION("Use DMSetSectionSF() (since v3.12)") PetscErrorCode DMSetDefaultSF(DM dm, PetscSF s) {return DMSetSectionSF(dm,s);}
229 PETSC_STATIC_INLINE PETSC_DEPRECATED_FUNCTION("Use DMCreateSectionSF() (since v3.12)") PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection l, PetscSection g) {return DMCreateSectionSF(dm,l,g);}
230 PETSC_EXTERN PetscErrorCode DMGetPointSF(DM, PetscSF *);
231 PETSC_EXTERN PetscErrorCode DMSetPointSF(DM, PetscSF);
232 
233 PETSC_EXTERN PetscErrorCode DMGetDefaultConstraints(DM, PetscSection *, Mat *);
234 PETSC_EXTERN PetscErrorCode DMSetDefaultConstraints(DM, PetscSection, Mat);
235 
236 PETSC_EXTERN PetscErrorCode DMGetOutputDM(DM, DM *);
237 PETSC_EXTERN PetscErrorCode DMGetOutputSequenceNumber(DM, PetscInt *, PetscReal *);
238 PETSC_EXTERN PetscErrorCode DMSetOutputSequenceNumber(DM, PetscInt, PetscReal);
239 PETSC_EXTERN PetscErrorCode DMOutputSequenceLoad(DM, PetscViewer, const char *, PetscInt, PetscReal *);
240 
241 PETSC_EXTERN PetscErrorCode DMGetNumFields(DM, PetscInt *);
242 PETSC_EXTERN PetscErrorCode DMSetNumFields(DM, PetscInt);
243 PETSC_EXTERN PetscErrorCode DMGetField(DM, PetscInt, DMLabel *, PetscObject *);
244 PETSC_EXTERN PetscErrorCode DMSetField(DM, PetscInt, DMLabel, PetscObject);
245 PETSC_EXTERN PetscErrorCode DMAddField(DM, DMLabel, PetscObject);
246 PETSC_EXTERN PetscErrorCode DMSetFieldAvoidTensor(DM, PetscInt, PetscBool);
247 PETSC_EXTERN PetscErrorCode DMGetFieldAvoidTensor(DM, PetscInt, PetscBool *);
248 PETSC_EXTERN PetscErrorCode DMClearFields(DM);
249 PETSC_EXTERN PetscErrorCode DMCopyFields(DM, DM);
250 PETSC_EXTERN PetscErrorCode DMGetAdjacency(DM, PetscInt, PetscBool *, PetscBool *);
251 PETSC_EXTERN PetscErrorCode DMSetAdjacency(DM, PetscInt, PetscBool, PetscBool);
252 PETSC_EXTERN PetscErrorCode DMGetBasicAdjacency(DM, PetscBool *, PetscBool *);
253 PETSC_EXTERN PetscErrorCode DMSetBasicAdjacency(DM, PetscBool, PetscBool);
254 
255 PETSC_EXTERN PetscErrorCode DMGetNumDS(DM, PetscInt *);
256 PETSC_EXTERN PetscErrorCode DMGetDS(DM, PetscDS *);
257 PETSC_EXTERN PetscErrorCode DMGetCellDS(DM, PetscInt, PetscDS *);
258 PETSC_EXTERN PetscErrorCode DMGetRegionDS(DM, DMLabel, IS *, PetscDS *);
259 PETSC_EXTERN PetscErrorCode DMSetRegionDS(DM, DMLabel, IS, PetscDS);
260 PETSC_EXTERN PetscErrorCode DMGetRegionNumDS(DM, PetscInt, DMLabel *, IS *, PetscDS *);
261 PETSC_EXTERN PetscErrorCode DMSetRegionNumDS(DM, PetscInt, DMLabel, IS, PetscDS);
262 PETSC_EXTERN PetscErrorCode DMFindRegionNum(DM, PetscDS, PetscInt *);
263 PETSC_EXTERN PetscErrorCode DMCreateDS(DM);
264 PETSC_EXTERN PetscErrorCode DMClearDS(DM);
265 PETSC_EXTERN PetscErrorCode DMCopyDS(DM, DM);
266 PETSC_EXTERN PetscErrorCode DMCopyDisc(DM, DM);
267 PETSC_EXTERN PetscErrorCode DMComputeExactSolution(DM, PetscReal, Vec, Vec);
268 PETSC_EXTERN PetscErrorCode DMGetNumAuxiliaryVec(DM, PetscInt *);
269 PETSC_EXTERN PetscErrorCode DMGetAuxiliaryVec(DM, DMLabel, PetscInt, Vec *);
270 PETSC_EXTERN PetscErrorCode DMSetAuxiliaryVec(DM, DMLabel, PetscInt, Vec);
271 PETSC_EXTERN PetscErrorCode DMGetAuxiliaryLabels(DM, DMLabel[], PetscInt[]);
272 PETSC_EXTERN PetscErrorCode DMCopyAuxiliaryVec(DM, DM);
273 
274 /*MC
275   DMInterpolationInfo - Structure for holding information about interpolation on a mesh
276 
277   Level: intermediate
278 
279   Synopsis:
280     comm   - The communicator
281     dim    - The spatial dimension of points
282     nInput - The number of input points
283     points - The input point coordinates
284     cells  - The cell containing each point
285     n      - The number of local points
286     coords - The point coordinates
287     dof    - The number of components to interpolate
288 
289 .seealso: DMInterpolationCreate(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
290 M*/
291 struct _DMInterpolationInfo {
292   MPI_Comm   comm;
293   PetscInt   dim;    /*1 The spatial dimension of points */
294   PetscInt   nInput; /* The number of input points */
295   PetscReal *points; /* The input point coordinates */
296   PetscInt  *cells;  /* The cell containing each point */
297   PetscInt   n;      /* The number of local points */
298   Vec        coords; /* The point coordinates */
299   PetscInt   dof;    /* The number of components to interpolate */
300 };
301 typedef struct _DMInterpolationInfo *DMInterpolationInfo;
302 
303 PETSC_EXTERN PetscErrorCode DMInterpolationCreate(MPI_Comm, DMInterpolationInfo *);
304 PETSC_EXTERN PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo, PetscInt);
305 PETSC_EXTERN PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo, PetscInt *);
306 PETSC_EXTERN PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo, PetscInt);
307 PETSC_EXTERN PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo, PetscInt *);
308 PETSC_EXTERN PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo, PetscInt, PetscReal[]);
309 PETSC_EXTERN PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo, DM, PetscBool, PetscBool);
310 PETSC_EXTERN PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo, Vec *);
311 PETSC_EXTERN PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo, Vec *);
312 PETSC_EXTERN PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo, Vec *);
313 PETSC_EXTERN PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo, DM, Vec, Vec);
314 PETSC_EXTERN PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *);
315 
316 PETSC_EXTERN PetscErrorCode DMCreateLabel(DM, const char []);
317 PETSC_EXTERN PetscErrorCode DMGetLabelValue(DM, const char[], PetscInt, PetscInt *);
318 PETSC_EXTERN PetscErrorCode DMSetLabelValue(DM, const char[], PetscInt, PetscInt);
319 PETSC_EXTERN PetscErrorCode DMClearLabelValue(DM, const char[], PetscInt, PetscInt);
320 PETSC_EXTERN PetscErrorCode DMGetLabelSize(DM, const char[], PetscInt *);
321 PETSC_EXTERN PetscErrorCode DMGetLabelIdIS(DM, const char[], IS *);
322 PETSC_EXTERN PetscErrorCode DMGetStratumSize(DM, const char [], PetscInt, PetscInt *);
323 PETSC_EXTERN PetscErrorCode DMGetStratumIS(DM, const char [], PetscInt, IS *);
324 PETSC_EXTERN PetscErrorCode DMSetStratumIS(DM, const char [], PetscInt, IS);
325 PETSC_EXTERN PetscErrorCode DMClearLabelStratum(DM, const char[], PetscInt);
326 PETSC_EXTERN PetscErrorCode DMGetLabelOutput(DM, const char[], PetscBool *);
327 PETSC_EXTERN PetscErrorCode DMSetLabelOutput(DM, const char[], PetscBool);
328 
329 PETSC_EXTERN PetscErrorCode DMGetNumLabels(DM, PetscInt *);
330 PETSC_EXTERN PetscErrorCode DMGetLabelName(DM, PetscInt, const char **);
331 PETSC_EXTERN PetscErrorCode DMHasLabel(DM, const char [], PetscBool *);
332 PETSC_EXTERN PetscErrorCode DMGetLabel(DM, const char *, DMLabel *);
333 PETSC_EXTERN PetscErrorCode DMSetLabel(DM, DMLabel);
334 PETSC_EXTERN PetscErrorCode DMGetLabelByNum(DM, PetscInt, DMLabel *);
335 PETSC_EXTERN PetscErrorCode DMAddLabel(DM, DMLabel);
336 PETSC_EXTERN PetscErrorCode DMRemoveLabel(DM, const char [], DMLabel *);
337 PETSC_EXTERN PetscErrorCode DMRemoveLabelBySelf(DM, DMLabel *, PetscBool);
338 PETSC_EXTERN PetscErrorCode DMCopyLabels(DM, DM, PetscCopyMode, PetscBool);
339 
340 PETSC_EXTERN PetscErrorCode DMAddBoundary(DM, DMBoundaryConditionType, const char[], DMLabel, PetscInt, const PetscInt[], PetscInt, PetscInt, const PetscInt[], void (*)(void), void (*)(void), void *, PetscInt *);
341 PETSC_EXTERN PetscErrorCode DMIsBoundaryPoint(DM, PetscInt, PetscBool *);
342 
343 PETSC_EXTERN PetscErrorCode DMProjectFunction(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void**,InsertMode,Vec);
344 PETSC_EXTERN PetscErrorCode DMProjectFunctionLocal(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void**,InsertMode,Vec);
345 PETSC_EXTERN PetscErrorCode DMProjectFunctionLabel(DM, PetscReal, DMLabel, PetscInt, const PetscInt[], PetscInt, const PetscInt[], PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **, InsertMode, Vec);
346 PETSC_EXTERN PetscErrorCode DMProjectFunctionLabelLocal(DM,PetscReal,DMLabel,PetscInt,const PetscInt[],PetscInt,const PetscInt[],PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
347 PETSC_EXTERN PetscErrorCode DMProjectFieldLocal(DM,PetscReal,Vec,void (**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],PetscInt,const PetscScalar[],PetscScalar[]),InsertMode,Vec);
348 PETSC_EXTERN PetscErrorCode DMProjectFieldLabelLocal(DM,PetscReal,DMLabel,PetscInt,const PetscInt[],PetscInt,const PetscInt[],Vec,void (**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],PetscInt,const PetscScalar[],PetscScalar[]),InsertMode,Vec);
349 PETSC_EXTERN PetscErrorCode DMProjectBdFieldLabelLocal(DM,PetscReal,DMLabel,PetscInt,const PetscInt[],PetscInt,const PetscInt[],Vec,void (**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],const PetscReal[],PetscInt,const PetscScalar[],PetscScalar[]),InsertMode,Vec);
350 PETSC_EXTERN PetscErrorCode DMComputeL2Diff(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
351 PETSC_EXTERN PetscErrorCode DMComputeL2GradientDiff(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **, Vec, const PetscReal [], PetscReal *);
352 PETSC_EXTERN PetscErrorCode DMComputeL2FieldDiff(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
353 PETSC_EXTERN PetscErrorCode DMComputeError(DM, Vec, PetscReal[], Vec *);
354 PETSC_EXTERN PetscErrorCode DMHasBasisTransform(DM,PetscBool*);
355 PETSC_EXTERN PetscErrorCode DMCopyTransform(DM, DM);
356 
357 PETSC_EXTERN PetscErrorCode DMGetCompatibility(DM,DM,PetscBool*,PetscBool*);
358 
359 PETSC_EXTERN PetscErrorCode DMMonitorSet(DM, PetscErrorCode (*)(DM, void *), void *, PetscErrorCode (*)(void**));
360 PETSC_EXTERN PetscErrorCode DMMonitorCancel(DM);
361 PETSC_EXTERN PetscErrorCode DMMonitorSetFromOptions(DM, const char[], const char[], const char[], PetscErrorCode (*)(DM, void *), PetscErrorCode (*)(DM, PetscViewerAndFormat *), PetscBool *);
362 PETSC_EXTERN PetscErrorCode DMMonitor(DM);
363 
364 PETSC_STATIC_INLINE PetscInt DMPolytopeTypeGetDim(DMPolytopeType ct)
365 {
366   switch (ct) {
367     case DM_POLYTOPE_POINT:
368       return 0;
369     case DM_POLYTOPE_SEGMENT:
370     case DM_POLYTOPE_POINT_PRISM_TENSOR:
371       return 1;
372     case DM_POLYTOPE_TRIANGLE:
373     case DM_POLYTOPE_QUADRILATERAL:
374     case DM_POLYTOPE_SEG_PRISM_TENSOR:
375       return 2;
376     case DM_POLYTOPE_TETRAHEDRON:
377     case DM_POLYTOPE_HEXAHEDRON:
378     case DM_POLYTOPE_TRI_PRISM:
379     case DM_POLYTOPE_TRI_PRISM_TENSOR:
380     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
381     case DM_POLYTOPE_PYRAMID:
382       return 3;
383     default: return -1;
384   }
385 }
386 
387 PETSC_STATIC_INLINE PetscInt DMPolytopeTypeGetConeSize(DMPolytopeType ct)
388 {
389   switch (ct) {
390     case DM_POLYTOPE_POINT:              return 0;
391     case DM_POLYTOPE_SEGMENT:            return 2;
392     case DM_POLYTOPE_POINT_PRISM_TENSOR: return 2;
393     case DM_POLYTOPE_TRIANGLE:           return 3;
394     case DM_POLYTOPE_QUADRILATERAL:      return 4;
395     case DM_POLYTOPE_SEG_PRISM_TENSOR:   return 4;
396     case DM_POLYTOPE_TETRAHEDRON:        return 4;
397     case DM_POLYTOPE_HEXAHEDRON:         return 6;
398     case DM_POLYTOPE_TRI_PRISM:          return 5;
399     case DM_POLYTOPE_TRI_PRISM_TENSOR:   return 5;
400     case DM_POLYTOPE_QUAD_PRISM_TENSOR:  return 6;
401     case DM_POLYTOPE_PYRAMID:            return 5;
402     default: return -1;
403   }
404 }
405 
406 PETSC_STATIC_INLINE PetscInt DMPolytopeTypeGetNumVertices(DMPolytopeType ct)
407 {
408   switch (ct) {
409     case DM_POLYTOPE_POINT:              return 1;
410     case DM_POLYTOPE_SEGMENT:            return 2;
411     case DM_POLYTOPE_POINT_PRISM_TENSOR: return 2;
412     case DM_POLYTOPE_TRIANGLE:           return 3;
413     case DM_POLYTOPE_QUADRILATERAL:      return 4;
414     case DM_POLYTOPE_SEG_PRISM_TENSOR:   return 4;
415     case DM_POLYTOPE_TETRAHEDRON:        return 4;
416     case DM_POLYTOPE_HEXAHEDRON:         return 8;
417     case DM_POLYTOPE_TRI_PRISM:          return 6;
418     case DM_POLYTOPE_TRI_PRISM_TENSOR:   return 6;
419     case DM_POLYTOPE_QUAD_PRISM_TENSOR:  return 8;
420     case DM_POLYTOPE_PYRAMID:            return 5;
421     default: return -1;
422   }
423 }
424 
425 PETSC_STATIC_INLINE DMPolytopeType DMPolytopeTypeSimpleShape(PetscInt dim, PetscBool simplex)
426 {
427   return dim == 0 ? DM_POLYTOPE_POINT :
428         (dim == 1 ? DM_POLYTOPE_SEGMENT :
429         (dim == 2 ? (simplex ? DM_POLYTOPE_TRIANGLE : DM_POLYTOPE_QUADRILATERAL) :
430         (dim == 3 ? (simplex ? DM_POLYTOPE_TETRAHEDRON : DM_POLYTOPE_HEXAHEDRON) : DM_POLYTOPE_UNKNOWN)));
431 }
432 
433 PETSC_STATIC_INLINE PetscInt DMPolytopeTypeGetNumArrangments(DMPolytopeType ct)
434 {
435   switch (ct) {
436     case DM_POLYTOPE_POINT:              return 1;
437     case DM_POLYTOPE_SEGMENT:            return 2;
438     case DM_POLYTOPE_POINT_PRISM_TENSOR: return 2;
439     case DM_POLYTOPE_TRIANGLE:           return 6;
440     case DM_POLYTOPE_QUADRILATERAL:      return 8;
441     case DM_POLYTOPE_SEG_PRISM_TENSOR:   return 4;
442     case DM_POLYTOPE_TETRAHEDRON:        return 24;
443     case DM_POLYTOPE_HEXAHEDRON:         return 48;
444     case DM_POLYTOPE_TRI_PRISM:          return 12;
445     case DM_POLYTOPE_TRI_PRISM_TENSOR:   return 12;
446     case DM_POLYTOPE_QUAD_PRISM_TENSOR:  return 16;
447     case DM_POLYTOPE_PYRAMID:            return 8;
448     default: return -1;
449   }
450 }
451 
452 /* An arrangement is a face order combined with an orientation for each face */
453 PETSC_STATIC_INLINE const PetscInt *DMPolytopeTypeGetArrangment(DMPolytopeType ct, PetscInt o)
454 {
455   static const PetscInt pntArr[1*2] = {0, 0};
456   /* a: swap */
457   static const PetscInt segArr[2*2*2] = {
458     1, 0,  0, 0, /* -1: a */
459     0, 0,  1, 0, /*  0: e */};
460   /* a: swap first two
461      b: swap last two */
462   static const PetscInt triArr[6*3*2] = {
463     0, -1,  2, -1,  1, -1, /* -3: b */
464     2, -1,  1, -1,  0, -1, /* -2: aba */
465     1, -1,  0, -1,  2, -1, /* -1: a */
466     0,  0,  1,  0,  2,  0, /*  0: identity */
467     1,  0,  2,  0,  0,  0, /*  1: ba */
468     2,  0,  0,  0,  1,  0, /*  2: ab */};
469   /* a: forward cyclic permutation
470      b: swap first and last pairs */
471   static const PetscInt quadArr[8*4*2] = {
472     1, -1,  0, -1,  3, -1,  2, -1, /* -4: b */
473     0, -1,  3, -1,  2, -1,  1, -1, /* -3: b a^3 = a b */
474     3, -1,  2, -1,  1, -1,  0, -1, /* -2: b a^2 = a^2 b */
475     2, -1,  1, -1,  0, -1,  3, -1, /* -1: b a   = a^3 b */
476     0,  0,  1,  0,  2,  0,  3,  0, /*  0: identity */
477     1,  0,  2,  0,  3,  0,  0,  0, /*  1: a */
478     2,  0,  3,  0,  0,  0,  1,  0, /*  2: a^2 */
479     3,  0,  0,  0,  1,  0,  2,  0, /*  3: a^3 */};
480   /* r: rotate 180
481      b: swap top and bottom segments */
482   static const PetscInt tsegArr[4*4*2] = {
483     1, -1,  0, -1,  3, -1,  2, -1, /* -2: r b */
484     0, -1,  1, -1,  3,  0,  2,  0, /* -1: r */
485     0,  0,  1,  0,  2,  0,  3,  0, /*  0: identity */
486     1,  0,  0,  0,  2, -1,  3, -1, /*  1: b */};
487   /* https://en.wikiversity.org/wiki/Symmetric_group_S4 */
488   static const PetscInt tetArr[24*4*2] = {
489     3, -2,  2, -3,  0, -1,  1, -1, /* -12: (1324)   p22 */
490     3, -1,  1, -3,  2, -1,  0, -1, /* -11: (14)     p21 */
491     3, -3,  0, -3,  1, -1,  2, -1, /* -10: (1234)   p18 */
492     2, -1,  3, -1,  1, -3,  0, -2, /*  -9: (1423)   p17 */
493     2, -3,  0, -1,  3, -2,  1, -3, /*  -8: (1342)   p13 */
494     2, -2,  1, -2,  0, -2,  3, -2, /*  -7: (24)     p14 */
495     1, -2,  0, -2,  2, -2,  3, -1, /*  -6: (34)     p6  */
496     1, -1,  3, -3,  0, -3,  2, -2, /*  -5: (1243)   p10 */
497     1, -3,  2, -1,  3, -1,  0, -3, /*  -4: (1432)   p9  */
498     0, -3,  1, -1,  3, -3,  2, -3, /*  -3: (12)     p1  */
499     0, -2,  2, -2,  1, -2,  3, -3, /*  -2: (23)     p2  */
500     0, -1,  3, -2,  2, -3,  1, -2, /*  -1: (13)     p5  */
501     0,  0,  1,  0,  2,  0,  3,  0, /*   0: ()       p0  */
502     0,  1,  3,  1,  1,  2,  2,  0, /*   1: (123)    p4  */
503     0,  2,  2,  1,  3,  0,  1,  2, /*   2: (132)    p3  */
504     1,  2,  0,  1,  3,  1,  2,  2, /*   3: (12)(34) p7  */
505     1,  0,  2,  0,  0,  0,  3,  1, /*   4: (243)    p8  */
506     1,  1,  3,  2,  2,  2,  0,  0, /*   5: (143)    p11 */
507     2,  1,  3,  0,  0,  2,  1,  0, /*   6: (13)(24) p16 */
508     2,  2,  1,  1,  3,  2,  0,  2, /*   7: (142)    p15 */
509     2,  0,  0,  0,  1,  0,  3,  2, /*   8: (234)    p12 */
510     3,  2,  2,  2,  1,  1,  0,  1, /*   9: (14)(23) p23 */
511     3,  0,  0,  2,  2,  1,  1,  1, /*  10: (134)    p19 */
512     3,  1,  1,  2,  0,  1,  2,  1  /*  11: (124)    p20 */};
513   /* Each rotation determines a permutation of the four diagonals, and this defines the isomorphism with S_4 */
514   static const PetscInt hexArr[48*6*2] = {
515     2, -3,  3, -2,  4, -2,  5, -3,  1, -3,  0, -1, /* -24: reflect bottom and use -3 on top */
516     4, -2,  5, -2,  0, -1,  1, -4,  3, -2,  2, -3, /* -23: reflect bottom and use -3 on top */
517     5, -3,  4, -1,  1, -2,  0, -3,  3, -4,  2, -1, /* -22: reflect bottom and use -3 on top */
518     3, -1,  2, -4,  4, -4,  5, -1,  0, -4,  1, -4, /* -21: reflect bottom and use -3 on top */
519     3, -3,  2, -2,  5, -1,  4, -4,  1, -1,  0, -3, /* -20: reflect bottom and use -3 on top */
520     4, -4,  5, -4,  1, -4,  0, -1,  2, -4,  3, -1, /* -19: reflect bottom and use -3 on top */
521     2, -1,  3, -4,  5, -3,  4, -2,  0, -2,  1, -2, /* -18: reflect bottom and use -3 on top */
522     5, -1,  4, -3,  0, -3,  1, -2,  2, -2,  3, -3, /* -17: reflect bottom and use -3 on top */
523     4, -3,  5, -1,  3, -2,  2, -4,  1, -4,  0, -4, /* -16: reflect bottom and use -3 on top */
524     5, -4,  4, -4,  3, -4,  2, -2,  0, -3,  1, -1, /* -15: reflect bottom and use -3 on top */
525     3, -4,  2, -1,  1, -1,  0, -4,  4, -4,  5, -4, /* -14: reflect bottom and use -3 on top */
526     2, -2,  3, -3,  0, -2,  1, -3,  4, -2,  5, -2, /* -13: reflect bottom and use -3 on top */
527     1, -3,  0, -1,  4, -1,  5, -4,  3, -1,  2, -4, /* -12: reflect bottom and use -3 on top */
528     1, -1,  0, -3,  5, -4,  4, -1,  2, -1,  3, -4, /* -11: reflect bottom and use -3 on top */
529     5, -2,  4, -2,  2, -2,  3, -4,  1, -2,  0, -2, /* -10: reflect bottom and use -3 on top */
530     1, -2,  0, -2,  2, -1,  3, -1,  4, -1,  5, -3, /*  -9: reflect bottom and use -3 on top */
531     4, -1,  5, -3,  2, -4,  3, -2,  0, -1,  1, -3, /*  -8: reflect bottom and use -3 on top */
532     3, -2,  2, -3,  0, -4,  1, -1,  5, -1,  4, -3, /*  -7: reflect bottom and use -3 on top */
533     1, -4,  0, -4,  3, -1,  2, -1,  5, -4,  4, -4, /*  -6: reflect bottom and use -3 on top */
534     2, -4,  3, -1,  1, -3,  0, -2,  5, -3,  4, -1, /*  -5: reflect bottom and use -3 on top */
535     0, -4,  1, -4,  4, -3,  5, -2,  2, -3,  3, -2, /*  -4: reflect bottom and use -3 on top */
536     0, -3,  1, -1,  3, -3,  2, -3,  4, -3,  5, -1, /*  -3: reflect bottom and use -3 on top */
537     0, -2,  1, -2,  5, -2,  4, -3,  3, -3,  2, -2, /*  -2: reflect bottom and use -3 on top */
538     0, -1,  1, -3,  2, -3,  3, -3,  5, -2,  4, -2, /*  -1: reflect bottom and use -3 on top */
539     0,  0,  1,  0,  2,  0,  3,  0,  4,  0,  5,  0, /*   0: identity */
540     0,  1,  1,  3,  5,  3,  4,  0,  2,  0,  3,  1, /*   1: 90  rotation about z */
541     0,  2,  1,  2,  3,  0,  2,  0,  5,  3,  4,  1, /*   2: 180 rotation about z */
542     0,  3,  1,  1,  4,  0,  5,  3,  3,  0,  2,  1, /*   3: 270 rotation about z */
543     2,  3,  3,  2,  1,  0,  0,  3,  4,  3,  5,  1, /*   4: 90  rotation about x */
544     1,  3,  0,  1,  3,  2,  2,  2,  4,  2,  5,  2, /*   5: 180 rotation about x */
545     3,  1,  2,  0,  0,  1,  1,  2,  4,  1,  5,  3, /*   6: 270 rotation about x */
546     4,  0,  5,  0,  2,  1,  3,  3,  1,  1,  0,  3, /*   7: 90  rotation about y */
547     1,  1,  0,  3,  2,  2,  3,  2,  5,  1,  4,  3, /*   8: 180 rotation about y */
548     5,  1,  4,  3,  2,  3,  3,  1,  0,  0,  1,  0, /*   9: 270 rotation about y */
549     1,  0,  0,  0,  5,  1,  4,  2,  3,  2,  2,  3, /*  10: 180 rotation about x+y */
550     1,  2,  0,  2,  4,  2,  5,  1,  2,  2,  3,  3, /*  11: 180 rotation about x-y */
551     2,  1,  3,  0,  0,  3,  1,  0,  5,  0,  4,  0, /*  12: 180 rotation about y+z */
552     3,  3,  2,  2,  1,  2,  0,  1,  5,  2,  4,  2, /*  13: 180 rotation about y-z */
553     5,  3,  4,  1,  3,  1,  2,  3,  1,  3,  0,  1, /*  14: 180 rotation about z+x */
554     4,  2,  5,  2,  3,  3,  2,  1,  0,  2,  1,  2, /*  15: 180 rotation about z-x */
555     5,  0,  4,  0,  0,  0,  1,  3,  3,  1,  2,  0, /*  16: 120 rotation about x+y+z (v0v6) */
556     2,  0,  3,  1,  5,  0,  4,  3,  1,  0,  0,  0, /*  17: 240 rotation about x+y+z (v0v6) */
557     4,  3,  5,  1,  1,  1,  0,  2,  3,  3,  2,  2, /*  18: 120 rotation about x+y-z (v4v2) */
558     3,  2,  2,  3,  5,  2,  4,  1,  0,  1,  1,  3, /*  19: 240 rotation about x+y-z (v4v2) */
559     3,  0,  2,  1,  4,  1,  5,  2,  1,  2,  0,  2, /*  20: 120 rotation about x-y+z (v1v5) */
560     5,  2,  4,  2,  1,  3,  0,  0,  2,  3,  3,  2, /*  21: 240 rotation about x-y+z (v1v5) */
561     4,  1,  5,  3,  0,  2,  1,  1,  2,  1,  3,  0, /*  22: 120 rotation about x-y-z (v7v3) */
562     2,  2,  3,  3,  4,  3,  5,  0,  0,  3,  1,  1, /*  23: 240 rotation about x-y-z (v7v3) */
563   };
564   static const PetscInt tripArr[12*5*2] = {
565     1, -3,  0, -1,  3, -1,  4, -1,  2, -1, /* -6: reflect bottom and top */
566     1, -1,  0, -3,  4, -1,  2, -1,  3, -1, /* -5: reflect bottom and top */
567     1, -2,  0, -2,  2, -1,  3, -1,  4, -1, /* -4: reflect bottom and top */
568     0, -3,  1, -1,  3, -3,  2, -3,  4, -3, /* -3: reflect bottom and top */
569     0, -2,  1, -2,  4, -3,  3, -3,  2, -3, /* -2: reflect bottom and top */
570     0, -1,  1, -3,  2, -3,  4, -3,  3, -3, /* -1: reflect bottom and top */
571     0,  0,  1,  0,  2,  0,  3,  0,  4,  0, /*  0: identity */
572     0,  1,  1,  2,  4,  0,  2,  0,  3,  0, /*  1: 120 rotation about z */
573     0,  2,  1,  1,  3,  0,  4,  0,  2,  0, /*  2: 240 rotation about z */
574     1,  1,  0,  2,  2,  2,  4,  2,  3,  2, /*  3: 180 rotation about y of 0 */
575     1,  0,  0,  0,  4,  2,  3,  2,  2,  2, /*  4: 180 rotation about y of 1 */
576     1,  2,  0,  1,  3,  2,  2,  2,  4,  2, /*  5: 180 rotation about y of 2 */
577   };
578   /* a: rotate 120 about z
579      b: swap top and bottom segments
580      r: reflect */
581   static const PetscInt ttriArr[12*5*2] = {
582     1, -3,  0, -3,  2, -2,  4, -2,  3, -2, /* -6: r b a^2 */
583     1, -2,  0, -2,  4, -2,  3, -2,  2, -2, /* -5: r b a */
584     1, -1,  0, -1,  3, -2,  2, -2,  4, -2, /* -4: r b */
585     0, -3,  1, -3,  2, -1,  4, -1,  3, -1, /* -3: r a^2 */
586     0, -2,  1, -2,  4, -1,  3, -1,  2, -1, /* -2: r a */
587     0, -1,  1, -1,  3, -1,  2, -1,  4, -1, /* -1: r */
588     0,  0,  1,  0,  2,  0,  3,  0,  4,  0, /*  0: identity */
589     0,  1,  1,  1,  3,  0,  4,  0,  2,  0, /*  1: a */
590     0,  2,  1,  2,  4,  0,  2,  0,  3,  0, /*  2: a^2 */
591     1,  0,  0,  0,  2,  1,  3,  1,  4,  1, /*  3: b */
592     1,  1,  0,  1,  3,  1,  4,  1,  2,  1, /*  4: b a */
593     1,  2,  0,  2,  4,  1,  2,  1,  3,  1, /*  5: b a^2 */
594   };
595   /* a: rotate 90 about z
596      b: swap top and bottom segments
597      r: reflect */
598   static const PetscInt tquadArr[16*6*2] = {
599     1, -4,  0, -4,  3, -2,  2, -2,  5, -2,  4, -2, /* -8: r b a^3 */
600     1, -3,  0, -3,  2, -2,  5, -2,  4, -2,  3, -2, /* -7: r b a^2 */
601     1, -2,  0, -2,  5, -2,  4, -2,  3, -2,  2, -2, /* -6: r b a */
602     1, -1,  0, -1,  4, -2,  3, -2,  2, -2,  5, -2, /* -5: r b */
603     0, -4,  1, -4,  3, -1,  2, -1,  5, -1,  4, -1, /* -4: r a^3 */
604     0, -3,  1, -3,  2, -1,  5, -1,  4, -1,  3, -1, /* -3: r a^2 */
605     0, -2,  1, -2,  5, -1,  4, -1,  3, -1,  2, -1, /* -2: r a */
606     0, -1,  1, -1,  4, -1,  3, -1,  2, -1,  5, -1, /* -1: r */
607     0,  0,  1,  0,  2,  0,  3,  0,  4,  0,  5,  0, /*  0: identity */
608     0,  1,  1,  1,  3,  0,  4,  0,  5,  0,  2,  0, /*  1: a */
609     0,  2,  1,  2,  4,  0,  5,  0,  2,  0,  3,  0, /*  2: a^2 */
610     0,  3,  1,  3,  5,  0,  2,  0,  3,  0,  4,  0, /*  3: a^3 */
611     1,  0,  0,  0,  2,  1,  3,  1,  4,  1,  5,  1, /*  4: b */
612     1,  1,  0,  1,  3,  1,  4,  1,  5,  1,  2,  1, /*  5: b a */
613     1,  2,  0,  2,  4,  1,  5,  1,  2,  1,  3,  1, /*  6: b a^2 */
614     1,  3,  0,  3,  5,  1,  2,  1,  3,  1,  4,  1, /*  7: b a^3 */
615   };
616   static const PetscInt pyrArr[8*5*2] = {
617     0, -4,  2, -3,  1, -3,  4, -3,  3, -3, /* -4: Reflect bottom face */
618     0, -3,  3, -3,  2, -3,  1, -3,  4, -3, /* -3: Reflect bottom face */
619     0, -2,  4, -3,  3, -3,  2, -3,  1, -3, /* -2: Reflect bottom face */
620     0, -1,  1, -3,  4, -3,  3, -3,  2, -3, /* -1: Reflect bottom face */
621     0,  0,  1,  0,  2,  0,  3,  0,  4,  0, /*  0: identity */
622     0,  1,  4,  0,  1,  0,  2,  0,  3,  0, /*  1:  90 rotation about z */
623     0,  2,  3,  0,  4,  0,  1,  0,  2,  0, /*  2: 180 rotation about z */
624     0,  3,  2,  0,  3,  0,  4,  0,  1,  0, /*  3: 270 rotation about z */
625   };
626   switch (ct) {
627     case DM_POLYTOPE_POINT:              return pntArr;
628     case DM_POLYTOPE_SEGMENT:            return &segArr[(o+1)*2*2];
629     case DM_POLYTOPE_POINT_PRISM_TENSOR: return &segArr[(o+1)*2*2];
630     case DM_POLYTOPE_TRIANGLE:           return &triArr[(o+3)*3*2];
631     case DM_POLYTOPE_QUADRILATERAL:      return &quadArr[(o+4)*4*2];
632     case DM_POLYTOPE_SEG_PRISM_TENSOR:   return &tsegArr[(o+2)*4*2];
633     case DM_POLYTOPE_TETRAHEDRON:        return &tetArr[(o+12)*4*2];
634     case DM_POLYTOPE_HEXAHEDRON:         return &hexArr[(o+24)*6*2];
635     case DM_POLYTOPE_TRI_PRISM:          return &tripArr[(o+6)*5*2];
636     case DM_POLYTOPE_TRI_PRISM_TENSOR:   return &ttriArr[(o+6)*5*2];
637     case DM_POLYTOPE_QUAD_PRISM_TENSOR:  return &tquadArr[(o+8)*6*2];
638     case DM_POLYTOPE_PYRAMID:            return &pyrArr[(o+4)*5*2];
639     default: return NULL;
640   }
641 }
642 
643 /* A vertex arrangment is a vertex order */
644 PETSC_STATIC_INLINE const PetscInt *DMPolytopeTypeGetVertexArrangment(DMPolytopeType ct, PetscInt o)
645 {
646   static const PetscInt pntVerts[1]    = {0};
647   static const PetscInt segVerts[2*2]  = {
648     1, 0,
649     0, 1};
650   static const PetscInt triVerts[6*3]  = {
651     1, 0, 2,
652     0, 2, 1,
653     2, 1, 0,
654     0, 1, 2,
655     1, 2, 0,
656     2, 0, 1};
657   static const PetscInt quadVerts[8*4]  = {
658     2, 1, 0, 3,
659     1, 0, 3, 2,
660     0, 3, 2, 1,
661     3, 2, 1, 0,
662     0, 1, 2, 3,
663     1, 2, 3, 0,
664     2, 3, 0, 1,
665     3, 0, 1, 2};
666   static const PetscInt tsegVerts[4*4]  = {
667     3, 2, 1, 0,
668     1, 0, 3, 2,
669     0, 1, 2, 3,
670     2, 3, 0, 1};
671   static const PetscInt tetVerts[24*4] = {
672     2, 3, 1, 0, /* -12: (1324)   p22 */
673     3, 1, 2, 0, /* -11: (14)     p21 */
674     1, 2, 3, 0, /* -10: (1234)   p18 */
675     3, 2, 0, 1, /*  -9: (1423)   p17 */
676     2, 0, 3, 1, /*  -8: (1342)   p13 */
677     0, 3, 2, 1, /*  -7: (24)     p14 */
678     0, 1, 3, 2, /*  -6: (34)     p6  */
679     1, 3, 0, 2, /*  -5: (1243)   p10 */
680     3, 0, 1, 2, /*  -4: (1432    p9  */
681     1, 0, 2, 3, /*  -3: (12)     p1  */
682     0, 2, 1, 3, /*  -2: (23)     p2  */
683     2, 1, 0, 3, /*  -1: (13)     p5  */
684     0, 1, 2, 3, /*   0: ()       p0  */
685     1, 2, 0, 3, /*   1: (123)    p4  */
686     2, 0, 1, 3, /*   2: (132)    p3  */
687     1, 0, 3, 2, /*   3: (12)(34) p7  */
688     0, 3, 1, 2, /*   4: (243)    p8  */
689     3, 1, 0, 2, /*   5: (143)    p11 */
690     2, 3, 0, 1, /*   6: (13)(24) p16 */
691     3, 0, 2, 1, /*   7: (142)    p15 */
692     0, 2, 3, 1, /*   8: (234)    p12 */
693     3, 2, 1, 0, /*   9: (14)(23) p23 */
694     2, 1, 3, 0, /*  10: (134)    p19 */
695     1, 3, 2, 0  /*  11: (124)    p20 */};
696   static const PetscInt hexVerts[48*8] = {
697     3,  0,  4,  5,  2,  6,  7,  1, /* -24: reflected 23 */
698     3,  5,  6,  2,  0,  1,  7,  4, /* -23: reflected 22 */
699     4,  0,  1,  7,  5,  6,  2,  3, /* -22: reflected 21 */
700     6,  7,  1,  2,  5,  3,  0,  4, /* -21: reflected 20 */
701     1,  2,  6,  7,  0,  4,  5,  3, /* -20: reflected 19 */
702     6,  2,  3,  5,  7,  4,  0,  1, /* -19: reflected 18 */
703     4,  5,  3,  0,  7,  1,  2,  6, /* -18: reflected 17 */
704     1,  7,  4,  0,  2,  3,  5,  6, /* -17: reflected 16 */
705     2,  3,  5,  6,  1,  7,  4,  0, /* -16: reflected 15 */
706     7,  4,  0,  1,  6,  2,  3,  5, /* -15: reflected 14 */
707     7,  1,  2,  6,  4,  5,  3,  0, /* -14: reflected 13 */
708     0,  4,  5,  3,  1,  2,  6,  7, /* -13: reflected 12 */
709     5,  4,  7,  6,  3,  2,  1,  0, /* -12: reflected 11 */
710     7,  6,  5,  4,  1,  0,  3,  2, /* -11: reflected 10 */
711     0,  1,  7,  4,  3,  5,  6,  2, /* -10: reflected  9 */
712     4,  7,  6,  5,  0,  3,  2,  1, /*  -9: reflected  8 */
713     5,  6,  2,  3,  4,  0,  1,  7, /*  -8: reflected  7 */
714     2,  6,  7,  1,  3,  0,  4,  5, /*  -7: reflected  6 */
715     6,  5,  4,  7,  2,  1,  0,  3, /*  -6: reflected  5 */
716     5,  3,  0,  4,  6,  7,  1,  2, /*  -5: reflected  4 */
717     2,  1,  0,  3,  6,  5,  4,  7, /*  -4: reflected  3 */
718     1,  0,  3,  2,  7,  6,  5,  4, /*  -3: reflected  2 */
719     0,  3,  2,  1,  4,  7,  6,  5, /*  -2: reflected  1 */
720     3,  2,  1,  0,  5,  4,  7,  6, /*  -1: reflected  0 */
721     0,  1,  2,  3,  4,  5,  6,  7, /*   0: identity */
722     1,  2,  3,  0,  7,  4,  5,  6, /*   1: 90  rotation about z */
723     2,  3,  0,  1,  6,  7,  4,  5, /*   2: 180 rotation about z */
724     3,  0,  1,  2,  5,  6,  7,  4, /*   3: 270 rotation about z */
725     4,  0,  3,  5,  7,  6,  2,  1, /*   4: 90  rotation about x */
726     7,  4,  5,  6,  1,  2,  3,  0, /*   5: 180 rotation about x */
727     1,  7,  6,  2,  0,  3,  5,  4, /*   6: 270 rotation about x */
728     3,  2,  6,  5,  0,  4,  7,  1, /*   7: 90  rotation about y */
729     5,  6,  7,  4,  3,  0,  1,  2, /*   8: 180 rotation about y */
730     4,  7,  1,  0,  5,  3,  2,  6, /*   9: 270 rotation about y */
731     4,  5,  6,  7,  0,  1,  2,  3, /*  10: 180 rotation about x+y */
732     6,  7,  4,  5,  2,  3,  0,  1, /*  11: 180 rotation about x-y */
733     3,  5,  4,  0,  2,  1,  7,  6, /*  12: 180 rotation about y+z */
734     6,  2,  1,  7,  5,  4,  0,  3, /*  13: 180 rotation about y-z */
735     1,  0,  4,  7,  2,  6,  5,  3, /*  14: 180 rotation about z+x */
736     6,  5,  3,  2,  7,  1,  0,  4, /*  15: 180 rotation about z-x */
737     0,  4,  7,  1,  3,  2,  6,  5, /*  16: 120 rotation about x+y+z (v0v6) */
738     0,  3,  5,  4,  1,  7,  6,  2, /*  17: 240 rotation about x+y+z (v0v6) */
739     5,  3,  2,  6,  4,  7,  1,  0, /*  18: 120 rotation about x+y-z (v4v2) */
740     7,  6,  2,  1,  4,  0,  3,  5, /*  19: 240 rotation about x+y-z (v4v2) */
741     2,  1,  7,  6,  3,  5,  4,  0, /*  20: 120 rotation about x-y+z (v1v5) */
742     7,  1,  0,  4,  6,  5,  3,  2, /*  21: 240 rotation about x-y+z (v1v5) */
743     2,  6,  5,  3,  1,  0,  4,  7, /*  22: 120 rotation about x-y-z (v7v3) */
744     5,  4,  0,  3,  6,  2,  1,  7, /*  23: 240 rotation about x-y-z (v7v3) */
745   };
746   static const PetscInt tripVerts[12*6] = {
747     4,  3,  5,  2,  1,  0, /* -6: reflect bottom and top */
748     5,  4,  3,  1,  0,  2, /* -5: reflect bottom and top */
749     3,  5,  4,  0,  2,  1, /* -4: reflect bottom and top */
750     1,  0,  2,  5,  4,  3, /* -3: reflect bottom and top */
751     0,  2,  1,  3,  5,  4, /* -2: reflect bottom and top */
752     2,  1,  0,  4,  3,  5, /* -1: reflect bottom and top */
753     0,  1,  2,  3,  4,  5, /*  0: identity */
754     1,  2,  0,  5,  3,  4, /*  1: 120 rotation about z */
755     2,  0,  1,  4,  5,  3, /*  2: 240 rotation about z */
756     4,  5,  3,  2,  0,  1, /*  3: 180 rotation about y of 0 */
757     3,  4,  5,  0,  1,  2, /*  4: 180 rotation about y of 1 */
758     5,  3,  4,  1,  2,  0, /*  5: 180 rotation about y of 2 */
759   };
760   static const PetscInt ttriVerts[12*6] = {
761     4,  3,  5,  1,  0,  2, /* -6: r b a^2 */
762     3,  5,  4,  0,  2,  1, /* -5: r b a */
763     5,  4,  3,  2,  1,  0, /* -4: r b */
764     1,  0,  2,  4,  3,  5, /* -3: r a^2 */
765     0,  2,  1,  3,  5,  4, /* -2: r a */
766     2,  1,  0,  5,  4,  3, /* -1: r */
767     0,  1,  2,  3,  4,  5, /*  0: identity */
768     1,  2,  0,  4,  5,  3, /*  1: a */
769     2,  0,  1,  5,  3,  4, /*  2: a^2 */
770     3,  4,  5,  0,  1,  2, /*  3: b */
771     4,  5,  3,  1,  2,  0, /*  4: b a */
772     5,  3,  4,  2,  0,  1, /*  5: b a^2 */
773   };
774   /* a: rotate 90 about z
775      b: swap top and bottom segments
776      r: reflect */
777   static const PetscInt tquadVerts[16*8] = {
778     6,  5,  4,  7,  2,  1,  0,  3, /* -8: r b a^3 */
779     5,  4,  7,  6,  1,  0,  3,  2, /* -7: r b a^2 */
780     4,  7,  6,  5,  0,  3,  2,  1, /* -6: r b a */
781     7,  6,  5,  4,  3,  2,  1,  0, /* -5: r b */
782     2,  1,  0,  3,  6,  5,  4,  7, /* -4: r a^3 */
783     1,  0,  3,  2,  5,  4,  7,  6, /* -3: r a^2 */
784     0,  3,  2,  1,  4,  7,  6,  5, /* -2: r a */
785     3,  2,  1,  0,  7,  6,  5,  4, /* -1: r */
786     0,  1,  2,  3,  4,  5,  6,  7, /*  0: identity */
787     1,  2,  3,  0,  5,  6,  7,  4, /*  1: a */
788     2,  3,  0,  1,  6,  7,  4,  5, /*  2: a^2 */
789     3,  0,  1,  2,  7,  4,  5,  6, /*  3: a^3 */
790     4,  5,  6,  7,  0,  1,  2,  3, /*  4: b */
791     5,  6,  7,  4,  1,  2,  3,  0, /*  5: b a */
792     6,  7,  4,  5,  2,  3,  0,  1, /*  6: b a^2 */
793     7,  4,  5,  6,  3,  0,  1,  2, /*  7: b a^3 */
794   };
795   static const PetscInt pyrVerts[8*5] = {
796     2,  1,  0,  3,  4, /* -4: Reflect bottom face */
797     1,  0,  3,  2,  4, /* -3: Reflect bottom face */
798     0,  3,  2,  1,  4, /* -2: Reflect bottom face */
799     3,  2,  1,  0,  4, /* -1: Reflect bottom face */
800     0,  1,  2,  3,  4, /*  0: identity */
801     1,  2,  3,  0,  4, /*  1:  90 rotation about z */
802     2,  3,  0,  1,  4, /*  2: 180 rotation about z */
803     3,  0,  1,  2,  4, /*  3: 270 rotation about z */
804   };
805   switch (ct) {
806     case DM_POLYTOPE_POINT:              return pntVerts;
807     case DM_POLYTOPE_SEGMENT:            return &segVerts[(o+1)*2];
808     case DM_POLYTOPE_POINT_PRISM_TENSOR: return &segVerts[(o+1)*2];
809     case DM_POLYTOPE_TRIANGLE:           return &triVerts[(o+3)*3];
810     case DM_POLYTOPE_QUADRILATERAL:      return &quadVerts[(o+4)*4];
811     case DM_POLYTOPE_SEG_PRISM_TENSOR:   return &tsegVerts[(o+2)*4];
812     case DM_POLYTOPE_TETRAHEDRON:        return &tetVerts[(o+12)*4];
813     case DM_POLYTOPE_HEXAHEDRON:         return &hexVerts[(o+24)*8];
814     case DM_POLYTOPE_TRI_PRISM:          return &tripVerts[(o+6)*6];
815     case DM_POLYTOPE_TRI_PRISM_TENSOR:   return &ttriVerts[(o+6)*6];
816     case DM_POLYTOPE_QUAD_PRISM_TENSOR:  return &tquadVerts[(o+8)*8];
817     case DM_POLYTOPE_PYRAMID:            return &pyrVerts[(o+4)*5];
818     default: return NULL;
819   }
820 }
821 
822 /* This is orientation o1 acting on orientation o2 */
823 PETSC_STATIC_INLINE PetscInt DMPolytopeTypeComposeOrientation(DMPolytopeType ct, PetscInt o1, PetscInt o2)
824 {
825   static const PetscInt segMult[2*2] = {
826      0, -1,
827     -1,  0};
828   static const PetscInt triMult[6*6] = {
829      0,  2,  1, -3, -1, -2,
830      1,  0,  2, -2, -3, -1,
831      2,  1,  0, -1, -2, -3,
832     -3, -2, -1,  0,  1,  2,
833     -2, -1, -3,  1,  2,  0,
834     -1, -3, -2,  2,  0,  1};
835   static const PetscInt quadMult[8*8] = {
836      0,  3,  2,  1, -4, -1, -2, -3,
837      1,  0,  3,  2, -3, -4, -1, -2,
838      2,  1,  0,  3, -2, -3, -4, -1,
839      3,  2,  1,  0, -1, -2, -3, -4,
840     -4, -3, -2, -1,  0,  1,  2,  3,
841     -3, -2, -1, -4,  1,  2,  3,  0,
842     -2, -1, -4, -3,  2,  3,  0,  1,
843     -1, -4, -3, -2,  3,  0,  1,  2};
844   static const PetscInt tsegMult[4*4] = {
845      0,  1, -2, -1,
846      1,  0, -1, -2,
847     -2, -1,  0,  1,
848     -1, -2,  1,  0};
849   static const PetscInt tetMult[24*24] = {
850     3, 2, 7, 0, 5, 10, 9, 8, 1, 6, 11, 4, -12, -7, -5, -9, -10, -2, -6, -1, -11, -3, -4, -8,
851     4, 0, 8, 1, 3, 11, 10, 6, 2, 7, 9, 5, -11, -9, -4, -8, -12, -1, -5, -3, -10, -2, -6, -7,
852     5, 1, 6, 2, 4, 9, 11, 7, 0, 8, 10, 3, -10, -8, -6, -7, -11, -3, -4, -2, -12, -1, -5, -9,
853     0, 8, 4, 3, 11, 1, 6, 2, 10, 9, 5, 7, -9, -4, -11, -12, -1, -8, -3, -10, -5, -6, -7, -2,
854     1, 6, 5, 4, 9, 2, 7, 0, 11, 10, 3, 8, -8, -6, -10, -11, -3, -7, -2, -12, -4, -5, -9, -1,
855     2, 7, 3, 5, 10, 0, 8, 1, 9, 11, 4, 6, -7, -5, -12, -10, -2, -9, -1, -11, -6, -4, -8, -3,
856     6, 5, 1, 9, 2, 4, 0, 11, 7, 3, 8, 10, -6, -10, -8, -3, -7, -11, -12, -4, -2, -9, -1, -5,
857     7, 3, 2, 10, 0, 5, 1, 9, 8, 4, 6, 11, -5, -12, -7, -2, -9, -10, -11, -6, -1, -8, -3, -4,
858     8, 4, 0, 11, 1, 3, 2, 10, 6, 5, 7, 9, -4, -11, -9, -1, -8, -12, -10, -5, -3, -7, -2, -6,
859     9, 11, 10, 6, 8, 7, 3, 5, 4, 0, 2, 1, -3, -1, -2, -6, -4, -5, -9, -7, -8, -12, -10, -11,
860     10, 9, 11, 7, 6, 8, 4, 3, 5, 1, 0, 2, -2, -3, -1, -5, -6, -4, -8, -9, -7, -11, -12, -10,
861     11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, -1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12,
862     -12, -11, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
863     -11, -10, -12, -8, -7, -9, -5, -4, -6, -2, -1, -3, 1, 2, 0, 4, 5, 3, 7, 8, 6, 10, 11, 9,
864     -10, -12, -11, -7, -9, -8, -4, -6, -5, -1, -3, -2, 2, 0, 1, 5, 3, 4, 8, 6, 7, 11, 9, 10,
865     -9, -5, -1, -12, -2, -4, -3, -11, -7, -6, -8, -10, 3, 10, 8, 0, 7, 11, 9, 4, 2, 6, 1, 5,
866     -8, -4, -3, -11, -1, -6, -2, -10, -9, -5, -7, -12, 4, 11, 6, 1, 8, 9, 10, 5, 0, 7, 2, 3,
867     -7, -6, -2, -10, -3, -5, -1, -12, -8, -4, -9, -11, 5, 9, 7, 2, 6, 10, 11, 3, 1, 8, 0, 4,
868     -3, -8, -4, -6, -11, -1, -9, -2, -10, -12, -5, -7, 6, 4, 11, 9, 1, 8, 0, 10, 5, 3, 7, 2,
869     -2, -7, -6, -5, -10, -3, -8, -1, -12, -11, -4, -9, 7, 5, 9, 10, 2, 6, 1, 11, 3, 4, 8, 0,
870     -1, -9, -5, -4, -12, -2, -7, -3, -11, -10, -6, -8, 8, 3, 10, 11, 0, 7, 2, 9, 4, 5, 6, 1,
871     -6, -2, -7, -3, -5, -10, -12, -8, -1, -9, -11, -4, 9, 7, 5, 6, 10, 2, 3, 1, 11, 0, 4, 8,
872     -5, -1, -9, -2, -4, -12, -11, -7, -3, -8, -10, -6, 10, 8, 3, 7, 11, 0, 4, 2, 9, 1, 5, 6,
873     -4, -3, -8, -1, -6, -11, -10, -9, -2, -7, -12, -5, 11, 6, 4, 8, 9, 1, 5, 0, 10, 2, 3, 7,
874     };
875   static const PetscInt hexMult[48*48] = {
876     18, 2, 5, 22, 21, 8, 16, 0, 13, 6, 11, 3, 15, 9, 4, 23, 12, 1, 19, 10, 7, 20, 14, 17, -24, -10, -20, -16, -12, -21, -4, -5, -18, -13, -15, -8, -2, -11, -14, -7, -3, -22, -6, -17, -19, -9, -1, -23,
877     8, 20, 19, 2, 5, 23, 0, 17, 11, 1, 15, 7, 13, 4, 10, 18, 3, 14, 21, 9, 12, 22, 6, 16, -23, -13, -17, -7, -8, -19, -16, -12, -22, -2, -14, -5, -10, -15, -11, -4, -20, -9, -21, -3, -6, -18, -24, -1,
878     2, 17, 23, 8, 0, 19, 5, 20, 1, 11, 9, 14, 12, 6, 3, 16, 10, 7, 22, 15, 13, 21, 4, 18, -22, -14, -19, -5, -15, -17, -10, -2, -23, -12, -13, -7, -16, -8, -4, -11, -24, -3, -18, -9, -1, -21, -20, -6,
879     21, 5, 2, 16, 18, 0, 22, 8, 4, 12, 3, 11, 14, 7, 13, 20, 6, 10, 17, 1, 9, 23, 15, 19, -21, -8, -18, -15, -4, -24, -12, -14, -20, -7, -16, -10, -11, -2, -5, -13, -6, -19, -3, -23, -22, -1, -9, -17,
880     16, 8, 0, 21, 22, 2, 18, 5, 12, 4, 1, 10, 9, 15, 6, 19, 13, 11, 23, 3, 14, 17, 7, 20, -20, -16, -24, -10, -2, -18, -11, -7, -21, -14, -8, -15, -12, -4, -13, -5, -9, -23, -1, -19, -17, -3, -6, -22,
881     5, 19, 20, 0, 8, 17, 2, 23, 10, 3, 7, 15, 6, 12, 11, 22, 1, 9, 16, 14, 4, 18, 13, 21, -19, -5, -22, -14, -16, -23, -8, -11, -17, -4, -7, -13, -15, -10, -12, -2, -21, -6, -20, -1, -9, -24, -18, -3,
882     22, 0, 8, 18, 16, 5, 21, 2, 6, 13, 10, 1, 7, 14, 12, 17, 4, 3, 20, 11, 15, 19, 9, 23, -18, -15, -21, -8, -11, -20, -2, -13, -24, -5, -10, -16, -4, -12, -7, -14, -1, -17, -9, -22, -23, -6, -3, -19,
883     0, 23, 17, 5, 2, 20, 8, 19, 3, 10, 14, 9, 4, 13, 1, 21, 11, 15, 18, 7, 6, 16, 12, 22, -17, -7, -23, -13, -10, -22, -15, -4, -19, -11, -5, -14, -8, -16, -2, -12, -18, -1, -24, -6, -3, -20, -21, -9,
884     10, 13, 6, 1, 11, 12, 3, 4, 8, 0, 22, 18, 19, 23, 5, 15, 2, 21, 9, 16, 17, 7, 20, 14, -16, -24, -10, -20, -23, -8, -19, -6, -15, -3, -21, -18, -22, -17, -9, -1, -14, -12, -7, -4, -11, -13, -5, -2,
885     1, 4, 12, 10, 3, 6, 11, 13, 0, 8, 16, 21, 17, 20, 2, 14, 5, 18, 7, 22, 19, 9, 23, 15, -15, -21, -8, -18, -17, -10, -22, -3, -16, -6, -24, -20, -19, -23, -1, -9, -5, -4, -13, -12, -2, -7, -14, -11,
886     14, 10, 3, 9, 7, 1, 15, 11, 17, 23, 0, 5, 16, 22, 20, 6, 19, 8, 12, 2, 21, 4, 18, 13, -14, -19, -5, -22, -3, -13, -9, -20, -7, -21, -23, -17, -6, -1, -24, -18, -12, -16, -2, -8, -10, -4, -11, -15,
887     7, 3, 10, 15, 14, 11, 9, 1, 20, 19, 5, 0, 18, 21, 17, 4, 23, 2, 13, 8, 22, 6, 16, 12, -13, -17, -7, -23, -9, -14, -3, -24, -5, -18, -22, -19, -1, -6, -20, -21, -2, -10, -12, -15, -16, -11, -4, -8,
888     13, 14, 15, 12, 4, 9, 6, 7, 21, 22, 23, 20, 2, 0, 18, 3, 16, 17, 1, 19, 8, 11, 5, 10, -12, -9, -11, -6, -21, -4, -24, -22, -2, -23, -3, -1, -20, -18, -19, -17, -16, -14, -15, -13, -5, -8, -10, -7,
889     6, 9, 7, 4, 12, 14, 13, 15, 16, 18, 17, 19, 0, 2, 22, 1, 21, 23, 3, 20, 5, 10, 8, 11, -11, -6, -12, -9, -20, -2, -18, -17, -4, -19, -1, -3, -21, -24, -23, -22, -8, -7, -10, -5, -13, -16, -15, -14,
890     3, 12, 4, 11, 1, 13, 10, 6, 2, 5, 21, 16, 23, 19, 0, 9, 8, 22, 15, 18, 20, 14, 17, 7, -10, -20, -16, -24, -22, -15, -17, -1, -8, -9, -18, -21, -23, -19, -3, -6, -13, -2, -5, -11, -4, -14, -7, -12,
891     20, 16, 18, 23, 17, 21, 19, 22, 14, 15, 4, 6, 3, 1, 7, 0, 9, 12, 2, 13, 11, 5, 10, 8, -9, -11, -6, -12, -14, -3, -13, -10, -1, -8, -2, -4, -7, -5, -16, -15, -23, -20, -22, -18, -24, -19, -17, -21,
892     11, 6, 13, 3, 10, 4, 1, 12, 5, 2, 18, 22, 20, 17, 8, 7, 0, 16, 14, 21, 23, 15, 19, 9, -8, -18, -15, -21, -19, -16, -23, -9, -10, -1, -20, -24, -17, -22, -6, -3, -7, -11, -14, -2, -12, -5, -13, -4,
893     9, 11, 1, 14, 15, 3, 7, 10, 23, 17, 2, 8, 21, 18, 19, 13, 20, 5, 4, 0, 16, 12, 22, 6, -7, -23, -13, -17, -1, -5, -6, -21, -14, -20, -19, -22, -9, -3, -18, -24, -11, -8, -4, -16, -15, -2, -12, -10,
894     19, 21, 22, 17, 23, 16, 20, 18, 9, 7, 12, 13, 1, 3, 15, 2, 14, 4, 0, 6, 10, 8, 11, 5, -6, -12, -9, -11, -7, -1, -5, -15, -3, -16, -4, -2, -14, -13, -8, -10, -19, -21, -17, -24, -18, -23, -22, -20,
895     15, 1, 11, 7, 9, 10, 14, 3, 19, 20, 8, 2, 22, 16, 23, 12, 17, 0, 6, 5, 18, 13, 21, 4, -5, -22, -14, -19, -6, -7, -1, -18, -13, -24, -17, -23, -3, -9, -21, -20, -4, -15, -11, -10, -8, -12, -2, -16,
896     4, 15, 14, 6, 13, 7, 12, 9, 18, 16, 20, 23, 5, 8, 21, 11, 22, 19, 10, 17, 0, 3, 2, 1, -4, -1, -2, -3, -24, -12, -21, -19, -11, -17, -6, -9, -18, -20, -22, -23, -15, -5, -16, -7, -14, -10, -8, -13,
897     17, 18, 16, 19, 20, 22, 23, 21, 7, 9, 6, 4, 10, 11, 14, 5, 15, 13, 8, 12, 1, 0, 3, 2, -3, -4, -1, -2, -13, -9, -14, -16, -6, -15, -12, -11, -5, -7, -10, -8, -22, -24, -23, -21, -20, -17, -19, -18,
898     12, 7, 9, 13, 6, 15, 4, 14, 22, 21, 19, 17, 8, 5, 16, 10, 18, 20, 11, 23, 2, 1, 0, 3, -2, -3, -4, -1, -18, -11, -20, -23, -12, -22, -9, -6, -24, -21, -17, -19, -10, -13, -8, -14, -7, -15, -16, -5,
899     23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, -1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12, -13, -14, -15, -16, -17, -18, -19, -20, -21, -22, -23, -24,
900     -24, -23, -22, -21, -20, -19, -18, -17, -16, -15, -14, -13, -12, -11, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
901     -13, -8, -10, -14, -7, -16, -5, -15, -23, -22, -20, -18, -9, -6, -17, -11, -19, -21, -12, -24, -3, -2, -1, -4, 1, 2, 3, 0, 17, 10, 19, 22, 11, 21, 8, 5, 23, 20, 16, 18, 9, 12, 7, 13, 6, 14, 15, 4,
902     -18, -19, -17, -20, -21, -23, -24, -22, -8, -10, -7, -5, -11, -12, -15, -6, -16, -14, -9, -13, -2, -1, -4, -3, 2, 3, 0, 1, 12, 8, 13, 15, 5, 14, 11, 10, 4, 6, 9, 7, 21, 23, 22, 20, 19, 16, 18, 17,
903     -5, -16, -15, -7, -14, -8, -13, -10, -19, -17, -21, -24, -6, -9, -22, -12, -23, -20, -11, -18, -1, -4, -3, -2, 3, 0, 1, 2, 23, 11, 20, 18, 10, 16, 5, 8, 17, 19, 21, 22, 14, 4, 15, 6, 13, 9, 7, 12,
904     -16, -2, -12, -8, -10, -11, -15, -4, -20, -21, -9, -3, -23, -17, -24, -13, -18, -1, -7, -6, -19, -14, -22, -5, 4, 21, 13, 18, 5, 6, 0, 17, 12, 23, 16, 22, 2, 8, 20, 19, 3, 14, 10, 9, 7, 11, 1, 15,
905     -20, -22, -23, -18, -24, -17, -21, -19, -10, -8, -13, -14, -2, -4, -16, -3, -15, -5, -1, -7, -11, -9, -12, -6, 5, 11, 8, 10, 6, 0, 4, 14, 2, 15, 3, 1, 13, 12, 7, 9, 18, 20, 16, 23, 17, 22, 21, 19,
906     -10, -12, -2, -15, -16, -4, -8, -11, -24, -18, -3, -9, -22, -19, -20, -14, -21, -6, -5, -1, -17, -13, -23, -7, 6, 22, 12, 16, 0, 4, 5, 20, 13, 19, 18, 21, 8, 2, 17, 23, 10, 7, 3, 15, 14, 1, 11, 9,
907     -12, -7, -14, -4, -11, -5, -2, -13, -6, -3, -19, -23, -21, -18, -9, -8, -1, -17, -15, -22, -24, -16, -20, -10, 7, 17, 14, 20, 18, 15, 22, 8, 9, 0, 19, 23, 16, 21, 5, 2, 6, 10, 13, 1, 11, 4, 12, 3,
908     -21, -17, -19, -24, -18, -22, -20, -23, -15, -16, -5, -7, -4, -2, -8, -1, -10, -13, -3, -14, -12, -6, -11, -9, 8, 10, 5, 11, 13, 2, 12, 9, 0, 7, 1, 3, 6, 4, 15, 14, 22, 19, 21, 17, 23, 18, 16, 20,
909     -4, -13, -5, -12, -2, -14, -11, -7, -3, -6, -22, -17, -24, -20, -1, -10, -9, -23, -16, -19, -21, -15, -18, -8, 9, 19, 15, 23, 21, 14, 16, 0, 7, 8, 17, 20, 22, 18, 2, 5, 12, 1, 4, 10, 3, 13, 6, 11,
910     -7, -10, -8, -5, -13, -15, -14, -16, -17, -19, -18, -20, -1, -3, -23, -2, -22, -24, -4, -21, -6, -11, -9, -12, 10, 5, 11, 8, 19, 1, 17, 16, 3, 18, 0, 2, 20, 23, 22, 21, 7, 6, 9, 4, 12, 15, 14, 13,
911     -14, -15, -16, -13, -5, -10, -7, -8, -22, -23, -24, -21, -3, -1, -19, -4, -17, -18, -2, -20, -9, -12, -6, -11, 11, 8, 10, 5, 20, 3, 23, 21, 1, 22, 2, 0, 19, 17, 18, 16, 15, 13, 14, 12, 4, 7, 9, 6,
912     -8, -4, -11, -16, -15, -12, -10, -2, -21, -20, -6, -1, -19, -22, -18, -5, -24, -3, -14, -9, -23, -7, -17, -13, 12, 16, 6, 22, 8, 13, 2, 23, 4, 17, 21, 18, 0, 5, 19, 20, 1, 9, 11, 14, 15, 10, 3, 7,
913     -15, -11, -4, -10, -8, -2, -16, -12, -18, -24, -1, -6, -17, -23, -21, -7, -20, -9, -13, -3, -22, -5, -19, -14, 13, 18, 4, 21, 2, 12, 8, 19, 6, 20, 22, 16, 5, 0, 23, 17, 11, 15, 1, 7, 9, 3, 10, 14,
914     -2, -5, -13, -11, -4, -7, -12, -14, -1, -9, -17, -22, -18, -21, -3, -15, -6, -19, -8, -23, -20, -10, -24, -16, 14, 20, 7, 17, 16, 9, 21, 2, 15, 5, 23, 19, 18, 22, 0, 8, 4, 3, 12, 11, 1, 6, 13, 10,
915     -11, -14, -7, -2, -12, -13, -4, -5, -9, -1, -23, -19, -20, -24, -6, -16, -3, -22, -10, -17, -18, -8, -21, -15, 15, 23, 9, 19, 22, 7, 18, 5, 14, 2, 20, 17, 21, 16, 8, 0, 13, 11, 6, 3, 10, 12, 4, 1,
916     -1, -24, -18, -6, -3, -21, -9, -20, -4, -11, -15, -10, -5, -14, -2, -22, -12, -16, -19, -8, -7, -17, -13, -23, 16, 6, 22, 12, 9, 21, 14, 3, 18, 10, 4, 13, 7, 15, 1, 11, 17, 0, 23, 5, 2, 19, 20, 8,
917     -23, -1, -9, -19, -17, -6, -22, -3, -7, -14, -11, -2, -8, -15, -13, -18, -5, -4, -21, -12, -16, -20, -10, -24, 17, 14, 20, 7, 10, 19, 1, 12, 23, 4, 9, 15, 3, 11, 6, 13, 0, 16, 8, 21, 22, 5, 2, 18,
918     -6, -20, -21, -1, -9, -18, -3, -24, -11, -4, -8, -16, -7, -13, -12, -23, -2, -10, -17, -15, -5, -19, -14, -22, 18, 4, 21, 13, 15, 22, 7, 10, 16, 3, 6, 12, 14, 9, 11, 1, 20, 5, 19, 0, 8, 23, 17, 2,
919     -17, -9, -1, -22, -23, -3, -19, -6, -13, -5, -2, -11, -10, -16, -7, -20, -14, -12, -24, -4, -15, -18, -8, -21, 19, 15, 23, 9, 1, 17, 10, 6, 20, 13, 7, 14, 11, 3, 12, 4, 8, 22, 0, 18, 16, 2, 5, 21,
920     -22, -6, -3, -17, -19, -1, -23, -9, -5, -13, -4, -12, -15, -8, -14, -21, -7, -11, -18, -2, -10, -24, -16, -20, 20, 7, 17, 14, 3, 23, 11, 13, 19, 6, 15, 9, 10, 1, 4, 12, 5, 18, 2, 22, 21, 0, 8, 16,
921     -3, -18, -24, -9, -1, -20, -6, -21, -2, -12, -10, -15, -13, -7, -4, -17, -11, -8, -23, -16, -14, -22, -5, -19, 21, 13, 18, 4, 14, 16, 9, 1, 22, 11, 12, 6, 15, 7, 3, 10, 23, 2, 17, 8, 0, 20, 19, 5,
922     -9, -21, -20, -3, -6, -24, -1, -18, -12, -2, -16, -8, -14, -5, -11, -19, -4, -15, -22, -10, -13, -23, -7, -17, 22, 12, 16, 6, 7, 18, 15, 11, 21, 1, 13, 4, 9, 14, 10, 3, 19, 8, 20, 2, 5, 17, 23, 0,
923     -19, -3, -6, -23, -22, -9, -17, -1, -14, -7, -12, -4, -16, -10, -5, -24, -13, -2, -20, -11, -8, -21, -15, -18, 23, 9, 19, 15, 11, 20, 3, 4, 17, 12, 14, 7, 1, 10, 13, 6, 2, 21, 5, 16, 18, 8, 0, 22,
924     };
925   static const PetscInt tripMult[12*12] = {
926     1, 0, 2, 3, 5, 4, -6, -4, -5, -2, -3, -1,
927     0, 2, 1, 4, 3, 5, -5, -6, -4, -3, -1, -2,
928     2, 1, 0, 5, 4, 3, -4, -5, -6, -1, -2, -3,
929     4, 3, 5, 0, 2, 1, -3, -1, -2, -5, -6, -4,
930     3, 5, 4, 1, 0, 2, -2, -3, -1, -6, -4, -5,
931     5, 4, 3, 2, 1, 0, -1, -2, -3, -4, -5, -6,
932     -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5,
933     -4, -6, -5, -2, -1, -3, 1, 2, 0, 5, 3, 4,
934     -5, -4, -6, -1, -3, -2, 2, 0, 1, 4, 5, 3,
935     -3, -2, -1, -6, -5, -4, 3, 4, 5, 0, 1, 2,
936     -1, -3, -2, -5, -4, -6, 4, 5, 3, 2, 0, 1,
937     -2, -1, -3, -4, -6, -5, 5, 3, 4, 1, 2, 0,
938   };
939   static const PetscInt ttriMult[12*12] = {
940     0, 2, 1, 3, 5, 4, -6, -4, -5, -3, -1, -2,
941     1, 0, 2, 4, 3, 5, -5, -6, -4, -2, -3, -1,
942     2, 1, 0, 5, 4, 3, -4, -5, -6, -1, -2, -3,
943     3, 5, 4, 0, 2, 1, -3, -1, -2, -6, -4, -5,
944     4, 3, 5, 1, 0, 2, -2, -3, -1, -5, -6, -4,
945     5, 4, 3, 2, 1, 0, -1, -2, -3, -4, -5, -6,
946     -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5,
947     -5, -4, -6, -2, -1, -3, 1, 2, 0, 4, 5, 3,
948     -4, -6, -5, -1, -3, -2, 2, 0, 1, 5, 3, 4,
949     -3, -2, -1, -6, -5, -4, 3, 4, 5, 0, 1, 2,
950     -2, -1, -3, -5, -4, -6, 4, 5, 3, 1, 2, 0,
951     -1, -3, -2, -4, -6, -5, 5, 3, 4, 2, 0, 1,
952   };
953   static const PetscInt tquadMult[16*16] = {
954     0, 3, 2, 1, 4, 7, 6, 5, -8, -5, -6, -7, -4, -1, -2, -3,
955     1, 0, 3, 2, 5, 4, 7, 6, -7, -8, -5, -6, -3, -4, -1, -2,
956     2, 1, 0, 3, 6, 5, 4, 7, -6, -7, -8, -5, -2, -3, -4, -1,
957     3, 2, 1, 0, 7, 6, 5, 4, -5, -6, -7, -8, -1, -2, -3, -4,
958     4, 7, 6, 5, 0, 3, 2, 1, -4, -1, -2, -3, -8, -5, -6, -7,
959     5, 4, 7, 6, 1, 0, 3, 2, -3, -4, -1, -2, -7, -8, -5, -6,
960     6, 5, 4, 7, 2, 1, 0, 3, -2, -3, -4, -1, -6, -7, -8, -5,
961     7, 6, 5, 4, 3, 2, 1, 0, -1, -2, -3, -4, -5, -6, -7, -8,
962     -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7,
963     -7, -6, -5, -8, -3, -2, -1, -4, 1, 2, 3, 0, 5, 6, 7, 4,
964     -6, -5, -8, -7, -2, -1, -4, -3, 2, 3, 0, 1, 6, 7, 4, 5,
965     -5, -8, -7, -6, -1, -4, -3, -2, 3, 0, 1, 2, 7, 4, 5, 6,
966     -4, -3, -2, -1, -8, -7, -6, -5, 4, 5, 6, 7, 0, 1, 2, 3,
967     -3, -2, -1, -4, -7, -6, -5, -8, 5, 6, 7, 4, 1, 2, 3, 0,
968     -2, -1, -4, -3, -6, -5, -8, -7, 6, 7, 4, 5, 2, 3, 0, 1,
969     -1, -4, -3, -2, -5, -8, -7, -6, 7, 4, 5, 6, 3, 0, 1, 2,
970   };
971   static const PetscInt pyrMult[8*8] = {
972     0, 3, 2, 1, -4, -1, -2, -3,
973     1, 0, 3, 2, -3, -4, -1, -2,
974     2, 1, 0, 3, -2, -3, -4, -1,
975     3, 2, 1, 0, -1, -2, -3, -4,
976     -4, -3, -2, -1, 0, 1, 2, 3,
977     -3, -2, -1, -4, 1, 2, 3, 0,
978     -2, -1, -4, -3, 2, 3, 0, 1,
979     -1, -4, -3, -2, 3, 0, 1, 2,
980   };
981   switch (ct) {
982     case DM_POLYTOPE_POINT:              return 0;
983     case DM_POLYTOPE_SEGMENT:
984     case DM_POLYTOPE_POINT_PRISM_TENSOR: return segMult[(o1+1)*2+o2+1];
985     case DM_POLYTOPE_TRIANGLE:           return triMult[(o1+3)*6+o2+3];
986     case DM_POLYTOPE_QUADRILATERAL:      return quadMult[(o1+4)*8+o2+4];
987     case DM_POLYTOPE_SEG_PRISM_TENSOR:   return tsegMult[(o1+2)*4+o2+2];
988     case DM_POLYTOPE_TETRAHEDRON:        return tetMult[(o1+12)*24+o2+12];
989     case DM_POLYTOPE_HEXAHEDRON:         return hexMult[(o1+24)*48+o2+24];
990     case DM_POLYTOPE_TRI_PRISM:          return tripMult[(o1+6)*12+o2+6];
991     case DM_POLYTOPE_TRI_PRISM_TENSOR:   return ttriMult[(o1+6)*12+o2+6];
992     case DM_POLYTOPE_QUAD_PRISM_TENSOR:  return tquadMult[(o1+8)*16+o2+8];
993     case DM_POLYTOPE_PYRAMID:            return pyrMult[(o1+4)*8+o2+4];
994     default: return 0;
995   }
996 }
997 
998 /* This is orientation o1 acting on orientation o2^{-1} */
999 PETSC_STATIC_INLINE PetscInt DMPolytopeTypeComposeOrientationInv(DMPolytopeType ct, PetscInt o1, PetscInt o2)
1000 {
1001   static const PetscInt triInv[6]    = {-3, -2, -1, 0, 2, 1};
1002   static const PetscInt quadInv[8]   = {-4, -3, -2, -1, 0, 3, 2, 1};
1003   static const PetscInt tetInv[24]   = {-9, -11, -4, -12, -5, -7, -6, -8, -10, -3, -2, -1, 0, 2, 1, 3, 8, 10, 6, 11, 4, 9, 5, 7};
1004   static const PetscInt hexInv[48]   = {-17, -18, -20, -19, -22, -21, -23, -24, -15, -16, -14, -13, -11, -12, -10, -9, -8, -5, -6, -7, -4, -3, -2, -1,
1005                                           0,   3,   2,   1,   6,   5,   4,   9,   8,   7,  10,  11,  12,  13,  14, 15, 17, 16, 19, 18, 21, 20, 23, 22};
1006   static const PetscInt tripInv[12]  = {-5, -6, -4, -3, -2, -1, 0, 2, 1, 3, 4, 5};
1007   static const PetscInt ttriInv[12]  = {-6, -5, -4, -3, -2, -1, 0, 2, 1, 3, 5, 4};
1008   static const PetscInt tquadInv[16] = {-8, -7, -6, -5, -4, -3, -2, -1, 0, 3, 2, 1, 4, 7, 6, 5};
1009   static const PetscInt pyrInv[8]    = {-4, -3, -2, -1, 0, 3, 2, 1};
1010   switch (ct) {
1011     case DM_POLYTOPE_POINT:              return 0;
1012     case DM_POLYTOPE_SEGMENT:
1013     case DM_POLYTOPE_POINT_PRISM_TENSOR: return DMPolytopeTypeComposeOrientation(ct, o1, o2);
1014     case DM_POLYTOPE_TRIANGLE:           return DMPolytopeTypeComposeOrientation(ct, o1, triInv[o2+3]);
1015     case DM_POLYTOPE_QUADRILATERAL:      return DMPolytopeTypeComposeOrientation(ct, o1, quadInv[o2+4]);
1016     case DM_POLYTOPE_SEG_PRISM_TENSOR:   return DMPolytopeTypeComposeOrientation(ct, o1, o2);
1017     case DM_POLYTOPE_TETRAHEDRON:        return DMPolytopeTypeComposeOrientation(ct, o1, tetInv[o2+12]);
1018     case DM_POLYTOPE_HEXAHEDRON:         return DMPolytopeTypeComposeOrientation(ct, o1, hexInv[o2+24]);
1019     case DM_POLYTOPE_TRI_PRISM:          return DMPolytopeTypeComposeOrientation(ct, o1, tripInv[o2+6]);
1020     case DM_POLYTOPE_TRI_PRISM_TENSOR:   return DMPolytopeTypeComposeOrientation(ct, o1, ttriInv[o2+6]);
1021     case DM_POLYTOPE_QUAD_PRISM_TENSOR:  return DMPolytopeTypeComposeOrientation(ct, o1, tquadInv[o2+8]);
1022     case DM_POLYTOPE_PYRAMID:            return DMPolytopeTypeComposeOrientation(ct, o1, pyrInv[o2+4]);
1023     default: return 0;
1024   }
1025 }
1026 
1027 PETSC_EXTERN PetscErrorCode DMPolytopeMatchOrientation(DMPolytopeType, const PetscInt[], const PetscInt[], PetscInt *, PetscBool *);
1028 PETSC_EXTERN PetscErrorCode DMPolytopeMatchVertexOrientation(DMPolytopeType, const PetscInt[], const PetscInt[], PetscInt *, PetscBool *);
1029 PETSC_EXTERN PetscErrorCode DMPolytopeGetOrientation(DMPolytopeType, const PetscInt[], const PetscInt[], PetscInt *);
1030 PETSC_EXTERN PetscErrorCode DMPolytopeGetVertexOrientation(DMPolytopeType, const PetscInt[], const PetscInt[], PetscInt *);
1031 PETSC_EXTERN PetscErrorCode DMPolytopeInCellTest(DMPolytopeType, const PetscReal[], PetscBool *);
1032 
1033 #endif
1034