ConformingVirtualElementOrder1.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_VIRTUALELEMENT_CONFORMINGVIRTUALELEMENTORDER1_HPP_
#define GEOS_VIRTUALELEMENT_CONFORMINGVIRTUALELEMENTORDER1_HPP_
 
#include "FiniteElementBase.hpp"
#include "codingUtilities/traits.hpp"
 
namespace geos
{
namespace finiteElement
{
template< localIndex MAXCELLNODES, localIndex MAXFACENODES >
class ConformingVirtualElementOrder1_impl : public FiniteElementBase_impl< MAXCELLNODES, MAXFACENODES, 1 >
{
public:
  using InputNodeCoords = arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD >;
  template< typename SUBREGION_TYPE >
  using InputCellToNodeMap = traits::ViewTypeConst< typename SUBREGION_TYPE::NodeMapType >;
  template< typename SUBREGION_TYPE >
  using InputCellToFaceMap = traits::ViewTypeConst< typename SUBREGION_TYPE::FaceMapType >;
  using InputFaceToNodeMap = ArrayOfArraysView< localIndex const >;
  using InputFaceToEdgeMap = ArrayOfArraysView< localIndex const >;
  using InputEdgeToNodeMap = arrayView2d< localIndex const >;
 
  static constexpr localIndex maxSupportPoints = MAXCELLNODES;
  static constexpr localIndex numNodes = MAXCELLNODES;
  static constexpr localIndex numQuadraturePoints = 1;
 
  struct StackVariables
  {
 
    localIndex numSupportPoints;
    real64 quadratureWeight;
    real64 basisFunctionsIntegralMean[MAXCELLNODES];
    real64 stabilizationMatrix[MAXCELLNODES][MAXCELLNODES];
    real64 basisDerivativesIntegralMean[MAXCELLNODES][3];
  };
 
  template< typename SUBREGION_TYPE >
  struct MeshData
  {
    InputNodeCoords nodesCoords;
    InputCellToNodeMap< SUBREGION_TYPE > cellToNodeMap;
    InputCellToFaceMap< SUBREGION_TYPE > cellToFaceMap;
    InputFaceToNodeMap faceToNodeMap;
    InputFaceToEdgeMap faceToEdgeMap;
    InputEdgeToNodeMap edgeToNodeMap;
    arrayView2d< real64 const > faceCenters;
    arrayView2d< real64 const > faceNormals;
    arrayView1d< real64 const > faceAreas;
    arrayView2d< real64 const > cellCenters;
    arrayView1d< real64 const > cellVolumes;
  };
 
  GEOS_HOST_DEVICE
  inline static
  localIndex getNumQuadraturePoints()
  {
    return numQuadraturePoints;
  }
 
  GEOS_HOST_DEVICE
  inline static
  localIndex getMaxSupportPoints()
  {
    return maxSupportPoints;
  }
 
  GEOS_HOST_DEVICE
  inline
  static localIndex getNumSupportPoints( StackVariables const & stack )
  {
    return stack.numSupportPoints;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static
  real64 calcJacobian( localIndex const q,
                       real64 const (&X)[numNodes][3],
                       StackVariables const & stack,
                       real64 (& J)[3][3] )
  {
    GEOS_UNUSED_VAR( q );
    GEOS_UNUSED_VAR( X );
    for( localIndex i = 0; i < 3; ++i )
    {
      for( localIndex j = 0; j < 3; ++j )
      {
        J[i][j] = ( i == j ) ? 1.0 : 0.0;
      }
    }
    return stack.quadratureWeight; // * 1.0;
  }
 
  GEOS_HOST_DEVICE
  inline
  static real64 calcGradN( localIndex const q,
                           real64 const (&X)[maxSupportPoints][3],
                           StackVariables const & stack,
                           real64 ( & gradN )[maxSupportPoints][3] )
  {
    for( localIndex i = 0; i < stack.numSupportPoints; ++i )
    {
      gradN[i][0] = stack.basisDerivativesIntegralMean[i][0];
      gradN[i][1] = stack.basisDerivativesIntegralMean[i][1];
      gradN[i][2] = stack.basisDerivativesIntegralMean[i][2];
    }
    return transformedQuadratureWeight( q, X, stack );
  }
 
  GEOS_HOST_DEVICE
  inline
  static void calcN( localIndex const q,
                     StackVariables const & stack,
                     real64 ( & N )[maxSupportPoints] )
  {
    GEOS_UNUSED_VAR( q );
    for( localIndex i = 0; i < stack.numSupportPoints; ++i )
    {
      N[i] = stack.basisFunctionsIntegralMean[i];
    }
  }
 
  template< localIndex NUMDOFSPERTRIALSUPPORTPOINT, localIndex MAXSUPPORTPOINTS, bool UPPER >
  GEOS_HOST_DEVICE
  inline
  static void addGradGradStabilization( StackVariables const & stack,
                                        real64 ( & matrix )
                                        [MAXSUPPORTPOINTS * NUMDOFSPERTRIALSUPPORTPOINT]
                                        [MAXSUPPORTPOINTS * NUMDOFSPERTRIALSUPPORTPOINT],
                                        real64 const & scaleFactor )
  {
    for( localIndex i = 0; i < stack.numSupportPoints; ++i )
    {
      localIndex startCol = (UPPER) ? i : 0;
      for( localIndex j = startCol; j < stack.numSupportPoints; ++j )
      {
        real64 const contribution = scaleFactor * stack.stabilizationMatrix[i][j];
        for( localIndex d = 0; d < NUMDOFSPERTRIALSUPPORTPOINT; ++d )
        {
          matrix[i*NUMDOFSPERTRIALSUPPORTPOINT + d][j*NUMDOFSPERTRIALSUPPORTPOINT + d] += contribution;
        }
      }
    }
  }
 
  template< localIndex NUMDOFSPERTRIALSUPPORTPOINT, localIndex MAXSUPPORTPOINTS >
  GEOS_HOST_DEVICE
  inline
  static void addEvaluatedGradGradStabilization( StackVariables const & stack,
                                                 real64 const ( &dofs )[MAXSUPPORTPOINTS][NUMDOFSPERTRIALSUPPORTPOINT],
                                                 real64 ( & targetVector )[MAXSUPPORTPOINTS][NUMDOFSPERTRIALSUPPORTPOINT],
                                                 real64 const scaleFactor )
  {
    for( localIndex d = 0; d < NUMDOFSPERTRIALSUPPORTPOINT; ++d )
    {
      for( localIndex i = 0; i < stack.numSupportPoints; ++i )
      {
        for( localIndex j = 0; j < stack.numSupportPoints; ++j )
        {
          targetVector[i][d] += scaleFactor * stack.stabilizationMatrix[i][j] * dofs[j][d];
        }
      }
    }
  }
 
  GEOS_HOST_DEVICE
  inline
  static real64 transformedQuadratureWeight( localIndex const q,
                                             real64 const ( &X )[maxSupportPoints][3],
                                             StackVariables const & stack )
  {
    GEOS_UNUSED_VAR( q, X );
    return stack.quadratureWeight;
  }
 
  template< typename SUBREGION_TYPE >
  inline
  static void fillMeshData( NodeManager const & nodeManager,
                            EdgeManager const & edgeManager,
                            FaceManager const & faceManager,
                            SUBREGION_TYPE const & cellSubRegion,
                            MeshData< SUBREGION_TYPE > & meshData
                            )
  {
    meshData.nodesCoords = nodeManager.referencePosition();
    meshData.cellToNodeMap = cellSubRegion.nodeList().toViewConst();
    meshData.cellToFaceMap = cellSubRegion.faceList().toViewConst();
    meshData.faceToNodeMap = faceManager.nodeList().toViewConst();
    meshData.faceToEdgeMap = faceManager.edgeList().toViewConst();
    meshData.edgeToNodeMap = edgeManager.nodeList().toViewConst();
    meshData.faceCenters = faceManager.faceCenter();
    meshData.faceNormals = faceManager.faceNormal();
    meshData.faceAreas = faceManager.faceArea();
    meshData.cellCenters = cellSubRegion.getElementCenter();
    meshData.cellVolumes = cellSubRegion.getElementVolume();
  }
 
  template< typename SUBREGION_TYPE >
  GEOS_HOST_DEVICE
  inline
  static void setupStack( localIndex const & cellIndex,
                          MeshData< SUBREGION_TYPE > const & meshData,
                          StackVariables & stack )
  {
    real64 const cellCenter[3] { meshData.cellCenters( cellIndex, 0 ),
                                 meshData.cellCenters( cellIndex, 1 ),
                                 meshData.cellCenters( cellIndex, 2 ) };
    real64 const cellVolume = meshData.cellVolumes( cellIndex );
 
    computeProjectors< SUBREGION_TYPE >( cellIndex,
                                         meshData.nodesCoords,
                                         meshData.cellToNodeMap,
                                         meshData.cellToFaceMap,
                                         meshData.faceToNodeMap,
                                         meshData.faceToEdgeMap,
                                         meshData.edgeToNodeMap,
                                         meshData.faceCenters,
                                         meshData.faceNormals,
                                         meshData.faceAreas,
                                         cellCenter,
                                         cellVolume,
                                         stack.numSupportPoints,
                                         stack.quadratureWeight,
                                         stack.basisFunctionsIntegralMean,
                                         stack.stabilizationMatrix,
                                         stack.basisDerivativesIntegralMean );
  }
 
  GEOS_HOST_DEVICE
  inline static
  localIndex getNumSupportPoints()
  {
    GEOS_ERROR( "VEM functions have to be called with the StackVariables syntax" );
    return 0;
  }
 
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static void getSamplingPointCoordInParentSpace( int const & linearIndex,
                                                  real64 (& samplingPointCoord)[3] )
  {
    GEOS_UNUSED_VAR( linearIndex, samplingPointCoord );
    GEOS_ERROR( "Element type not supported." );
  }
 
  GEOS_HOST_DEVICE
  inline
  static void calcN( localIndex const q,
                     real64 ( & N )[maxSupportPoints] )
  {
    GEOS_ERROR( "VEM functions have to be called with the StackVariables syntax" );
    GEOS_UNUSED_VAR( q );
    for( localIndex i = 0; i < maxSupportPoints; ++i )
    {
      N[i] = 0.0;
    }
  }
 
  GEOS_HOST_DEVICE
  inline
  static void calcN( real64 const (&coords)[3],
                     real64 (& N)[maxSupportPoints] )
  {
    GEOS_ERROR( "VEM functions have to be called with the StackVariables syntax" );
    GEOS_UNUSED_VAR( coords );
    for( localIndex i = 0; i < maxSupportPoints; ++i )
    {
      N[i] = 0.0;
    }
  }
 
  GEOS_HOST_DEVICE
  inline
  static real64 invJacobianTransformation( int const q,
                                           real64 const ( &X )[numNodes][3],
                                           real64 ( & J )[3][3] )
  {
    GEOS_ERROR( "No reference element map is defined for VEM classes" );
    GEOS_UNUSED_VAR( q, X );
    for( localIndex i = 0; i < 3; ++i )
    {
      for( localIndex j = 0; j < 3; ++j )
      {
        J[i][j] = 0.0;
      }
    }
    return 0.0;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static real64 calcJacobian( localIndex const q,
                              real64 const (&X)[numNodes][3],
                              real64 (& J)[3][3] )
  {
    GEOS_ERROR( "VEM functions have to be called with the StackVariables syntax" );
    GEOS_UNUSED_VAR( q, X );
    for( localIndex i = 0; i < 3; ++i )
    {
      for( localIndex j = 0; j < 3; ++j )
      {
        J[i][j] = 0.0;
      }
    }
    return 0.0;
  }
 
  GEOS_HOST_DEVICE
  inline
  static real64 calcGradN( localIndex const q,
                           real64 const ( &X )[maxSupportPoints][3],
                           real64 ( & gradN )[maxSupportPoints][3] )
  {
    GEOS_ERROR( "VEM functions have to be called with the StackVariables syntax" );
    GEOS_UNUSED_VAR( q, X );
    for( localIndex i = 0; i < maxSupportPoints; ++i )
    {
      for( localIndex j = 0; j < 3; ++j )
      {
        gradN[i][j] = 0.0;
      }
    }
    return 0.0;
  }
 
  GEOS_HOST_DEVICE
  real64 transformedQuadratureWeight( localIndex const q,
                                      real64 const ( &X )[maxSupportPoints][3] ) const
  {
    GEOS_ERROR( "VEM functions have to be called with the StackVariables syntax" );
    GEOS_UNUSED_VAR( q, X );
    return 0.0;
  }
 
private:
 
  GEOS_HOST_DEVICE
  static void
    computeFaceIntegrals( InputNodeCoords const & nodesCoords,
                          localIndex const (&faceToNodes)[MAXFACENODES],
                          localIndex const (&faceToEdges)[MAXFACENODES],
                          localIndex const & numFaceVertices,
                          real64 const & faceArea,
                          real64 const (&faceCenter)[3],
                          real64 const (&faceNormal)[3],
                          InputEdgeToNodeMap const & edgeToNodes,
                          real64 const & invCellDiameter,
                          real64 const (&cellCenter)[3],
                          real64 ( &basisIntegrals )[MAXFACENODES],
                          real64 ( &threeDMonomialIntegrals )[3] );
 
  template< typename SUBREGION_TYPE >
  GEOS_HOST_DEVICE
  inline
  static void
    computeProjectors( localIndex const & cellIndex,
                       InputNodeCoords const & nodesCoords,
                       InputCellToNodeMap< SUBREGION_TYPE > const & cellToNodeMap,
                       InputCellToFaceMap< SUBREGION_TYPE > const & elementToFaceMap,
                       InputFaceToNodeMap const & faceToNodeMap,
                       InputFaceToEdgeMap const & faceToEdgeMap,
                       InputEdgeToNodeMap const & edgeToNodeMap,
                       arrayView2d< real64 const > const faceCenters,
                       arrayView2d< real64 const > const faceNormals,
                       arrayView1d< real64 const > const faceAreas,
                       real64 const (&cellCenter)[3],
                       real64 const & cellVolume,
                       localIndex & numSupportPoints,
                       real64 & quadratureWeight,
                       real64 ( &basisFunctionsIntegralMean )[MAXCELLNODES],
                       real64 ( &stabilizationMatrix )[MAXCELLNODES][MAXCELLNODES],
                       real64 ( &basisDerivativesIntegralMean )[MAXCELLNODES][3]
                       );
 
  template< localIndex DIMENSION, typename POINT_COORDS_TYPE >
  GEOS_HOST_DEVICE
  inline
  static real64 computeDiameter( POINT_COORDS_TYPE points,
                                 localIndex const & numPoints )
  {
    real64 diameter = 0;
    for( localIndex numPoint = 0; numPoint < numPoints; ++numPoint )
    {
      for( localIndex numOthPoint = 0; numOthPoint < numPoint; ++numOthPoint )
      {
        real64 candidateDiameter = 0.0;
        for( localIndex i = 0; i < DIMENSION; ++i )
        {
          real64 coordDiff = points[numPoint][i] - points[numOthPoint][i];
          candidateDiameter += coordDiff * coordDiff;
        }
        if( diameter < candidateDiameter )
        {
          diameter = candidateDiameter;
        }
      }
    }
    return LvArray::math::sqrt< real64 >( diameter );
  }
 
  template< localIndex DIMENSION, typename POINT_COORDS_TYPE, typename POINT_SELECTION_TYPE >
  GEOS_HOST_DEVICE
  inline
  static real64 computeDiameter( POINT_COORDS_TYPE points,
                                 POINT_SELECTION_TYPE selectedPoints,
                                 localIndex const & numSelectedPoints )
  {
    real64 diameter = 0;
    for( localIndex numPoint = 0; numPoint < numSelectedPoints; ++numPoint )
    {
      for( localIndex numOthPoint = 0; numOthPoint < numPoint; ++numOthPoint )
      {
        real64 candidateDiameter = 0.0;
        for( localIndex i = 0; i < DIMENSION; ++i )
        {
          real64 coordDiff = points[selectedPoints[numPoint]][i] -
                             points[selectedPoints[numOthPoint]][i];
          candidateDiameter += coordDiff * coordDiff;
        }
        if( diameter < candidateDiameter )
        {
          diameter = candidateDiameter;
        }
      }
    }
    return LvArray::math::sqrt< real64 >( diameter );
  }
};
 
 
 
template< localIndex MAXCELLNODES, localIndex MAXFACENODES >
class ConformingVirtualElementOrder1 final : public ConformingVirtualElementOrder1_impl< MAXCELLNODES, MAXFACENODES >,
  public FiniteElementBase
{
public:
  using ImplType = ConformingVirtualElementOrder1_impl< MAXCELLNODES, MAXFACENODES >;
 
  constexpr static localIndex numNodes = ConformingVirtualElementOrder1_impl< MAXCELLNODES, MAXFACENODES >::numNodes;
 
  constexpr static localIndex numFaces = ConformingVirtualElementOrder1_impl< MAXCELLNODES, MAXFACENODES >::numFaces;
 
  constexpr static localIndex maxSupportPoints = numNodes;
 
  constexpr static localIndex numQuadraturePoints = ConformingVirtualElementOrder1_impl< MAXCELLNODES, MAXFACENODES >::numQuadraturePoints;
 
  ConformingVirtualElementOrder1_impl< MAXCELLNODES, MAXFACENODES > * getImpl()
  {
    return static_cast< ConformingVirtualElementOrder1_impl< MAXCELLNODES, MAXFACENODES > * >(this);
  }
 
  const ConformingVirtualElementOrder1_impl< MAXCELLNODES, MAXFACENODES > * getImpl() const
  {
    return static_cast< const ConformingVirtualElementOrder1_impl< MAXCELLNODES, MAXFACENODES > * >(this);
  }
 
  ConformingVirtualElementOrder1():
    FiniteElementBase( numNodes,
                       maxSupportPoints,
                       numQuadraturePoints )
  {}
 
  virtual ~ConformingVirtualElementOrder1() override final = default;
 
 
};
 
using H1_Tetrahedron_VEM_Gauss1 = ConformingVirtualElementOrder1< 4, 3 >;
#if !defined( GEOS_USE_HIP )
using H1_Hexahedron_VEM_Gauss1 = ConformingVirtualElementOrder1< 8, 4 >;
#endif
using H1_Pyramid_VEM_Gauss1 = ConformingVirtualElementOrder1< 5, 4 >;
#if !defined( GEOS_USE_HIP )
using H1_Wedge_VEM_Gauss1 = ConformingVirtualElementOrder1< 6, 4 >;
#endif
using H1_Prism5_VEM_Gauss1 = ConformingVirtualElementOrder1< 10, 5 >;
using H1_Prism6_VEM_Gauss1 = ConformingVirtualElementOrder1< 12, 6 >;
using H1_Prism7_VEM_Gauss1 = ConformingVirtualElementOrder1< 14, 7 >;
using H1_Prism8_VEM_Gauss1 = ConformingVirtualElementOrder1< 16, 8 >;
using H1_Prism9_VEM_Gauss1 = ConformingVirtualElementOrder1< 18, 9 >;
using H1_Prism10_VEM_Gauss1 = ConformingVirtualElementOrder1< 20, 10 >;
#if !defined( GEOS_USE_HIP )
using H1_Prism11_VEM_Gauss1 = ConformingVirtualElementOrder1< 22, 11 >;
#endif
 
 
using H1_Tetrahedron_VEM_Gauss1_impl = ConformingVirtualElementOrder1_impl< 4, 3 >;
#if !defined( GEOS_USE_HIP )
using H1_Hexahedron_VEM_Gauss1_impl = ConformingVirtualElementOrder1_impl< 8, 4 >;
#endif
using H1_Pyramid_VEM_Gauss1_impl = ConformingVirtualElementOrder1_impl< 5, 4 >;
#if !defined( GEOS_USE_HIP )
using H1_Wedge_VEM_Gauss1_impl = ConformingVirtualElementOrder1_impl< 6, 4 >;
#endif
using H1_Prism5_VEM_Gauss1_impl = ConformingVirtualElementOrder1_impl< 10, 5 >;
using H1_Prism6_VEM_Gauss1_impl = ConformingVirtualElementOrder1_impl< 12, 6 >;
using H1_Prism7_VEM_Gauss1_impl = ConformingVirtualElementOrder1_impl< 14, 7 >;
using H1_Prism8_VEM_Gauss1_impl = ConformingVirtualElementOrder1_impl< 16, 8 >;
using H1_Prism9_VEM_Gauss1_impl = ConformingVirtualElementOrder1_impl< 18, 9 >;
using H1_Prism10_VEM_Gauss1_impl = ConformingVirtualElementOrder1_impl< 20, 10 >;
#if !defined( GEOS_USE_HIP )
using H1_Prism11_VEM_Gauss1_impl = ConformingVirtualElementOrder1_impl< 22, 11 >;
#endif
 
} // namespace finiteElement
} // namespace geos
 
#include "ConformingVirtualElementOrder1_impl.hpp"
 
#endif // GEOS_VIRTUALELEMENT_CONFORMINGVIRTUALELEMENTORDER1_HPP_