Line data    Source code

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