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 396 : FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS,
6 : FieldRHS>::FiniteElementSpace(UniformCartesian<T, Dim>& mesh,
7 : ElementType& ref_element,
8 : const QuadratureType& quadrature)
9 396 : : mesh_m(mesh)
10 : , ref_element_m(ref_element)
11 396 : , quadrature_m(quadrature) {
12 : assert(mesh.Dimension == Dim && "Mesh dimension does not match the dimension of the space");
13 :
14 396 : nr_m = mesh_m.getGridsize();
15 396 : hr_m = mesh_m.getMeshSpacing();
16 396 : 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 396 : }
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 24 : KOKKOS_FUNCTION size_t FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType,
44 : FieldLHS, FieldRHS>::numElements() const {
45 24 : Vector<size_t, Dim> cells_per_dim = nr_m - 1u;
46 :
47 24 : size_t num_elements = 1;
48 84 : for (size_t d = 0; d < Dim; ++d) {
49 60 : num_elements *= cells_per_dim[d];
50 : }
51 :
52 24 : return num_elements;
53 24 : }
54 :
55 : template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
56 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
57 : KOKKOS_FUNCTION size_t
58 20 : FiniteElementSpace<T, Dim, NumElementDOFs, ElementType, QuadratureType, FieldLHS,
59 : FieldRHS>::numElementsInDim(const size_t& dim) const {
60 20 : 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 446368 : 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 446368 : size_t index = vertex_index;
71 :
72 : // Create a vector to store the vertex indices in each dimension for the corresponding
73 : // vertex.
74 446368 : indices_t vertex_indices;
75 :
76 : // This is the number of vertices in each dimension.
77 446368 : 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 446368 : size_t remaining_number_of_vertices = 1;
83 1777696 : for (const size_t num_vertices : vertices_per_dim) {
84 1331328 : remaining_number_of_vertices *= num_vertices;
85 : }
86 :
87 1777696 : for (int d = Dim - 1; d >= 0; --d) {
88 1331328 : remaining_number_of_vertices /= vertices_per_dim[d];
89 1331328 : vertex_indices[d] = index / remaining_number_of_vertices;
90 1331328 : index -= vertex_indices[d] * remaining_number_of_vertices;
91 : }
92 :
93 446368 : return vertex_indices;
94 446368 : };
95 :
96 : template <typename T, unsigned Dim, unsigned NumElementDOFs, typename ElementType,
97 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
98 : KOKKOS_FUNCTION size_t
99 320 : 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 320 : ippl::Vector<size_t, Dim> vec(1);
105 896 : for (size_t d = 1; d < dim; ++d) {
106 1408 : for (size_t d2 = d; d2 < Dim; ++d2) {
107 832 : vec[d2] *= nr_m[d - 1];
108 : }
109 : }
110 :
111 : // return the dot product between the vertex ndindex and vec.
112 320 : return vertexNDIndex.dot(vec);
113 320 : }
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 1968 : 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 1968 : size_t index = element_index;
124 :
125 : // Create a vector to store the element indices in each dimension for the corresponding
126 : // element.
127 1968 : 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 1968 : 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 1968 : size_t remaining_number_of_cells = 1;
137 7452 : for (const size_t num_cells : cells_per_dim) {
138 5484 : remaining_number_of_cells *= num_cells;
139 : }
140 :
141 7452 : for (int d = Dim - 1; d >= 0; --d) {
142 5484 : remaining_number_of_cells /= cells_per_dim[d];
143 5484 : element_nd_index[d] = (index / remaining_number_of_cells);
144 5484 : index -= (element_nd_index[d]) * remaining_number_of_cells;
145 : }
146 :
147 1968 : return element_nd_index;
148 1968 : }
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 7866 : 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 7866 : 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 7866 : Vector<size_t, Dim> cells_per_dim = nr_m - 1;
163 :
164 7866 : size_t remaining_number_of_cells = 1;
165 :
166 29610 : for (unsigned int d = 0; d < Dim; ++d) {
167 21744 : element_index += ndindex[d] * remaining_number_of_cells;
168 21744 : remaining_number_of_cells *= cells_per_dim[d];
169 : }
170 :
171 7866 : return element_index;
172 7866 : }
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 48 : 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 48 : const Vector<size_t, Dim> num_vertices = nr_m;
183 :
184 48 : size_t smallest_vertex_index = 0;
185 176 : for (int d = Dim - 1; d >= 0; --d) {
186 128 : size_t temp_index = element_nd_index[d];
187 240 : for (int i = d; i >= 1; --i) {
188 112 : temp_index *= num_vertices[i];
189 : }
190 128 : smallest_vertex_index += temp_index;
191 : }
192 :
193 : // Vector to store the vertex indices for the element
194 48 : vertex_indices_t vertex_indices;
195 48 : vertex_indices[0] = smallest_vertex_index;
196 48 : 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 128 : for (size_t d = 1; d < Dim; ++d) {
214 304 : for (size_t i = 0; i < static_cast<unsigned>(1 << d); ++i) {
215 224 : size_t size = 1;
216 576 : for (size_t j = 0; j < d; ++j) {
217 352 : size *= num_vertices[j];
218 : }
219 224 : vertex_indices[i + (1 << d)] = vertex_indices[i] + size;
220 : }
221 : }
222 :
223 48 : return vertex_indices;
224 48 : }
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 56 : 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 56 : indices_list_t vertex_nd_indices;
235 :
236 56 : 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 328 : for (size_t i = 0; i < (1 << Dim); ++i) {
257 272 : vertex_nd_indices[i] = smallest_vertex_nd_index;
258 944 : for (size_t j = 0; j < Dim; ++j) {
259 672 : vertex_nd_indices[i][j] += (i >> j) & 1;
260 : }
261 : }
262 :
263 56 : return vertex_nd_indices;
264 56 : }
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 52 : 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 52 : vertex_points_t vertex_points;
275 :
276 : // get all the NDIndices for the vertices of this element
277 52 : indices_list_t vertex_nd_indices = this->getElementMeshVertexNDIndices(elementNDIndex);
278 :
279 : // get the coordinates of the vertices of this element
280 308 : for (size_t i = 0; i < vertex_nd_indices.dim; ++i) {
281 256 : NDIndex<Dim> temp_ndindex;
282 896 : for (size_t d = 0; d < Dim; ++d) {
283 640 : temp_ndindex[d] = Index(vertex_nd_indices[i][d], vertex_nd_indices[i][d]);
284 640 : vertex_points[i][d] = (temp_ndindex[d].first() * this->hr_m[d]) + this->origin_m[d];
285 : }
286 : }
287 52 : return vertex_points;
288 52 : }
289 :
290 : } // namespace ippl
|
