Branch data     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