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