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 446752 : 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 446752 : size_t index = vertex_index;
71 :
72 : // Create a vector to store the vertex indices in each dimension for the corresponding
73 : // vertex.
74 446752 : indices_t vertex_indices;
75 :
76 : // This is the number of vertices in each dimension.
77 446752 : 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 446752 : size_t remaining_number_of_vertices = 1;
83 1778720 : for (const size_t num_vertices : vertices_per_dim) {
84 1331968 : remaining_number_of_vertices *= num_vertices;
85 : }
86 :
87 1778720 : for (int d = Dim - 1; d >= 0; --d) {
88 1331968 : remaining_number_of_vertices /= vertices_per_dim[d];
89 1331968 : vertex_indices[d] = index / remaining_number_of_vertices;
90 1331968 : index -= vertex_indices[d] * remaining_number_of_vertices;
91 : }
92 :
93 446752 : return vertex_indices;
94 446752 : };
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 1268 : 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 1268 : size_t index = element_index;
124 :
125 : // Create a vector to store the element indices in each dimension for the corresponding
126 : // element.
127 1268 : 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 1268 : 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 1268 : size_t remaining_number_of_cells = 1;
137 4758 : for (const size_t num_cells : cells_per_dim) {
138 3490 : remaining_number_of_cells *= num_cells;
139 : }
140 :
141 4758 : for (int d = Dim - 1; d >= 0; --d) {
142 3490 : remaining_number_of_cells /= cells_per_dim[d];
143 3490 : element_nd_index[d] = (index / remaining_number_of_cells);
144 3490 : index -= (element_nd_index[d]) * remaining_number_of_cells;
145 : }
146 :
147 1268 : return element_nd_index;
148 1268 : }
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 30 : 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 30 : indices_list_t vertex_nd_indices;
235 :
236 30 : 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 174 : for (size_t i = 0; i < (1 << Dim); ++i) {
257 144 : vertex_nd_indices[i] = smallest_vertex_nd_index;
258 496 : for (size_t j = 0; j < Dim; ++j) {
259 352 : vertex_nd_indices[i][j] += (i >> j) & 1;
260 : }
261 : }
262 :
263 30 : return vertex_nd_indices;
264 30 : }
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 28 : 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 28 : vertex_points_t vertex_points;
275 :
276 : // get all the NDIndices for the vertices of this element
277 28 : indices_list_t vertex_nd_indices = this->getElementMeshVertexNDIndices(elementNDIndex);
278 :
279 : // get the coordinates of the vertices of this element
280 164 : for (size_t i = 0; i < vertex_nd_indices.dim; ++i) {
281 136 : NDIndex<Dim> temp_ndindex;
282 472 : for (size_t d = 0; d < Dim; ++d) {
283 336 : temp_ndindex[d] = Index(vertex_nd_indices[i][d], vertex_nd_indices[i][d]);
284 336 : vertex_points[i][d] = (temp_ndindex[d].first() * this->hr_m[d]) + this->origin_m[d];
285 : }
286 : }
287 28 : return vertex_points;
288 28 : }
289 :
290 : } // namespace ippl
|
