Line data    Source code

       1              : namespace ippl {
       2              :     template <typename T>
       3              :     KOKKOS_FUNCTION typename HexahedralElement<T>::vertex_points_t
       4       256032 :     HexahedralElement<T>::getLocalVertices() const {
       5              :         // For the ordering of local vertices, see section 3.3.1:
       6              :         // https://amas.web.psi.ch/people/aadelmann/ETH-Accel-Lecture-1/projectscompleted/phys/bachelor_thesis_buehlluk.pdf
       7      2304288 :         HexahedralElement::vertex_points_t vertices = {
       8              :             {0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {1.0, 1.0, 0.0}, {0.0, 1.0, 0.0},
       9              :             {0.0, 0.0, 1.0}, {1.0, 0.0, 1.0}, {1.0, 1.0, 1.0}, {0.0, 1.0, 1.0}};
      10              : 
      11       512032 :         return vertices;
      12       512064 :     }
      13              : 
      14              :     template <typename T>
      15              :     KOKKOS_FUNCTION typename HexahedralElement<T>::point_t
      16           16 :     HexahedralElement<T>::getTransformationJacobian(
      17              :         const HexahedralElement<T>::vertex_points_t& global_vertices) const {
      18           16 :         HexahedralElement::point_t jacobian;
      19              : 
      20           16 :         jacobian[0] = (global_vertices[1][0] - global_vertices[0][0]);
      21           16 :         jacobian[1] = (global_vertices[2][1] - global_vertices[0][1]);
      22           16 :         jacobian[2] = (global_vertices[4][2] - global_vertices[0][2]);
      23              : 
      24           16 :         return jacobian;
      25              :     }
      26              : 
      27              :     template <typename T>
      28              :     KOKKOS_FUNCTION typename HexahedralElement<T>::point_t
      29           48 :     HexahedralElement<T>::getInverseTransformationJacobian(
      30              :         const HexahedralElement<T>::vertex_points_t& global_vertices) const {
      31           48 :         HexahedralElement::point_t inv_jacobian;
      32              : 
      33           48 :         inv_jacobian[0] = 1.0 / (global_vertices[1][0] - global_vertices[0][0]);
      34           48 :         inv_jacobian[1] = 1.0 / (global_vertices[2][1] - global_vertices[0][1]);
      35           48 :         inv_jacobian[2] = 1.0 / (global_vertices[4][2] - global_vertices[0][2]);
      36              : 
      37           48 :         return inv_jacobian;
      38              :     }
      39              : 
      40              :     template <typename T>
      41           48 :     KOKKOS_FUNCTION typename HexahedralElement<T>::point_t HexahedralElement<T>::globalToLocal(
      42              :         const HexahedralElement<T>::vertex_points_t& global_vertices,
      43              :         const HexahedralElement<T>::point_t& global_point) const {
      44              :         // This is actually not a matrix, but an IPPL vector that represents a diagonal matrix
      45           48 :         const HexahedralElement<T>::point_t glob2loc_matrix =
      46              :             getInverseTransformationJacobian(global_vertices);
      47              : 
      48           48 :         HexahedralElement<T>::point_t local_point =
      49           96 :             glob2loc_matrix * (global_point - global_vertices[0]);
      50              : 
      51           48 :         return local_point;
      52           48 :     }
      53              : 
      54              :     template <typename T>
      55           16 :     KOKKOS_FUNCTION typename HexahedralElement<T>::point_t HexahedralElement<T>::localToGlobal(
      56              :         const HexahedralElement<T>::vertex_points_t& global_vertices,
      57              :         const HexahedralElement<T>::point_t& local_point) const {
      58              :         // This is actually not a matrix but an IPPL vector that represents a diagonal matrix
      59           16 :         const HexahedralElement<T>::point_t loc2glob_matrix =
      60              :             getTransformationJacobian(global_vertices);
      61              : 
      62           16 :         HexahedralElement<T>::point_t global_point =
      63           32 :             (loc2glob_matrix * local_point) + global_vertices[0];
      64              : 
      65           16 :         return global_point;
      66           16 :     }
      67              : 
      68              :     template <typename T>
      69            0 :     KOKKOS_FUNCTION T HexahedralElement<T>::getDeterminantOfTransformationJacobian(
      70              :         const HexahedralElement<T>::vertex_points_t& global_vertices) const {
      71            0 :         T determinant = 1.0;
      72              : 
      73              :         // Since the jacobian is a diagonal matrix in our case the determinant is the product of the
      74              :         // diagonal elements
      75            0 :         for (const T& jacobian_val : getTransformationJacobian(global_vertices)) {
      76            0 :             determinant *= jacobian_val;
      77              :         }
      78              : 
      79            0 :         return determinant;
      80              :     }
      81              : 
      82              :     template <typename T>
      83              :     KOKKOS_FUNCTION typename HexahedralElement<T>::point_t
      84            0 :     HexahedralElement<T>::getInverseTransposeTransformationJacobian(
      85              :         const HexahedralElement<T>::vertex_points_t& global_vertices) const {
      86              :         // Simply return the inverse transformation jacobian since it is a diagonal matrix
      87            0 :         return getInverseTransformationJacobian(global_vertices);
      88              :     }
      89              : 
      90              :     template <typename T>
      91       256000 :     KOKKOS_FUNCTION bool HexahedralElement<T>::isPointInRefElement(
      92              :         const Vector<T, 3>& point) const {
      93              :         // check if the local coordinates are inside the reference element
      94              : 
      95      1024000 :         for (size_t d = 0; d < 3; d++) {
      96       768000 :             if (point[d] > 1.0 || point[d] < 0.0) {
      97              :                 // The global coordinates are outside of the support.
      98            0 :                 return false;
      99              :             }
     100              :         }
     101              : 
     102       256000 :         return true;
     103              :     }
     104              : 
     105              : }  // namespace ippl