wmudggeometry.h

Go to the documentation of this file.

#ifndef WMUDGGEOMETRY_H
#define WMUDGGEOMETRY_H
 
// Wm includes
#include "basis/wmbasisarrayset.h"
#include "lib/geometry/unstructured_geometry/wmunstructuredgeometry.h"
#include "lib/wmunstructuredpatch.h"
 
// STL includes
#include <memory>
 
class face_t
{
    friend class WmUDGGeometry;
 
public:
    int get_inside_element_index() const
    {
        return _insideElementIndex;
    }
    int get_outside_element_index() const
    {
        return _outsideElementIndex;
    }
    int get_inside_face_index() const
    {
        return _insideFaceIndex;
    }
    int get_outside_face_index() const
    {
        return _outsideFaceIndex;
    }
    int get_inside_face_orientation() const
    {
        return _insideFaceOrientation;
    }
    int get_outside_face_orientation() const
    {
        return _outsideFaceOrientation;
    }
    const real* get_norm() const
    {
        return _norm;
    }
    int get_global_id() const
    {
        return _globalIndex;
    }
    const int* get_inside_face_global_node_indexes() const
    {
        return _insideFaceGlobalNodeIndexes;
    }
    const int* get_outside_face_global_node_indexes() const
    {
        return _outsideFaceGlobalNodeIndexes;
    }
    const int* get_inside_face_local_node_indexes() const
    {
        return _insideFaceLocalNodeIndexes;
    }
    const int* get_outside_face_local_node_indexes() const
    {
        return _outsideFaceLocalNodeIndexes;
    }
    bool is_boundary_face() const
    {
        return _isBoundaryFace;
    }
    bool is_patch_boundary_face() const
    {
        return _isPatchBoundaryFace;
    }
 
    void print(std::stringstream& ss) const;
 
protected:
    face_t(const int insideElementIndex,
           const int outsideElementIndex,
           const int insideFaceIndex,
           const int outsideFaceIndex,
           const WmUnstructuredGeometry& geometry,
           const WmBasisArraySet& dgBasisSet);
 
    void setup(const WmBasisArraySet& dgBasisSet);
 
    int _insideElementIndex;
    int _outsideElementIndex;
    int _insideFaceIndex;
    int _outsideFaceIndex;
    int _insideFaceOrientation;
    int _outsideFaceOrientation;
    const real* _norm;
    int _globalIndex;
 
    int* _insideFaceGlobalNodeIndexes;
    int* _outsideFaceGlobalNodeIndexes;
    int* _insideFaceLocalNodeIndexes;
    int* _outsideFaceLocalNodeIndexes;
 
    bool _isBoundaryFace;
    bool _isPatchBoundaryFace;
 
private:
    std::vector<int> _nodeIndexAllocation;
};
 
class WmUDGGeometry
{
public:
    typedef std::map<std::string, std::vector<int>> boundaryFacesMap_t;
    typedef std::map<std::string, std::vector<std::unique_ptr<face_t>>> bcFaceMap_t;
 
    WmUDGGeometry(const WmUnstructuredPatch& patch,
                  const std::vector<std::string>& boundaryNames,
                  const WmBasisArraySet& dgBasisSet,
                  const bool useLimiter = false);
 
    ~WmUDGGeometry();
 
    const boundaryFacesMap_t& getBoundaryFaces() const
    {
        return _boundaryFaceIndexes;
    }
 
    int getNumFaces() const
    {
        return _numFaces;
    }
    int getNumFaces_Interior() const
    {
        return _numInteriorFaces;
    }
    int getNumFaces_Boundary() const
    {
        return _numBoundaryFaces;
    }
    int getNumFaces_Shared() const
    {
        return _numSharedFaces;
    }
 
    int getNumElements() const
    {
        return _numElements_PEG;
    }
    int getNumElements_Inclusive() const
    {
        return _numElements_PIPBCEG;
    }
 
    int getNumElements_Interpatch_Elements_Group() const
    {
        return _numElements_PIPEG; // right now only giving up to 1 layer?
    }
 
    const real* getGeometryArray_Face() const
    {
        return &(_faceGeometry[0]);
    }
 
    const real* getGeometryArray_Element() const
    {
        return &(_elementGeometry[0]);
    }
    const real* getGeometryArray_Jacobian() const
    {
        return &(_jacobianGeometry[0]);
    }
    const real* getGeometryArray_Limiter() const
    {
        return &(_limiterGeometry[0]);
    }
 
    int getGeometryArray_Face_BlockSize() const
    {
        return _faceGeometryBlockSize;
    }
 
    int getGeometryArray_Element_BlockSize() const
    {
        return _elementGeometryBlockSize;
    }
 
    int getGeometryArray_Jacobian_BlockSize() const
    {
        return _jacobianGeometryBlockSize;
    }
 
    int getGeometryArray_Limiter_BlockSize() const
    {
        return _limiterGeometryBlockSize;
    }
 
    int getGeometryArray_Face_Size() const
    {
        return _numFaces * _faceGeometryBlockSize;
    }
 
    int getGeometryArray_Element_Size() const
    {
        return _numElements_PIPEG * _elementGeometryBlockSize;
    }
 
    int getGeometryArray_Jacobian_Size() const
    {
        return _numElements_PIPEG * _jacobianGeometryBlockSize;
    }
 
    int getGeometryArray_Limiter_Size() const
    {
        return _numElements_PIPEG * _limiterGeometryBlockSize;
    }
 
    const int* getReconstructionArray_StencilIndexes() const
    {
        return &(_reconstruction_stencilIndexes[0]);
    }
 
    const real* getReconstructionArray_Weights() const
    {
        return &(_reconstruction_w[0]);
    }
 
    const real* getReconstructionArray_Reconstruction() const
    {
        return &(_reconstruction_lInverse[0]);
    }
 
    const real* getReconstructionArray_Correction() const
    {
        return &(_reconstruction_b[0]);
    }
 
    int getReconstructionArray_StencilIndexes_BlockSize() const
    {
        return _reconstruction_stencilIndexes_blockSize;
    }
 
    int getReconstructionArray_Weights_BlockSize() const
    {
        return _reconstruction_w_blockSize;
    }
 
    int getReconstructionArray_Reconstruction_BlockSize() const
    {
        return _reconstruction_lInverse_blockSize;
    }
 
    int getReconstructionArray_Correction_BlockSize() const
    {
        return _reconstruction_b_blockSize;
    }
 
    int getReconstructionArray_StencilIndexes_Size() const
    {
        return _reconstruction_stencilIndexes.size();
    }
 
    int getReconstructionArray_Weights_Size() const
    {
        return _reconstruction_w.size();
    }
 
    int getReconstructionArray_Reconstruction_Size() const
    {
        return _reconstruction_lInverse.size();
    }
 
    int getReconstructionArray_Correction_Size() const
    {
        return _reconstruction_b.size();
    }
 
    const int* getFaceNodeIndexes_Interconnected() const
    {
        return &(_interconnectedFaceNodeIndexes[0]);
    }
 
    // returns interconnected elements by face (inside,outside)
    const int* getFaceElementIndexes_Interconnected() const
    {
        return &(_interconnectedFaceElementIndexes[0]);
    }
 
    // returns interconnected elements by face (inside,outside)
    // but outside cannot be -1 -> instead increases numbering after _numElements_PIPEG
    const int* getFaceElementIndexes_InterconnectedNoMinusOne() const
    {
        return &(_interconnectedFaceElementIndexesNoMinusOne[0]);
    }
 
    // this yields local node numbers for global faces
    const int* getFaceNodeIndexes_Composite() const
    {
        return &(_compositeFaceNodeIndexes[0]);
    }
 
    // this yields global node numbers for global faces
    const int* getGlobalFaceNodeIndexes_Composite() const
    {
        return &(_compositeGlobalFaceNodeIndexes[0]);
    }
 
    int getFaceElementIndexes_Interconnected_BlockSize() const
    {
        return _interconnectedFaceElementIndexesBlockSize;
    }
 
    int getFaceNodeIndexes_Interconnected_BlockSize() const
    {
        return _interconnectedFaceNodeIndexesBlockSize;
    }
    int getFaceNodeIndexes_Composite_BlockSize() const
    {
        return _compositeFaceNodeIndexesBlockSize;
    }
 
    int getFaceElementIndexes_Interconnected_Size() const
    {
        return _numFaces * _interconnectedFaceElementIndexesBlockSize;
    }
 
    int getFaceNodeIndexes_Interconnected_Size() const
    {
        return _numFaces * _interconnectedFaceNodeIndexesBlockSize;
    }
 
    int getFaceNodeIndexes_Composite_Size() const
    {
        return _numElements_PIPEG * _compositeFaceNodeIndexesBlockSize;
    }
 
    // Some Indexing functions
    int getElementIndexes_NumInLayers(int lowerLayer, int upperLayer) const;
    int getFaceIndexes_NumInLayers(int lowerLayer, int upperLayer) const;
 
    int getElementIndexes_StartIndexForLayer(int layer) const;
    int getFaceIndexes_StartIndexForLayer(int layer) const;
 
    bool hasBoundary(const std::string& boundary) const;
 
    // iman addition
    int isBoundaryFace(int faceIndex) const
    {
        return _isBoundaryFaceArray[faceIndex];
    }
 
    // iman addition
    int isPatchBoundaryFace(int faceIndex) const
    {
        return _isPatchBoundaryFaceArray[faceIndex];
    }
    // iman addition
    const real* getCommonNormal(int faceIndex) const
    {
        return &_commonNormals[faceIndex * 3];
    }
 
    const WmUnstructuredPatch& getPatch() const
    {
        return _patch;
    }
 
    const WmUnstructuredGeometry* get_unstructured_geometry() const
    {
        return _uG;
    }
 
    const face_t& get_face(const int global_face_index) const
    {
        return *_all_faces[global_face_index];
    }
 
    int get_global_face_from_element_local_face(int global_element_index,
                                                int local_face_index) const
    {
        return _globalFaceNumbersOfGlobalElementNumbers[global_element_index *
                                                            _uG->getNumFacesPerElement() +
                                                        local_face_index];
    }
 
    // print the face data block for a given global face index
    // // globalFaceGeometry[effectiveDx, (nx,ny,nz), ((p0x, p0y, p0z), (p1x, p1y,
    // p1z), ...)] - per face
    void printFaceInformation(int global_face_index) const
    {
        std::stringstream ss;
        const int numDims = _basisSet.getNumDims();
        const int numFaceNodes = _basisSet.getNumFaceNodes();
        const real* fGeoArray =
            &(_faceGeometry.at(global_face_index * _faceGeometryBlockSize));
 
        ss << "global_face_index = " << global_face_index << " on patch "
           << getPatch().getPatchIndex() << std::endl;
        ss << "effective_dx = " << fGeoArray[0] << std::endl;
        ss << "normal = (";
        for (int i = 0; i < numDims; i++)
        {
            ss << fGeoArray[1 + i] << ", ";
        }
        ss << ")" << std::endl;
 
        ss << "nodes = (";
        for (int i = 0; i < numFaceNodes; i++)
        {
            ss << "(";
            for (int j = 0; j < numDims; j++)
            {
                ss << fGeoArray[1 + numDims + i * numDims + j] << ", ";
            }
            ss << "), ";
        }
        ss << ")" << std::endl;
 
        std::cout << ss.str();
    }
 
    // print the Element data block for a given global element index
    // globalElementGeometry[(dx,dy,dz), ((p0x, p0y, p0z), (p1x, p1y, p1z),
    // ...)] - per element (inclusive)
 
    void printElementInformation(int global_element_index)
    {
        std::stringstream ss;
 
        const int numDims = _basisSet.getNumDims();
        const int numElementNodes = _basisSet.getNumTotalNodes();
        const real* eGeomArray =
            &(_elementGeometry[0]) + global_element_index * _elementGeometryBlockSize;
 
        ss << "global_element_index = " << global_element_index << " on patch "
           << getPatch().getPatchIndex() << std::endl;
        ss << "(dx,dy,dv) = (";
        for (int i = 0; i < numDims; i++)
        {
            ss << eGeomArray[i] << ", ";
        }
        ss << ")" << std::endl;
 
        ss << "nodes = (";
        for (int i = 0; i < numElementNodes; i++)
        {
            ss << "(";
            for (int j = 0; j < numDims; j++)
            {
                ss << eGeomArray[numDims + i * numDims + j] << ", ";
            }
            ss << "), ";
        }
        ss << ")" << std::endl;
 
        std::cout << ss.str();
    }
 
    int getNumL1Elements() const
    {
        return _l1_elements.size();
    }
 
    // The following 2 functions go hand in hand to find the global node index of some
    // position within an element
 
    void getElementGlobalNodePositionsAndIndexes(const int element_index,
                                                 int* global_node_indexes,
                                                 real* global_node_positions) const;
 
    int getGlobalNodeIndexAtPosition(const real* global_node_positions,
                                     const int* global_node_indexes,
                                     const real position[3]) const;
 
    void printDGNeighborhood(std::stringstream& ss) const;
 
    void printElementwiseFaceInsideOutsideNodes(std::stringstream& ss) const;
 
    const int* getDGNeighborhood() const
    {
        return _dg_neighborhood.data();
    }
 
    const int* getLocalDGNeighborhood(const int element_index) const
    {
        return _dg_neighborhood.data() + element_index * _uG->getNumFacesPerElement();
    }
 
    int getNeighborLocalFaceIndex(const int element_index,
                                  const int local_face_index) const
    {
        return _dg_neighbor_local_face[element_index * _uG->getNumFacesPerElement() +
                                       local_face_index];
    }
 
    int getFaceOrientation(const int element_index, const int local_face_index) const
    {
        return _dg_orientations[element_index * _uG->getNumFacesPerElement() +
                                local_face_index];
    }
 
    const int* getFaceNodeGlobalInsideOutsideIndexes(const int element_index,
                                                     const int local_face_index,
                                                     const int local_node_index) const
    {
        const int numFacesPerElement = _uG->getNumFacesPerElement();
        const int numFaceNodes = _basisSet.getNumFaceNodes();
        return _inside_outside_global_node_indexes.data() +
               element_index * numFacesPerElement * numFaceNodes * 2 +
               local_face_index * numFaceNodes * 2 + local_node_index * 2;
    }
 
    void printCompositeFaceNodeIndexesArray(std::stringstream& ss) const;
 
    std::shared_ptr<WxRange> ownElementRange() const;
 
    template<typename Func> void forEachNode(Func f)
    {
        WxRange scope = *ownElementRange();
        int N_p = _basisSet.getNumTotalNodes();
        int ndims = _basisSet.getNumDims();
 
        std::vector<real> x(ndims);
 
        const real* globalElementGeometryData = getGeometryArray_Element();
        const int element_geometry_block_size = getGeometryArray_Element_BlockSize();
 
        for (int elementIndex = scope.lower(0); elementIndex <= scope.upper(0);
             elementIndex++)
        {
            for (int localNodeIndex = 0; localNodeIndex < N_p; localNodeIndex++)
            {
                const int globalNodeIndex = elementIndex * N_p + localNodeIndex;
 
                x[0] =
                    globalElementGeometryData[element_geometry_block_size * elementIndex +
                                              ndims * (localNodeIndex + 1) + 0];
                if (ndims > 1)
                {
                    x[1] = globalElementGeometryData[element_geometry_block_size *
                                                         elementIndex +
                                                     ndims * (localNodeIndex + 1) + 1];
                }
                if (ndims > 2)
                {
                    x[2] = globalElementGeometryData[element_geometry_block_size *
                                                         elementIndex +
                                                     ndims * (localNodeIndex + 1) + 2];
                }
 
                f(globalNodeIndex, x);
            }
        }
    }
 
protected:
    const WmUnstructuredPatch& _patch;
    const WmUnstructuredGeometry* _uG;
    const WmBasisArraySet& _basisSet;
 
    std::vector<std::string> _boundaryNames;
    boundaryFacesMap_t _boundaryFaceIndexes;
 
    // vectors and maps of face objects for the various types
    std::vector<std::shared_ptr<face_t>> _interior_faces;
    std::vector<std::shared_ptr<face_t>> _all_faces;
    bcFaceMap_t _boundary_faces;
    bcFaceMap_t _shared_faces;
 
    int _numFaces;
    int _numBoundaryFaces;
    int _numSharedFaces;
    int _numInteriorFaces;
 
    // Number of elements in patch
    int _numElements_PEG;
 
    // Number of elements in patch, interpatch elements group
    int _numElements_PIPEG;
 
    // Number of elements in patch, interpatch, boundary condition elements
    // group
    int _numElements_PIPBCEG;
 
    int _faceGeometryBlockSize;
    int _elementGeometryBlockSize;
    int _jacobianGeometryBlockSize;
    int _limiterGeometryBlockSize;
 
    std::vector<real> _faceGeometry;
    // std::vector<std::unique_ptr<face_t>> _face_t_array;
    std::vector<real> _elementGeometry;
    std::vector<real> _jacobianGeometry;
    std::vector<real> _limiterGeometry;
 
    int _reconstruction_stencilIndexes_blockSize;
    int _reconstruction_w_blockSize;
    int _reconstruction_lInverse_blockSize;
    int _reconstruction_b_blockSize;
 
    std::vector<int> _reconstruction_stencilIndexes;
    std::vector<real> _reconstruction_w;
    std::vector<real> _reconstruction_lInverse;
    std::vector<real> _reconstruction_b;
 
    int _interconnectedFaceNodeIndexesBlockSize;
    int _interconnectedFaceElementIndexesBlockSize;
    int _compositeFaceNodeIndexesBlockSize;
 
    std::vector<int> _numFacesPerLayer;
 
    std::vector<int> _interconnectedFaceNodeIndexes;
    std::vector<int>
        _interconnectedFaceElementIndexes; // indexes of element numbers associated with
                                           // face, first is inside face, second is
                                           // outside face
 
    // this is same as _interconnectedFaceElementIndexes except outside boundary elements
    // are not -1 but have numbers added to end of the ther elements in the patch
    std::vector<int> _interconnectedFaceElementIndexesNoMinusOne;
 
    std::vector<int> _compositeFaceNodeIndexes; // local node indexes for global face
    std::vector<int>
        _compositeGlobalFaceNodeIndexes; // global node indexes for global face (added by
                                         // iman nov 25, 2017 for use by limiter)
    std::vector<int> _isBoundaryFaceArray; // iman addition aug 29 2017 for LDG to know if
                                           // face is on boundary or not
    std::vector<int>
        _isPatchBoundaryFaceArray; // iman addition sep 19 2017 for LDG to know if face is
                                   // on a patch boundary or not
    std::vector<real>
        _commonNormals; // iman addition sep 20 2017 for LDG common face normals
 
    std::vector<real>
        _globalFaceNumbersOfGlobalElementNumbers; // this gives the global face numbers of
                                                  // the element faces (should be in order
                                                  // of local faces from 0 to whatever).
                                                  // Think of it as neighborhood but
                                                  // instead of neighboring element
                                                  // numbers it gives the global face
                                                  // number between the neighbor and
                                                  // current element
 
    std::vector<int> _l1_elements; // element neighbours patch
 
    std::vector<real> _numerical_flux_multipliers;
 
    std::vector<int> _dg_neighborhood;
 
    std::vector<int> _dg_neighbor_local_face;
 
    std::vector<int> _dg_orientations;
 
    std::vector<int> _inside_outside_global_node_indexes;
 
private:
    // functions for constructor
    void constructFacesAndLayer1Vector();
 
    void setInteriorFaceIndexes();
 
    void setNumFacesPerLayer();
 
    void constructAllFacesVectorAndCountNumFaces();
 
    void setupLimiters(const bool useLimiter);
 
    void setupFaceGeometry();
 
    void setupElementGeometry(const bool useLimiter);
 
    void setupDGNeighborhoodLocalNeighborFacesAndOrientations();
 
    int boundaryFaceIndex(const int global_face_index) const;
 
    int boundaryElementIndex(const int global_face_index) const;
 
    void setupElementwiseFaceInsideOutsideNodes();
 
    //    WmUDGGeometry(){}
};
#endif // WMUDGGEOMETRY_H