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