LCOV - code coverage report
Current view: top level - src/FEM - FiniteElementSpace.hpp (source / functions) Coverage Total Hit
Test: final_report.info Lines: 92.3 % 104 96
Test Date: 2025-07-17 08:40:11 Functions: 66.1 % 109 72

            Line data    Source code
       1              : 
       2              : namespace ippl {
       3              :     template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
       4              :               typename QuadratureType, typename FieldLHS, typename FieldRHS>
       5          198 :     FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS,
       6              :                        FieldRHS>::FiniteElementSpace(UniformCartesian<T, Dim>& mesh,
       7              :                                                      ElementType& ref_element,
       8              :                                                      const QuadratureType& quadrature)
       9          198 :         : mesh_m(mesh)
      10              :         , ref_element_m(ref_element)
      11          198 :         , quadrature_m(quadrature) {
      12              :         assert(mesh.Dimension == Dim && "Mesh dimension does not match the dimension of the space");
      13              : 
      14          198 :         nr_m     = mesh_m.getGridsize();
      15          198 :         hr_m     = mesh_m.getMeshSpacing();
      16          198 :         origin_m = mesh_m.getOrigin();
      17              : 
      18              :         /*for (size_t d = 0; d < Dim; ++d) {
      19              :             assert(nr_m[d] > 1 && "Mesh has no cells in at least one dimension");
      20              :         }*/
      21          198 :     }
      22              : 
      23              :     template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
      24              :               typename QuadratureType, typename FieldLHS, typename FieldRHS>
      25            0 :     void FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS,
      26              :                        FieldRHS>::setMesh(UniformCartesian<T, Dim>& mesh)
      27              :     {
      28              :         assert(mesh.Dimension == Dim && "Mesh dimension does not match the dimension of the space");
      29              : 
      30            0 :         mesh_m = mesh;
      31              : 
      32            0 :         nr_m     = mesh_m.getGridsize();
      33            0 :         hr_m     = mesh_m.getMeshSpacing();
      34            0 :         origin_m = mesh_m.getOrigin();
      35              : 
      36            0 :         for (size_t d = 0; d < Dim; ++d) {
      37            0 :             assert(nr_m[d] > 1 && "Mesh has no cells in at least one dimension");
      38              :         }
      39            0 :     }
      40              : 
      41              :     template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
      42              :               typename QuadratureType, typename FieldLHS, typename FieldRHS>
      43           12 :     KOKKOS_FUNCTION size_t FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType,
      44              :                                               FieldLHS, FieldRHS>::numElements() const {
      45           12 :         Vector<size_t, Dim> cells_per_dim = nr_m - 1u;
      46              : 
      47           12 :         size_t num_elements = 1;
      48           42 :         for (size_t d = 0; d < Dim; ++d) {
      49           30 :             num_elements *= cells_per_dim[d];
      50              :         }
      51              : 
      52           12 :         return num_elements;
      53           12 :     }
      54              : 
      55              :     template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
      56              :               typename QuadratureType, typename FieldLHS, typename FieldRHS>
      57              :     KOKKOS_FUNCTION size_t
      58           10 :     FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS,
      59              :                        FieldRHS>::numElementsInDim(const size_t& dim) const {
      60           10 :         return nr_m[dim] - 1u;
      61              :     }
      62              : 
      63              :     template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
      64              :               typename QuadratureType, typename FieldLHS, typename FieldRHS>
      65              :     KOKKOS_FUNCTION typename FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType,
      66              :                                                 FieldLHS, FieldRHS>::indices_t
      67          444 :     FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS,
      68              :                        FieldRHS>::getMeshVertexNDIndex(const size_t& vertex_index) const {
      69              :         // Copy the vertex index to the index variable we can alter during the computation.
      70          444 :         size_t index = vertex_index;
      71              : 
      72              :         // Create a vector to store the vertex indices in each dimension for the corresponding
      73              :         // vertex.
      74          444 :         indices_t vertex_indices;
      75              : 
      76              :         // This is the number of vertices in each dimension.
      77          444 :         Vector<size_t, Dim> vertices_per_dim = nr_m;
      78              : 
      79              :         // The number_of_lower_dim_vertices is the product of the number of vertices per
      80              :         // dimension, it will get divided by the current dimensions number to get the index in
      81              :         // that dimension
      82          444 :         size_t remaining_number_of_vertices = 1;
      83         1176 :         for (const size_t num_vertices : vertices_per_dim) {
      84          732 :             remaining_number_of_vertices *= num_vertices;
      85              :         }
      86              : 
      87         1176 :         for (int d = Dim - 1; d >= 0; --d) {
      88          732 :             remaining_number_of_vertices /= vertices_per_dim[d];
      89          732 :             vertex_indices[d] = index / remaining_number_of_vertices;
      90          732 :             index -= vertex_indices[d] * remaining_number_of_vertices;
      91              :         }
      92              : 
      93          444 :         return vertex_indices;
      94          444 :     };
      95              : 
      96              :     template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
      97              :               typename QuadratureType, typename FieldLHS, typename FieldRHS>
      98              :     KOKKOS_FUNCTION size_t
      99          160 :     FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS, FieldRHS>::
     100              :         getMeshVertexIndex(
     101              :             const FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS,
     102              :                                      FieldRHS>::indices_t& vertexNDIndex) const {
     103              :         // Compute the vector to multiply the ndindex with
     104          160 :         ippl::Vector<size_t, Dim> vec(1);
     105          448 :         for (size_t d = 1; d < dim; ++d) {
     106          704 :             for (size_t d2 = d; d2 < Dim; ++d2) {
     107          416 :                 vec[d2] *= nr_m[d - 1];
     108              :             }
     109              :         }
     110              : 
     111              :         // return the dot product between the vertex ndindex and vec.
     112          160 :         return vertexNDIndex.dot(vec);
     113          160 :     }
     114              : 
     115              :     // implementation of function to retrieve the index of an element in each dimension
     116              :     template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
     117              :               typename QuadratureType, typename FieldLHS, typename FieldRHS>
     118              :     KOKKOS_FUNCTION typename FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType,
     119              :                                                 FieldLHS, FieldRHS>::indices_t
     120          812 :     FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS,
     121              :                        FieldRHS>::getElementNDIndex(const size_t& element_index) const {
     122              :         // Copy the element index to the index variable we can alter during the computation.
     123          812 :         size_t index = element_index;
     124              : 
     125              :         // Create a vector to store the element indices in each dimension for the corresponding
     126              :         // element.
     127          812 :         indices_t element_nd_index;
     128              : 
     129              :         // This is the number of cells in each dimension. It is one less than the number of
     130              :         // vertices in each dimension, which is in nr_m (mesh.getGridsize()).
     131          812 :         Vector<size_t, Dim> cells_per_dim = nr_m - 1;
     132              : 
     133              :         // The number_of_lower_dim_cells is the product of all the number of cells per
     134              :         // dimension, it will get divided by the current dimension's size to get the index in
     135              :         // that dimension
     136          812 :         size_t remaining_number_of_cells = 1;
     137         2982 :         for (const size_t num_cells : cells_per_dim) {
     138         2170 :             remaining_number_of_cells *= num_cells;
     139              :         }
     140              : 
     141         2982 :         for (int d = Dim - 1; d >= 0; --d) {
     142         2170 :             remaining_number_of_cells /= cells_per_dim[d];
     143         2170 :             element_nd_index[d] = (index / remaining_number_of_cells);
     144         2170 :             index -= (element_nd_index[d]) * remaining_number_of_cells;
     145              :         }
     146              : 
     147          812 :         return element_nd_index;
     148          812 :     }
     149              : 
     150              :     // implementation of function to retrieve the global index of an element given the ndindex
     151              :     template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
     152              :               typename QuadratureType, typename FieldLHS, typename FieldRHS>
     153              :     KOKKOS_FUNCTION size_t
     154         7794 :     FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS, FieldRHS>::
     155              :         getElementIndex(
     156              :             const FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS,
     157              :                                      FieldRHS>::indices_t& ndindex) const {
     158         7794 :         size_t element_index = 0;
     159              : 
     160              :         // This is the number of cells in each dimension. It is one less than the number of
     161              :         // vertices in each dimension, which is returned by Mesh::getGridsize().
     162         7794 :         Vector<size_t, Dim> cells_per_dim = nr_m - 1;
     163              : 
     164         7794 :         size_t remaining_number_of_cells = 1;
     165              : 
     166        29340 :         for (unsigned int d = 0; d < Dim; ++d) {
     167        21546 :             element_index += ndindex[d] * remaining_number_of_cells;
     168        21546 :             remaining_number_of_cells *= cells_per_dim[d];
     169              :         }
     170              : 
     171         7794 :         return element_index;
     172         7794 :     }
     173              : 
     174              :     template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
     175              :               typename QuadratureType, typename FieldLHS, typename FieldRHS>
     176              :     KOKKOS_FUNCTION typename FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType,
     177              :                                                 FieldLHS, FieldRHS>::vertex_indices_t
     178           24 :     FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS, FieldRHS>::
     179              :         getElementMeshVertexIndices(
     180              :             const FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS,
     181              :                                      FieldRHS>::indices_t& element_nd_index) const {
     182           24 :         const Vector<size_t, Dim> num_vertices = nr_m;
     183              : 
     184           24 :         size_t smallest_vertex_index = 0;
     185           88 :         for (int d = Dim - 1; d >= 0; --d) {
     186           64 :             size_t temp_index = element_nd_index[d];
     187          120 :             for (int i = d; i >= 1; --i) {
     188           56 :                 temp_index *= num_vertices[i];
     189              :             }
     190           64 :             smallest_vertex_index += temp_index;
     191              :         }
     192              : 
     193              :         // Vector to store the vertex indices for the element
     194           24 :         vertex_indices_t vertex_indices;
     195           24 :         vertex_indices[0] = smallest_vertex_index;
     196           24 :         vertex_indices[1] = vertex_indices[0] + 1;
     197              : 
     198              :         /*
     199              :         The following for loop computes the following computations:
     200              : 
     201              :         2D:
     202              :             vertex_indices[2] = vertex_indices[0] + num_vertices[0];
     203              :             vertex_indices[3] = vertex_indices[1] + num_vertices[0];
     204              :         3D:
     205              :             vertex_indices[4] = vertex_indices[0] + (num_vertices[0] * num_vertices[1]);
     206              :             vertex_indices[5] = vertex_indices[1] + (num_vertices[0] * num_vertices[1]);
     207              :             vertex_indices[6] = vertex_indices[2] + (num_vertices[0] * num_vertices[1]);
     208              :             vertex_indices[7] = vertex_indices[3] + (num_vertices[0] * num_vertices[1]);
     209              : 
     210              :         ...
     211              :         */
     212              : 
     213           64 :         for (size_t d = 1; d < Dim; ++d) {
     214          152 :             for (size_t i = 0; i < static_cast<unsigned>(1 << d); ++i) {
     215          112 :                 size_t size = 1;
     216          288 :                 for (size_t j = 0; j < d; ++j) {
     217          176 :                     size *= num_vertices[j];
     218              :                 }
     219          112 :                 vertex_indices[i + (1 << d)] = vertex_indices[i] + size;
     220              :             }
     221              :         }
     222              : 
     223           24 :         return vertex_indices;
     224           24 :     }
     225              : 
     226              :     template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
     227              :               typename QuadratureType, typename FieldLHS, typename FieldRHS>
     228              :     KOKKOS_FUNCTION typename FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType,
     229              :                                                 FieldLHS, FieldRHS>::indices_list_t
     230           12 :     FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS, FieldRHS>::
     231              :         getElementMeshVertexNDIndices(
     232              :             const FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS,
     233              :                                      FieldRHS>::indices_t& elementNDIndex) const {
     234           12 :         indices_list_t vertex_nd_indices;
     235              : 
     236           12 :         indices_t smallest_vertex_nd_index = elementNDIndex;
     237              : 
     238              :         // vertex_nd_indices[0] = smallest_vertex_nd_index;
     239              :         // vertex_nd_indices[1] = smallest_vertex_nd_index;
     240              :         // vertex_nd_indices[1][0] += 1;
     241              : 
     242              :         // vertex_nd_indices[2] = vertex_nd_indices[0];
     243              :         // vertex_nd_indices[2][1] += 1;
     244              :         // vertex_nd_indices[3] = vertex_nd_indices[1];
     245              :         // vertex_nd_indices[3][1] += 1;
     246              : 
     247              :         // vertex_nd_indices[4] = vertex_nd_indices[0];
     248              :         // vertex_nd_indices[4][2] += 1;
     249              :         // vertex_nd_indices[5] = vertex_nd_indices[1];
     250              :         // vertex_nd_indices[5][2] += 1;
     251              :         // vertex_nd_indices[6] = vertex_nd_indices[2];
     252              :         // vertex_nd_indices[6][2] += 1;
     253              :         // vertex_nd_indices[7] = vertex_nd_indices[3];
     254              :         // vertex_nd_indices[7][2] += 1;
     255              : 
     256           56 :         for (size_t i = 0; i < (1 << Dim); ++i) {
     257           44 :             vertex_nd_indices[i] = smallest_vertex_nd_index;
     258          136 :             for (size_t j = 0; j < Dim; ++j) {
     259           92 :                 vertex_nd_indices[i][j] += (i >> j) & 1;
     260              :             }
     261              :         }
     262              : 
     263           12 :         return vertex_nd_indices;
     264           12 :     }
     265              : 
     266              :     template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
     267              :               typename QuadratureType, typename FieldLHS, typename FieldRHS>
     268              :     KOKKOS_FUNCTION typename FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType,
     269              :                                                 FieldLHS, FieldRHS>::vertex_points_t
     270           10 :     FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS, FieldRHS>::
     271              :         getElementMeshVertexPoints(
     272              :             const FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS,
     273              :                                      FieldRHS>::indices_t& elementNDIndex) const {
     274           10 :         vertex_points_t vertex_points;
     275              : 
     276              :         // get all the NDIndices for the vertices of this element
     277           10 :         indices_list_t vertex_nd_indices = this->getElementMeshVertexNDIndices(elementNDIndex);
     278              : 
     279              :         // get the coordinates of the vertices of this element
     280           46 :         for (size_t i = 0; i < vertex_nd_indices.dim; ++i) {
     281           36 :             NDIndex<Dim> temp_ndindex;
     282          112 :             for (size_t d = 0; d < Dim; ++d) {
     283           76 :                 temp_ndindex[d]     = Index(vertex_nd_indices[i][d], vertex_nd_indices[i][d]);
     284           76 :                 vertex_points[i][d] = (temp_ndindex[d].first() * this->hr_m[d]) + this->origin_m[d];
     285              :             }
     286              :         }
     287           10 :         return vertex_points;
     288           10 :     }
     289              : 
     290              : }  // namespace ippl
        

Generated by: LCOV version 2.0-1