Actual source code: dmpleximpl.h

      1: #pragma once

      3: #include <petscmat.h>
      4: #include <petscdmplex.h>
      5: #include <petscdmplextransform.h>
      6: #include <petscbt.h>
      7: #include <petscsf.h>
      8: #include <petsc/private/dmimpl.h>

     10: #if defined(PETSC_HAVE_EXODUSII)
     11:   #include <exodusII.h>
     12: #endif

     14: PETSC_INTERN PetscBool  Plexcite;
     15: PETSC_INTERN const char PlexCitation[];

     17: PETSC_EXTERN PetscLogEvent DMPLEX_Interpolate;
     18: PETSC_EXTERN PetscLogEvent DMPLEX_Partition;
     19: PETSC_EXTERN PetscLogEvent DMPLEX_PartSelf;
     20: PETSC_EXTERN PetscLogEvent DMPLEX_PartLabelInvert;
     21: PETSC_EXTERN PetscLogEvent DMPLEX_PartLabelCreateSF;
     22: PETSC_EXTERN PetscLogEvent DMPLEX_PartStratSF;
     23: PETSC_EXTERN PetscLogEvent DMPLEX_CreatePointSF;
     24: PETSC_EXTERN PetscLogEvent DMPLEX_Distribute;
     25: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeMultistage;
     26: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeCones;
     27: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeLabels;
     28: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeSF;
     29: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeOverlap;
     30: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeField;
     31: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeData;
     32: PETSC_EXTERN PetscLogEvent DMPLEX_Migrate;
     33: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolateSF;
     34: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalBegin;
     35: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalEnd;
     36: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalBegin;
     37: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalEnd;
     38: PETSC_EXTERN PetscLogEvent DMPLEX_Stratify;
     39: PETSC_EXTERN PetscLogEvent DMPLEX_Symmetrize;
     40: PETSC_EXTERN PetscLogEvent DMPLEX_Preallocate;
     41: PETSC_EXTERN PetscLogEvent DMPLEX_ResidualFEM;
     42: PETSC_EXTERN PetscLogEvent DMPLEX_JacobianFEM;
     43: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolatorFEM;
     44: PETSC_EXTERN PetscLogEvent DMPLEX_InjectorFEM;
     45: PETSC_EXTERN PetscLogEvent DMPLEX_IntegralFEM;
     46: PETSC_EXTERN PetscLogEvent DMPLEX_CreateGmsh;
     47: PETSC_EXTERN PetscLogEvent DMPLEX_CreateBoxSFC;
     48: PETSC_EXTERN PetscLogEvent DMPLEX_RebalanceSharedPoints;
     49: PETSC_EXTERN PetscLogEvent DMPLEX_CreateFromFile;
     50: PETSC_EXTERN PetscLogEvent DMPLEX_CreateFromOptions;
     51: PETSC_EXTERN PetscLogEvent DMPLEX_BuildFromCellList;
     52: PETSC_EXTERN PetscLogEvent DMPLEX_BuildCoordinatesFromCellList;
     53: PETSC_EXTERN PetscLogEvent DMPLEX_LocatePoints;
     54: PETSC_EXTERN PetscLogEvent DMPLEX_TopologyView;
     55: PETSC_EXTERN PetscLogEvent DMPLEX_DistributionView;
     56: PETSC_EXTERN PetscLogEvent DMPLEX_LabelsView;
     57: PETSC_EXTERN PetscLogEvent DMPLEX_CoordinatesView;
     58: PETSC_EXTERN PetscLogEvent DMPLEX_SectionView;
     59: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalVectorView;
     60: PETSC_EXTERN PetscLogEvent DMPLEX_LocalVectorView;
     61: PETSC_EXTERN PetscLogEvent DMPLEX_TopologyLoad;
     62: PETSC_EXTERN PetscLogEvent DMPLEX_DistributionLoad;
     63: PETSC_EXTERN PetscLogEvent DMPLEX_LabelsLoad;
     64: PETSC_EXTERN PetscLogEvent DMPLEX_CoordinatesLoad;
     65: PETSC_EXTERN PetscLogEvent DMPLEX_SectionLoad;
     66: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalVectorLoad;
     67: PETSC_EXTERN PetscLogEvent DMPLEX_LocalVectorLoad;
     68: PETSC_EXTERN PetscLogEvent DMPLEX_MetricEnforceSPD;
     69: PETSC_EXTERN PetscLogEvent DMPLEX_MetricNormalize;
     70: PETSC_EXTERN PetscLogEvent DMPLEX_MetricAverage;
     71: PETSC_EXTERN PetscLogEvent DMPLEX_MetricIntersection;
     72: PETSC_EXTERN PetscLogEvent DMPLEX_Generate;
     73: PETSC_EXTERN PetscLogEvent DMPLEX_GetLocalOffsets;

     75: PETSC_EXTERN PetscLogEvent DMPLEX_RebalBuildGraph;
     76: PETSC_EXTERN PetscLogEvent DMPLEX_RebalRewriteSF;
     77: PETSC_EXTERN PetscLogEvent DMPLEX_RebalGatherGraph;
     78: PETSC_EXTERN PetscLogEvent DMPLEX_RebalPartition;
     79: PETSC_EXTERN PetscLogEvent DMPLEX_RebalScatterPart;
     80: PETSC_EXTERN PetscLogEvent DMPLEX_Uninterpolate;

     82: struct _n_PetscGridHash {
     83:   PetscInt     dim;
     84:   PetscReal    lower[3];    /* The lower-left corner */
     85:   PetscReal    upper[3];    /* The upper-right corner */
     86:   PetscReal    extent[3];   /* The box size */
     87:   PetscReal    h[3];        /* The subbox size */
     88:   PetscInt     n[3];        /* The number of subboxes */
     89:   PetscSection cellSection; /* Offsets for cells in each subbox*/
     90:   IS           cells;       /* List of cells in each subbox */
     91:   DMLabel      cellsSparse; /* Sparse storage for cell map */
     92: };

     94: typedef struct {
     95:   PetscBool isotropic;               /* Is the metric isotropic? */
     96:   PetscBool uniform;                 /* Is the metric uniform? */
     97:   PetscBool restrictAnisotropyFirst; /* Should anisotropy or normalization come first? */
     98:   PetscBool noInsert;                /* Should node insertion/deletion be turned off? */
     99:   PetscBool noSwap;                  /* Should facet swapping be turned off? */
    100:   PetscBool noMove;                  /* Should node movement be turned off? */
    101:   PetscBool noSurf;                  /* Should surface modification be turned off? */
    102:   PetscReal h_min, h_max;            /* Minimum/maximum tolerated metric magnitudes */
    103:   PetscReal a_max;                   /* Maximum tolerated anisotropy */
    104:   PetscReal targetComplexity;        /* Target metric complexity */
    105:   PetscReal p;                       /* Degree for L-p normalization methods */
    106:   PetscReal gradationFactor;         /* Maximum tolerated length ratio for adjacent edges */
    107:   PetscReal hausdorffNumber;         /* Max. distance between piecewise linear representation of boundary and reconstructed ideal boundary */
    108:   PetscInt  numIter;                 /* Number of ParMmg mesh adaptation iterations */
    109:   PetscInt  verbosity;               /* Level of verbosity for remesher (-1 = no output, 10 = maximum) */
    110: } DMPlexMetricCtx;

    112: /* Point Numbering in Plex:

    114:    Points are numbered contiguously by stratum. Strate are organized as follows:

    116:    First Stratum:  Cells [height 0]
    117:    Second Stratum: Vertices [depth 0]
    118:    Third Stratum:  Faces [height 1]
    119:    Fourth Stratum: Edges [depth 1]

    121:    We do this so that the numbering of a cell-vertex mesh does not change after interpolation. Within a given stratum,
    122:    we allow additional segregation of by cell type.
    123: */
    124: typedef struct {
    125:   PetscInt refct;

    127:   PetscSection coneSection;      /* Layout of cones (inedges for DAG) */
    128:   PetscInt    *cones;            /* Cone for each point */
    129:   PetscInt    *coneOrientations; /* Orientation of each cone point, means cone traversal should start on point 'o', and if negative start on -(o+1) and go in reverse */
    130:   PetscSection supportSection;   /* Layout of cones (inedges for DAG) */
    131:   PetscInt    *supports;         /* Cone for each point */

    133:   struct {                  // DMPolytopeType is an enum (usually size 4), but this needs frequent access
    134:     uint8_t value_as_uint8; // in a struct to guard for stronger typing
    135:   } *cellTypes;

    137:   /* Transformation */
    138:   DMPlexTransform tr;                                               /* Type of transform used to define an ephemeral mesh */
    139:   char           *transformType;                                    /* Type of transform for uniform cell refinement */
    140:   PetscBool       refinementUniform;                                /* Flag for uniform cell refinement */
    141:   PetscReal       refinementLimit;                                  /* Maximum volume for refined cell */
    142:   PetscErrorCode (*refinementFunc)(const PetscReal[], PetscReal *); /* Function giving the maximum volume for refined cell */

    144:   /* Interpolation */
    145:   DMPlexInterpolatedFlag interpolated;
    146:   DMPlexInterpolatedFlag interpolatedCollective;
    147:   PetscBool              interpolatePreferTensor; // When different orderings exist, prefer the tensor order

    149:   /* Ordering */
    150:   DMReorderDefaultFlag reorderDefault; /* Reorder the DM by default */

    152:   /* Distribution */
    153:   PetscBool distDefault;      /* Distribute the DM by default */
    154:   PetscInt  overlap;          /* Overlap of the partitions as passed to DMPlexDistribute() or DMPlexDistributeOverlap() */
    155:   PetscInt  numOvLabels;      /* The number of labels used for candidate overlap points */
    156:   DMLabel   ovLabels[16];     /* Labels used for candidate overlap points */
    157:   PetscInt  ovValues[16];     /* Label values used for candidate overlap points */
    158:   PetscInt  numOvExLabels;    /* The number of labels used for exclusion */
    159:   DMLabel   ovExLabels[16];   /* Labels used to exclude points from the overlap */
    160:   PetscInt  ovExValues[16];   /* Label values used to exclude points from the overlap */
    161:   char     *distributionName; /* Name of the specific parallel distribution of the DM */
    162:   MPI_Comm  nonempty_comm;    /* Communicator used for visualization when some processes do not have cells */

    164:   /* Hierarchy */
    165:   PetscBool regularRefinement; /* This flag signals that we are a regular refinement of coarseMesh */

    167:   /* Generation */
    168:   char            *tetgenOpts;
    169:   char            *triangleOpts;
    170:   PetscPartitioner partitioner;
    171:   PetscBool        partitionBalance; /* Evenly divide partition overlap when distributing */
    172:   PetscBool        remeshBd;

    174:   /* Submesh */
    175:   DMLabel          subpointMap;   /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */
    176:   IS               subpointIS;    /* IS holding point number in the enclosing mesh of every point in the submesh chart */
    177:   PetscObjectState subpointState; /* The state of subpointMap when the subpointIS was last created */

    179:   /* Labels and numbering */
    180:   PetscObjectState depthState;    /* State of depth label, so that we can determine if a user changes it */
    181:   PetscObjectState celltypeState; /* State of celltype label, so that we can determine if a user changes it */
    182:   IS               globalVertexNumbers;
    183:   IS               globalCellNumbers;

    185:   /* Constraints */
    186:   PetscSection anchorSection;          /* maps constrained points to anchor points */
    187:   IS           anchorIS;               /* anchors indexed by the above section */
    188:   PetscErrorCode (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
    189:   PetscErrorCode (*computeanchormatrix)(DM, PetscSection, PetscSection, Mat);

    191:   /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
    192:   PetscSection parentSection; /* dof == 1 if point has parent */
    193:   PetscInt    *parents;       /* point to parent */
    194:   PetscInt    *childIDs;      /* point to child ID */
    195:   PetscSection childSection;  /* inverse of parent section */
    196:   PetscInt    *children;      /* point to children */
    197:   DM           referenceTree; /* reference tree to which child ID's refer */
    198:   PetscErrorCode (*getchildsymmetry)(DM, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt *, PetscInt *);

    200:   /* MATIS support */
    201:   PetscSection subdomainSection;

    203:   /* Adjacency */
    204:   PetscBool useAnchors;                                                          /* Replace constrained points with their anchors in adjacency lists */
    205:   PetscErrorCode (*useradjacency)(DM, PetscInt, PetscInt *, PetscInt[], void *); /* User callback for adjacency */
    206:   void *useradjacencyctx;                                                        /* User context for callback */

    208:   // Periodicity
    209:   struct {
    210:     // Specified by the user
    211:     PetscInt num_face_sfs;          // number of face_sfs
    212:     PetscSF *face_sfs;              // root(donor faces) <-- leaf(local faces)
    213:     PetscScalar (*transform)[4][4]; // geometric transform
    214:     // Created eagerly (depends on points)
    215:     PetscSF composed_sf; // root(non-periodic global points) <-- leaf(local points)
    216:     IS     *periodic_points;
    217:   } periodic;

    219:   /* Projection */
    220:   PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
    221:   PetscInt activePoint;         /* current active point in iteration */

    223:   /* Output */
    224:   PetscInt  vtkCellHeight;          /* The height of cells for output, default is 0 */
    225:   PetscReal scale[NUM_PETSC_UNITS]; /* The scale for each SI unit */

    227:   /* Geometry */
    228:   PetscReal     minradius;       /* Minimum distance from cell centroid to face */
    229:   PetscBool     useHashLocation; /* Use grid hashing for point location */
    230:   PetscGridHash lbox;            /* Local box for searching */
    231:   PetscPointFn *coordFunc;       /* Function used to remap newly introduced vertices */

    233:   /* Neighbors */
    234:   PetscMPIInt *neighbors;

    236:   /* Metric */
    237:   DMPlexMetricCtx *metricCtx;

    239:   /* FEM */
    240:   PetscBool useCeed;      /* This should convert to a registration system when there are more FEM backends */
    241:   PetscBool useMatClPerm; /* Use the closure permutation when assembling matrices */

    243:   /* CAD */
    244:   PetscBool ignoreModel; /* If TRUE, Plex refinement will skip Snap-To-Geometry feature ignoring attached CAD geometry information */

    246:   // Transforms
    247:   PetscBool       saveTransform; // Flag to save transform, if this mesh was produced by a transform
    248:   DMPlexTransform transform;     // Transform producing this mesh (note that this will hold on to the original mesh too)

    250:   /* Debugging */
    251:   PetscBool printSetValues;
    252:   PetscInt  printAdj;
    253:   PetscInt  printFEM;
    254:   PetscInt  printFVM;
    255:   PetscInt  printL2;
    256:   PetscInt  printLocate;
    257:   PetscInt  printProject;
    258:   PetscReal printTol;
    259: } DM_Plex;

    261: PETSC_INTERN PetscErrorCode DMPlexCopy_Internal(DM, PetscBool, PetscBool, DM);
    262: PETSC_INTERN PetscErrorCode DMPlexReplace_Internal(DM, DM *);
    263: PETSC_INTERN PetscErrorCode DMPlexCopyEGADSInfo_Internal(DM, DM);

    265: PETSC_INTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM, PetscViewer);
    266: PETSC_INTERN PetscErrorCode VecView_Plex_Local(Vec, PetscViewer);
    267: PETSC_INTERN PetscErrorCode VecView_Plex_Native(Vec, PetscViewer);
    268: PETSC_INTERN PetscErrorCode VecView_Plex(Vec, PetscViewer);
    269: PETSC_INTERN PetscErrorCode VecLoad_Plex_Local(Vec, PetscViewer);
    270: PETSC_INTERN PetscErrorCode VecLoad_Plex_Native(Vec, PetscViewer);
    271: PETSC_INTERN PetscErrorCode VecLoad_Plex(Vec, PetscViewer);
    272: PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
    273: PETSC_INTERN PetscErrorCode DMPlexGetFieldTypes_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt **, PetscInt **, PetscViewerVTKFieldType **);
    274: PETSC_INTERN PetscErrorCode DMPlexRestoreFieldTypes_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt **, PetscInt **, PetscViewerVTKFieldType **);
    275: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM, PetscViewer);
    276: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject, PetscViewer);
    277: #if defined(PETSC_HAVE_HDF5)
    278: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
    279: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
    280: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
    281: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
    282: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
    283: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
    284: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
    285: PETSC_INTERN PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM, IS, PetscViewer);
    286: PETSC_INTERN PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM, PetscViewer);
    287: PETSC_INTERN PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM, IS, PetscViewer);
    288: PETSC_INTERN PetscErrorCode DMPlexSectionView_HDF5_Internal(DM, PetscViewer, DM);
    289: PETSC_INTERN PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM, PetscViewer, DM, Vec);
    290: PETSC_INTERN PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM, PetscViewer, DM, Vec);
    291: PETSC_INTERN PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM, PetscViewer, PetscSF *);
    292: PETSC_INTERN PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM, PetscViewer, PetscSF);
    293: PETSC_INTERN PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM, PetscViewer, PetscSF);
    294: PETSC_INTERN PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM, PetscViewer, DM, PetscSF, PetscSF *, PetscSF *);
    295: PETSC_INTERN PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM, PetscViewer, DM, PetscSF, Vec);
    296: PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
    297: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
    298: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM, PetscViewer);
    299: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
    300: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
    301: PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
    302: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
    303: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);
    304: #endif
    305: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_CGNS(Vec, PetscViewer);
    306: #if defined(PETSC_HAVE_CGNS)
    307: PETSC_EXTERN PetscErrorCode VecLoad_Plex_CGNS_Internal(Vec, PetscViewer);
    308: #endif

    310: PETSC_INTERN PetscErrorCode DMPlexClosurePoints_Private(DM, PetscInt, const PetscInt[], IS *);
    311: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(DM, PetscOptionItems);
    312: PETSC_INTERN PetscErrorCode DMSetFromOptions_Overlap_Plex(DM, PetscOptionItems, PetscInt *);
    313: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
    314: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM[]);
    315: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
    316: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM[]);
    317: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, Vec, DMLabel, DMLabel, DM *);
    318: PETSC_INTERN PetscErrorCode DMExtrude_Plex(DM, PetscInt, DM *);
    319: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
    320: PETSC_INTERN PetscErrorCode DMPlexInsertTimeDerivativeBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
    321: PETSC_INTERN PetscErrorCode DMPlexInsertBounds_Plex(DM, PetscBool, PetscReal, Vec);
    322: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, InsertMode, Vec);
    323: PETSC_INTERN PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM, PetscReal, DMLabel, PetscInt, const PetscInt[], PetscInt, const PetscInt[], PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, InsertMode, Vec);
    324: PETSC_INTERN PetscErrorCode DMProjectFieldLocal_Plex(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);
    325: PETSC_INTERN PetscErrorCode DMProjectFieldLabelLocal_Plex(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);
    326: PETSC_INTERN PetscErrorCode DMProjectBdFieldLabelLocal_Plex(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);
    327: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
    328: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, const PetscReal[], PetscReal *);
    329: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
    330: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);

    332: #if defined(PETSC_HAVE_EXODUSII)
    333: PETSC_INTERN PetscErrorCode DMView_PlexExodusII(DM, PetscViewer);
    334: PETSC_INTERN PetscErrorCode VecView_PlexExodusII_Internal(Vec, PetscViewer);
    335: PETSC_INTERN PetscErrorCode VecLoad_PlexExodusII_Internal(Vec, PetscViewer);
    336: #endif
    337: PETSC_INTERN PetscErrorCode DMView_PlexCGNS(DM, PetscViewer);
    338: PETSC_INTERN PetscErrorCode DMPlexCreateCGNSFromFile_Internal(MPI_Comm, const char *, PetscBool, DM *);
    339: PETSC_INTERN PetscErrorCode DMPlexCreateCGNS_Internal_Serial(MPI_Comm, PetscInt, PetscBool, DM *);
    340: PETSC_INTERN PetscErrorCode DMPlexCreateCGNS_Internal_Parallel(MPI_Comm, PetscInt, PetscBool, DM *);
    341: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM, PetscInt, PetscInt, PetscInt *);
    342: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM, PetscInt, PetscBool, PetscBool, PetscBool, PetscInt *, PetscInt *[]);
    343: PETSC_INTERN PetscErrorCode DMPlexGetMaxAdjacencySize_Internal(DM, PetscBool, PetscInt *);
    344: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM, DMPolytopeType, const PetscInt[], PetscInt *, const DMPolytopeType *[], const PetscInt *[], const PetscInt *[]);
    345: PETSC_INTERN PetscErrorCode DMPlexRestoreRawFaces_Internal(DM, DMPolytopeType, const PetscInt[], PetscInt *, const DMPolytopeType *[], const PetscInt *[], const PetscInt *[]);
    346: PETSC_INTERN PetscErrorCode DMPlexComputeCellType_Internal(DM, PetscInt, PetscInt, DMPolytopeType *);
    347: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], DMLabel, PetscInt, const PetscScalar[], InsertMode);
    348: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
    349: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM, PetscSection, PetscInt[], PetscInt[]);
    350: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM, DM, const char *, DM *);
    351: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM, DM, PetscSF, PetscInt *, Mat);
    352: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM, DM, PetscSF, PetscInt *, Mat);
    353: PETSC_INTERN PetscErrorCode DMPlexAnchorsModifyMat(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, const PetscScalar[], PetscInt *, PetscInt *, PetscInt *[], PetscScalar *[], PetscInt[], PetscBool);
    354: PETSC_INTERN PetscErrorCode DMPlexAnchorsModifyMat_Internal(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, PetscInt, PetscInt, const PetscScalar[], PetscInt *, PetscInt *, PetscInt *[], PetscScalar *[], PetscInt[], PetscBool, PetscBool);
    355: PETSC_INTERN PetscErrorCode DMPlexAnchorsGetSubMatModification(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, PetscInt *, PetscInt *, PetscInt *[], PetscInt[], PetscScalar *[]);
    356: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM, PetscInt, const PetscScalar[], PetscInt, PetscInt *);
    357: /* this is PETSC_EXTERN just because of src/dm/impls/plex/tests/ex18.c */
    358: PETSC_EXTERN PetscErrorCode DMPlexOrientInterface_Internal(DM);
    359: PETSC_INTERN PetscErrorCode DMPlexOrientCells_Internal(DM, IS, IS);

    361: /* Applications may use this function */
    362: PETSC_EXTERN PetscErrorCode DMPlexCreateNumbering_Plex(DM, PetscInt, PetscInt, PetscInt, PetscInt *, PetscSF, IS *);

    364: PETSC_INTERN PetscErrorCode DMPlexInterpolateInPlace_Internal(DM);
    365: PETSC_INTERN PetscErrorCode DMPlexCreateBoxMesh_Tensor_SFC_Internal(DM, PetscInt, const PetscInt[], const PetscReal[], const PetscReal[], const DMBoundaryType[], PetscBool);
    366: PETSC_INTERN PetscErrorCode DMPlexMigrateIsoperiodicFaceSF_Internal(DM, DM, PetscSF);
    367: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
    368: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, Vec, DMLabel, DMLabel, DM *);
    369: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, Vec, DMLabel, DMLabel, DM *);
    370: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM, Mat *);

    372: PETSC_INTERN PetscErrorCode DMPlexGetOverlap_Plex(DM, PetscInt *);
    373: PETSC_INTERN PetscErrorCode DMPlexSetOverlap_Plex(DM, DM, PetscInt);
    374: PETSC_INTERN PetscErrorCode DMPlexDistributeGetDefault_Plex(DM, PetscBool *);
    375: PETSC_INTERN PetscErrorCode DMPlexDistributeSetDefault_Plex(DM, PetscBool);
    376: PETSC_INTERN PetscErrorCode DMPlexReorderGetDefault_Plex(DM, DMReorderDefaultFlag *);
    377: PETSC_INTERN PetscErrorCode DMPlexReorderSetDefault_Plex(DM, DMReorderDefaultFlag);
    378: PETSC_INTERN PetscErrorCode DMPlexGetUseCeed_Plex(DM, PetscBool *);
    379: PETSC_INTERN PetscErrorCode DMPlexSetUseCeed_Plex(DM, PetscBool);
    380: PETSC_INTERN PetscErrorCode DMReorderSectionGetDefault_Plex(DM, DMReorderDefaultFlag *);
    381: PETSC_INTERN PetscErrorCode DMReorderSectionSetDefault_Plex(DM, DMReorderDefaultFlag);
    382: PETSC_INTERN PetscErrorCode DMReorderSectionGetType_Plex(DM, MatOrderingType *);
    383: PETSC_INTERN PetscErrorCode DMReorderSectionSetType_Plex(DM, MatOrderingType);

    385: #if 1
    386: static inline PetscInt DihedralInvert(PetscInt N, PetscInt a)
    387: {
    388:   return (a <= 0) ? a : (N - a);
    389: }

    391: static inline PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
    392: {
    393:   if (!N) return 0;
    394:   return (a >= 0) ? ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) : ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
    395: }

    397: static inline PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
    398: {
    399:   return DihedralCompose(N, DihedralInvert(N, a), b);
    400: }
    401: #else
    402: /* TODO
    403:    This is a reimplementation of the tensor dihedral symmetries using the new orientations.
    404:    These should be turned on when we convert to new-style orientations in p4est.
    405: */
    406: /* invert dihedral symmetry: return a^-1,
    407:  * using the representation described in
    408:  * DMPlexGetConeOrientation() */
    409: static inline PetscInt DihedralInvert(PetscInt N, PetscInt a)
    410: {
    411:   switch (N) {
    412:   case 0:
    413:     return 0;
    414:   case 2:
    415:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, 0, a);
    416:   case 4:
    417:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, 0, a);
    418:   case 8:
    419:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, 0, a);
    420:   default:
    421:     SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralInvert()");
    422:   }
    423:   return 0;
    424: }

    426: /* compose dihedral symmetry: return b * a,
    427:  * using the representation described in
    428:  * DMPlexGetConeOrientation() */
    429: static inline PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
    430: {
    431:   switch (N) {
    432:   case 0:
    433:     return 0;
    434:   case 2:
    435:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
    436:   case 4:
    437:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
    438:   case 8:
    439:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
    440:   default:
    441:     SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
    442:   }
    443:   return 0;
    444: }

    446: /* swap dihedral symmetries: return b * a^-1,
    447:  * using the representation described in
    448:  * DMPlexGetConeOrientation() */
    449: static inline PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
    450: {
    451:   switch (N) {
    452:   case 0:
    453:     return 0;
    454:   case 2:
    455:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
    456:   case 4:
    457:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
    458:   case 8:
    459:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
    460:   default:
    461:     SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
    462:   }
    463:   return 0;
    464: }
    465: #endif
    466: PETSC_INTERN PetscInt       DMPolytopeConvertNewOrientation_Internal(DMPolytopeType, PetscInt);
    467: PETSC_INTERN PetscInt       DMPolytopeConvertOldOrientation_Internal(DMPolytopeType, PetscInt);
    468: PETSC_INTERN PetscErrorCode DMPlexConvertOldOrientations_Internal(DM);

    470: PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode DMPlexComputeIntegral_Internal(DM, Vec, PetscInt, PetscInt, PetscScalar *, void *);
    471: PETSC_INTERN PetscErrorCode                             DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);

    473: /* Matvec with A in row-major storage, x and y can be aliased */
    474: static inline void DMPlex_Mult2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
    475: {
    476:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
    477:   y[0 * ldx]             = A[0] * z[0] + A[1] * z[1];
    478:   y[1 * ldx]             = A[2] * z[0] + A[3] * z[1];
    479:   (void)PetscLogFlops(6.0);
    480: }
    481: static inline void DMPlex_Mult3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
    482: {
    483:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
    484:   y[0 * ldx]             = A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
    485:   y[1 * ldx]             = A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
    486:   y[2 * ldx]             = A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
    487:   (void)PetscLogFlops(15.0);
    488: }
    489: static inline void DMPlex_MultTranspose2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
    490: {
    491:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
    492:   y[0 * ldx]             = A[0] * z[0] + A[2] * z[1];
    493:   y[1 * ldx]             = A[1] * z[0] + A[3] * z[1];
    494:   (void)PetscLogFlops(6.0);
    495: }
    496: static inline void DMPlex_MultTranspose3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
    497: {
    498:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
    499:   y[0 * ldx]             = PetscConj(A[0]) * z[0] + PetscConj(A[3]) * z[1] + PetscConj(A[6]) * z[2];
    500:   y[1 * ldx]             = PetscConj(A[1]) * z[0] + PetscConj(A[4]) * z[1] + PetscConj(A[7]) * z[2];
    501:   y[2 * ldx]             = PetscConj(A[2]) * z[0] + PetscConj(A[5]) * z[1] + PetscConj(A[8]) * z[2];
    502:   (void)PetscLogFlops(15.0);
    503: }
    504: static inline void DMPlex_Mult2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
    505: {
    506:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
    507:   y[0 * ldx]             = A[0] * z[0] + A[1] * z[1];
    508:   y[1 * ldx]             = A[2] * z[0] + A[3] * z[1];
    509:   (void)PetscLogFlops(6.0);
    510: }
    511: static inline void DMPlex_Mult3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
    512: {
    513:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
    514:   y[0 * ldx]             = A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
    515:   y[1 * ldx]             = A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
    516:   y[2 * ldx]             = A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
    517:   (void)PetscLogFlops(15.0);
    518: }
    519: static inline void DMPlex_MultAdd2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
    520: {
    521:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
    522:   y[0 * ldx] += A[0] * z[0] + A[1] * z[1];
    523:   y[1 * ldx] += A[2] * z[0] + A[3] * z[1];
    524:   (void)PetscLogFlops(6.0);
    525: }
    526: static inline void DMPlex_MultAdd3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
    527: {
    528:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
    529:   y[0 * ldx] += A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
    530:   y[1 * ldx] += A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
    531:   y[2 * ldx] += A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
    532:   (void)PetscLogFlops(15.0);
    533: }
    534: /*
    535:   A: packed, row-major m x n array
    536:   x: length m array
    537:   y: length n arra
    538:   ldx: the stride in x and y

    540:   A[i,j] = A[i * n + j]
    541:   A^T[j,i] = A[i,j]
    542: */
    543: static inline void DMPlex_MultTransposeReal_Internal(const PetscReal A[], PetscInt m, PetscInt n, PetscInt ldx, const PetscScalar x[], PetscScalar y[])
    544: {
    545:   PetscScalar z[3];
    546:   PetscInt    i, j;
    547:   for (i = 0; i < m; ++i) z[i] = x[i * ldx];
    548:   for (j = 0; j < n; ++j) {
    549:     const PetscInt l = j * ldx;
    550:     y[l]             = 0;
    551:     for (i = 0; i < m; ++i) y[l] += A[i * n + j] * z[i];
    552:   }
    553:   (void)PetscLogFlops(2 * m * n);
    554: }
    555: static inline void DMPlex_MultTranspose2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
    556: {
    557:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
    558:   y[0 * ldx]             = A[0] * z[0] + A[2] * z[1];
    559:   y[1 * ldx]             = A[1] * z[0] + A[3] * z[1];
    560:   (void)PetscLogFlops(6.0);
    561: }
    562: static inline void DMPlex_MultTranspose3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
    563: {
    564:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
    565:   y[0 * ldx]             = A[0] * z[0] + A[3] * z[1] + A[6] * z[2];
    566:   y[1 * ldx]             = A[1] * z[0] + A[4] * z[1] + A[7] * z[2];
    567:   y[2 * ldx]             = A[2] * z[0] + A[5] * z[1] + A[8] * z[2];
    568:   (void)PetscLogFlops(15.0);
    569: }

    571: static inline void DMPlex_MatMult2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
    572: {
    573: #define PLEX_DIM__ 2
    574:   PetscScalar z[PLEX_DIM__];
    575:   for (PetscInt j = 0; j < n; ++j) {
    576:     for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
    577:     DMPlex_Mult2D_Internal(A, 1, z, z);
    578:     for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
    579:   }
    580:   (void)PetscLogFlops(8.0 * n);
    581: #undef PLEX_DIM__
    582: }
    583: static inline void DMPlex_MatMult3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
    584: {
    585: #define PLEX_DIM__ 3
    586:   PetscScalar z[PLEX_DIM__];
    587:   for (PetscInt j = 0; j < n; ++j) {
    588:     for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
    589:     DMPlex_Mult3D_Internal(A, 1, z, z);
    590:     for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
    591:   }
    592:   (void)PetscLogFlops(8.0 * n);
    593: #undef PLEX_DIM__
    594: }
    595: static inline void DMPlex_MatMultTranspose2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
    596: {
    597: #define PLEX_DIM__ 2
    598:   PetscScalar z[PLEX_DIM__];
    599:   for (PetscInt j = 0; j < n; ++j) {
    600:     for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
    601:     DMPlex_MultTranspose2D_Internal(A, 1, z, z);
    602:     for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
    603:   }
    604:   (void)PetscLogFlops(8.0 * n);
    605: #undef PLEX_DIM__
    606: }
    607: static inline void DMPlex_MatMultTranspose3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
    608: {
    609: #define PLEX_DIM__ 3
    610:   PetscScalar z[PLEX_DIM__];
    611:   for (PetscInt j = 0; j < n; ++j) {
    612:     for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
    613:     DMPlex_MultTranspose3D_Internal(A, 1, z, z);
    614:     for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
    615:   }
    616:   (void)PetscLogFlops(8.0 * n);
    617: #undef PLEX_DIM__
    618: }

    620: static inline void DMPlex_MatMultLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
    621: {
    622:   for (PetscInt j = 0; j < m; ++j) DMPlex_MultTranspose2D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
    623:   (void)PetscLogFlops(8.0 * m);
    624: }
    625: static inline void DMPlex_MatMultLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
    626: {
    627:   for (PetscInt j = 0; j < m; ++j) DMPlex_MultTranspose3D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
    628:   (void)PetscLogFlops(8.0 * m);
    629: }
    630: static inline void DMPlex_MatMultTransposeLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
    631: {
    632:   for (PetscInt j = 0; j < m; ++j) DMPlex_Mult2D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
    633:   (void)PetscLogFlops(8.0 * m);
    634: }
    635: static inline void DMPlex_MatMultTransposeLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
    636: {
    637:   for (PetscInt j = 0; j < m; ++j) DMPlex_Mult3D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
    638:   (void)PetscLogFlops(8.0 * m);
    639: }

    641: static inline void DMPlex_PTAP2DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
    642: {
    643:   PetscScalar out[4];
    644:   PetscInt    i, j, k, l;
    645:   for (i = 0; i < 2; ++i) {
    646:     for (j = 0; j < 2; ++j) {
    647:       out[i * 2 + j] = 0.;
    648:       for (k = 0; k < 2; ++k) {
    649:         for (l = 0; l < 2; ++l) out[i * 2 + j] += P[k * 2 + i] * A[k * 2 + l] * P[l * 2 + j];
    650:       }
    651:     }
    652:   }
    653:   for (i = 0; i < 2 * 2; ++i) C[i] = out[i];
    654:   (void)PetscLogFlops(48.0);
    655: }
    656: static inline void DMPlex_PTAP3DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
    657: {
    658:   PetscScalar out[9];
    659:   PetscInt    i, j, k, l;
    660:   for (i = 0; i < 3; ++i) {
    661:     for (j = 0; j < 3; ++j) {
    662:       out[i * 3 + j] = 0.;
    663:       for (k = 0; k < 3; ++k) {
    664:         for (l = 0; l < 3; ++l) out[i * 3 + j] += P[k * 3 + i] * A[k * 3 + l] * P[l * 3 + j];
    665:       }
    666:     }
    667:   }
    668:   for (i = 0; i < 2 * 2; ++i) C[i] = out[i];
    669:   (void)PetscLogFlops(243.0);
    670: }
    671: /* TODO Fix for aliasing of A and C */
    672: static inline void DMPlex_PTAPReal_Internal(const PetscReal P[], PetscInt m, PetscInt n, const PetscScalar A[], PetscScalar C[])
    673: {
    674:   PetscInt i, j, k, l;
    675:   for (i = 0; i < n; ++i) {
    676:     for (j = 0; j < n; ++j) {
    677:       C[i * n + j] = 0.;
    678:       for (k = 0; k < m; ++k) {
    679:         for (l = 0; l < m; ++l) C[i * n + j] += P[k * n + i] * A[k * m + l] * P[l * n + j];
    680:       }
    681:     }
    682:   }
    683:   (void)PetscLogFlops(243.0);
    684: }

    686: static inline void DMPlex_Transpose2D_Internal(PetscScalar A[])
    687: {
    688:   PetscScalar tmp;
    689:   tmp  = A[1];
    690:   A[1] = A[2];
    691:   A[2] = tmp;
    692: }
    693: static inline void DMPlex_Transpose3D_Internal(PetscScalar A[])
    694: {
    695:   PetscScalar tmp;
    696:   tmp  = A[1];
    697:   A[1] = A[3];
    698:   A[3] = tmp;
    699:   tmp  = A[2];
    700:   A[2] = A[6];
    701:   A[6] = tmp;
    702:   tmp  = A[5];
    703:   A[5] = A[7];
    704:   A[7] = tmp;
    705: }

    707: static inline void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
    708: {
    709:   // Allow zero volume cells
    710:   const PetscReal invDet = detJ == 0 ? 1.0 : (PetscReal)1.0 / detJ;

    712:   invJ[0] = invDet * J[3];
    713:   invJ[1] = -invDet * J[1];
    714:   invJ[2] = -invDet * J[2];
    715:   invJ[3] = invDet * J[0];
    716:   (void)PetscLogFlops(5.0);
    717: }

    719: static inline void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
    720: {
    721:   // Allow zero volume cells
    722:   const PetscReal invDet = detJ == 0 ? 1.0 : (PetscReal)1.0 / detJ;

    724:   invJ[0 * 3 + 0] = invDet * (J[1 * 3 + 1] * J[2 * 3 + 2] - J[1 * 3 + 2] * J[2 * 3 + 1]);
    725:   invJ[0 * 3 + 1] = invDet * (J[0 * 3 + 2] * J[2 * 3 + 1] - J[0 * 3 + 1] * J[2 * 3 + 2]);
    726:   invJ[0 * 3 + 2] = invDet * (J[0 * 3 + 1] * J[1 * 3 + 2] - J[0 * 3 + 2] * J[1 * 3 + 1]);
    727:   invJ[1 * 3 + 0] = invDet * (J[1 * 3 + 2] * J[2 * 3 + 0] - J[1 * 3 + 0] * J[2 * 3 + 2]);
    728:   invJ[1 * 3 + 1] = invDet * (J[0 * 3 + 0] * J[2 * 3 + 2] - J[0 * 3 + 2] * J[2 * 3 + 0]);
    729:   invJ[1 * 3 + 2] = invDet * (J[0 * 3 + 2] * J[1 * 3 + 0] - J[0 * 3 + 0] * J[1 * 3 + 2]);
    730:   invJ[2 * 3 + 0] = invDet * (J[1 * 3 + 0] * J[2 * 3 + 1] - J[1 * 3 + 1] * J[2 * 3 + 0]);
    731:   invJ[2 * 3 + 1] = invDet * (J[0 * 3 + 1] * J[2 * 3 + 0] - J[0 * 3 + 0] * J[2 * 3 + 1]);
    732:   invJ[2 * 3 + 2] = invDet * (J[0 * 3 + 0] * J[1 * 3 + 1] - J[0 * 3 + 1] * J[1 * 3 + 0]);
    733:   (void)PetscLogFlops(37.0);
    734: }

    736: static inline void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
    737: {
    738:   *detJ = J[0] * J[3] - J[1] * J[2];
    739:   (void)PetscLogFlops(3.0);
    740: }

    742: static inline void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
    743: {
    744:   *detJ = (J[0 * 3 + 0] * (J[1 * 3 + 1] * J[2 * 3 + 2] - J[1 * 3 + 2] * J[2 * 3 + 1]) + J[0 * 3 + 1] * (J[1 * 3 + 2] * J[2 * 3 + 0] - J[1 * 3 + 0] * J[2 * 3 + 2]) + J[0 * 3 + 2] * (J[1 * 3 + 0] * J[2 * 3 + 1] - J[1 * 3 + 1] * J[2 * 3 + 0]));
    745:   (void)PetscLogFlops(12.0);
    746: }

    748: static inline void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
    749: {
    750:   *detJ = PetscRealPart(J[0]) * PetscRealPart(J[3]) - PetscRealPart(J[1]) * PetscRealPart(J[2]);
    751:   (void)PetscLogFlops(3.0);
    752: }

    754: static inline void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
    755: {
    756:   *detJ = (PetscRealPart(J[0 * 3 + 0]) * (PetscRealPart(J[1 * 3 + 1]) * PetscRealPart(J[2 * 3 + 2]) - PetscRealPart(J[1 * 3 + 2]) * PetscRealPart(J[2 * 3 + 1])) + PetscRealPart(J[0 * 3 + 1]) * (PetscRealPart(J[1 * 3 + 2]) * PetscRealPart(J[2 * 3 + 0]) - PetscRealPart(J[1 * 3 + 0]) * PetscRealPart(J[2 * 3 + 2])) + PetscRealPart(J[0 * 3 + 2]) * (PetscRealPart(J[1 * 3 + 0]) * PetscRealPart(J[2 * 3 + 1]) - PetscRealPart(J[1 * 3 + 1]) * PetscRealPart(J[2 * 3 + 0])));
    757:   (void)PetscLogFlops(12.0);
    758: }

    760: static inline void DMPlex_WaxpyD_Internal(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w)
    761: {
    762:   PetscInt d;
    763:   for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
    764: }

    766: static inline PetscReal DMPlex_DotD_Internal(PetscInt dim, const PetscScalar *x, const PetscReal *y)
    767: {
    768:   PetscReal sum = 0.0;
    769:   PetscInt  d;
    770:   for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d]) * y[d];
    771:   return sum;
    772: }

    774: static inline PetscReal DMPlex_DotRealD_Internal(PetscInt dim, const PetscReal *x, const PetscReal *y)
    775: {
    776:   PetscReal sum = 0.0;
    777:   PetscInt  d;
    778:   for (d = 0; d < dim; ++d) sum += x[d] * y[d];
    779:   return sum;
    780: }

    782: static inline PetscReal DMPlex_NormD_Internal(PetscInt dim, const PetscReal *x)
    783: {
    784:   PetscReal sum = 0.0;
    785:   PetscInt  d;
    786:   for (d = 0; d < dim; ++d) sum += x[d] * x[d];
    787:   return PetscSqrtReal(sum);
    788: }

    790: static inline PetscReal DMPlex_DistD_Internal(PetscInt dim, const PetscScalar *x, const PetscScalar *y)
    791: {
    792:   PetscReal sum = 0.0;
    793:   PetscInt  d;
    794:   for (d = 0; d < dim; ++d) sum += PetscRealPart(PetscConj(x[d] - y[d]) * (x[d] - y[d]));
    795:   return PetscSqrtReal(sum);
    796: }

    798: static inline PetscReal DMPlex_DistRealD_Internal(PetscInt dim, const PetscReal *x, const PetscReal *y)
    799: {
    800:   PetscReal sum = 0.0;
    801:   PetscInt  d;
    802:   for (d = 0; d < dim; ++d) sum += (x[d] - y[d]) * (x[d] - y[d]);
    803:   return PetscSqrtReal(sum);
    804: }

    806: PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM, PetscInt, PetscInt, PetscDualSpace *);
    807: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection, PetscBool, PetscInt, PetscInt, PetscInt *, PetscBool, const PetscInt[], const PetscInt[], PetscInt[]);
    808: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection, PetscBool, PetscInt, PetscInt, PetscInt[], PetscBool, const PetscInt ***, PetscInt, const PetscInt[], PetscInt[]);
    809: PETSC_INTERN PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM, PetscInt, PetscInt, PetscBool, PetscInt *, PetscInt *[]);
    810: PETSC_INTERN PetscErrorCode DMPlexMatSetClosure_Internal(DM, PetscSection, PetscSection, PetscBool, Mat, PetscInt, const PetscScalar[], InsertMode);

    812: PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode DMPlexGetAllCells_Internal(DM, IS *);
    813: PETSC_INTERN PetscErrorCode                             DMPlexGetAllFaces_Internal(DM, IS *);
    814: PETSC_EXTERN PetscErrorCode                             DMSNESGetFEGeom(DMField, IS, PetscQuadrature, PetscFEGeomMode, PetscFEGeom **);
    815: PETSC_EXTERN PetscErrorCode                             DMSNESRestoreFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
    816: PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM, PetscSection, IS, PetscReal, Vec, Vec, Vec, void *);
    817: PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM, PetscSection, PetscSection, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
    818: PETSC_INTERN PetscErrorCode                             DMCreateSubDomainDM_Plex(DM, DMLabel, PetscInt, IS *, DM *);
    819: PETSC_INTERN PetscErrorCode                             DMPlexBasisTransformPoint_Internal(DM, DM, Vec, PetscInt, PetscBool[], PetscBool, PetscScalar *);
    820: PETSC_INTERN PetscErrorCode                             DMPlexBasisTransformPointTensor_Internal(DM, DM, Vec, PetscInt, PetscBool, PetscInt, PetscScalar *);
    821: PETSC_INTERN PetscErrorCode                             DMPlexBasisTransformApplyReal_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscReal *, PetscReal *, void *);
    822: PETSC_INTERN PetscErrorCode                             DMPlexBasisTransformApply_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscScalar *, PetscScalar *, void *);
    823: PETSC_INTERN PetscErrorCode                             DMCreateNeumannOverlap_Plex(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **);
    824: PETSC_INTERN PetscErrorCode                             DMPlexMarkBoundaryFaces_Internal(DM, PetscInt, PetscInt, DMLabel, PetscBool);
    825: PETSC_INTERN PetscErrorCode                             DMPlexDistributeOverlap_Internal(DM, PetscInt, MPI_Comm, const char *, PetscSF *, DM *);

    827: PETSC_INTERN PetscErrorCode DMPlexInterpolateFaces_Internal(DM, PetscInt, DM);

    829: PETSC_INTERN PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM, DMLabel, PetscInt, PetscBool, PetscBool, DMLabel, DM);

    831: PETSC_INTERN PetscErrorCode DMPeriodicCoordinateSetUp_Internal(DM);

    833: /* Functions in the vtable */
    834: PETSC_INTERN PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
    835: PETSC_INTERN PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat);
    836: PETSC_INTERN PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat);
    837: PETSC_INTERN PetscErrorCode DMCreateMassMatrixLumped_Plex(DM, Vec *, Vec *);
    838: PETSC_INTERN PetscErrorCode DMCreateGradientMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat);
    839: PETSC_INTERN PetscErrorCode DMCreateLocalSection_Plex(DM dm);
    840: PETSC_INTERN PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm);
    841: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J);
    842: PETSC_INTERN PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
    843: PETSC_INTERN PetscErrorCode DMCreateCellCoordinateDM_Plex(DM dm, DM *cdm);
    844: PETSC_INTERN PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field);
    845: PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm);
    846: PETSC_INTERN PetscErrorCode DMSetUp_Plex(DM dm);
    847: PETSC_INTERN PetscErrorCode DMDestroy_Plex(DM dm);
    848: PETSC_INTERN PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
    849: PETSC_INTERN PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
    850: PETSC_INTERN PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm);
    851: PETSC_INTERN PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm);
    852: PETSC_INTERN PetscErrorCode DMCreateDomainDecompositionScatters_Plex(DM, PetscInt, DM *, VecScatter **, VecScatter **, VecScatter **);
    853: PETSC_INTERN PetscErrorCode DMCreateDomainDecomposition_Plex(DM, PetscInt *, char ***, IS **, IS **, DM **);
    854: PETSC_INTERN PetscErrorCode DMCreateSectionPermutation_Plex(DM dm, IS *permutation, PetscBT *blockStarts);

    856: // Coordinate mapping functions
    857: PETSC_INTERN void coordMap_identity(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[]);
    858: PETSC_INTERN void coordMap_shear(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[]);
    859: PETSC_INTERN void coordMap_flare(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[]);
    860: PETSC_INTERN void coordMap_annulus(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[]);
    861: PETSC_INTERN void coordMap_shell(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[]);
    862: PETSC_INTERN void coordMap_sinusoid(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[]);

    864: PETSC_INTERN PetscErrorCode DMSnapToGeomModel_EGADS(DM, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);

    866: // FIXME: DM with Attached CAD Models (STEP, IGES, BRep, EGADS, EGADSlite)
    867: PETSC_EXTERN PetscErrorCode DMPlexSnapToGeomModel(DM, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);

    869: // Coordinate <-> Reference mapping functions
    870: PETSC_INTERN PetscErrorCode DMPlexCoordinatesToReference_FE(DM, PetscFE, PetscInt, PetscInt, const PetscReal[], PetscReal[], Vec, PetscInt, PetscInt, PetscInt, PetscReal *);
    871: PETSC_INTERN PetscErrorCode DMPlexReferenceToCoordinates_FE(DM, PetscFE, PetscInt, PetscInt, const PetscReal[], PetscReal[], Vec, PetscInt, PetscInt);