IPPL (Independent Parallel Particle Layer)
IPPL
Loading...
Searching...
No Matches
NedelecSpace.hpp
Go to the documentation of this file.
1#include "NedelecSpace.h"
2namespace ippl {
3
4 // NedelecSpace constructor, which calls the FiniteElementSpace constructor,
5 // and decomposes the elements among ranks according to layout.
6 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
7 typename QuadratureType, typename FieldType>
8 NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
9 ::NedelecSpace(UniformCartesian<T, Dim>& mesh, ElementType& ref_element,
10 const QuadratureType& quadrature, const Layout_t& layout)
11 : FiniteElementSpace<T, Dim, getNedelecNumElementDOFs(Dim, Order),
12 ElementType, QuadratureType, FEMVector<T>, FEMVector<T>>
13 (mesh, ref_element, quadrature) {
14 // Assert that the dimension is either 2 or 3.
15 static_assert(Dim >= 2 && Dim <= 3,
16 "The Nedelec Finite Element space only supports 2D and3D meshes");
17
18 // Initialize the elementIndices view
19 initializeElementIndices(layout);
20 }
21
22 // NedelecSpace constructor, which calls the FiniteElementSpace constructor.
23 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
24 typename QuadratureType, typename FieldType>
25 NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
26 ::NedelecSpace(UniformCartesian<T, Dim>& mesh, ElementType& ref_element,
27 const QuadratureType& quadrature)
28 : FiniteElementSpace<T, Dim, getNedelecNumElementDOFs(Dim, Order),
29 ElementType, QuadratureType, FEMVector<T>, FEMVector<T>>
30 (mesh, ref_element, quadrature) {
31
32 // Assert that the dimension is either 2 or 3.
33 static_assert(Dim >= 2 && Dim <= 3,
34 "The Nedelec Finite Element space only supports 2D and 3D meshes");
35 }
36
37 // NedelecSpace initializer, to be made available to the FEMPoissonSolver
38 // such that we can call it from setRhs.
39 // Sets the correct mesh and decomposes the elements among ranks according to layout.
40 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
41 typename QuadratureType, typename FieldType>
42 void NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
43 ::initialize(UniformCartesian<T, Dim>& mesh, const Layout_t& layout) {
44
45 FiniteElementSpace<T, Dim, getNedelecNumElementDOFs(Dim, Order),
46 ElementType, QuadratureType, FEMVector<T>, FEMVector<T>>::setMesh(mesh);
47
48 // Initialize the elementIndices view
49 initializeElementIndices(layout);
50
51 // set the local DOF position vector
52 localDofPositions_m(0)(0) = 0.5;
53 localDofPositions_m(1)(1) = 0.5;
54 localDofPositions_m(2)(0) = 0.5; localDofPositions_m(2)(1) = 1;
55 localDofPositions_m(3)(0) = 1; localDofPositions_m(3)(1) = 0.5;
56 localDofPositions_m(4)(2) = 0.5;
57 localDofPositions_m(5)(0) = 1; localDofPositions_m(5)(2) = 0.5;
58 localDofPositions_m(6)(0) = 1; localDofPositions_m(6)(1) = 1;
59 localDofPositions_m(6)(2) = 0.5;
60 localDofPositions_m(7)(1) = 1; localDofPositions_m(7)(2) = 0.5;
61 localDofPositions_m(8)(0) = 0.5; localDofPositions_m(8)(2) = 1;
62 localDofPositions_m(9)(1) = 0.5; localDofPositions_m(9)(2) = 1;
63 localDofPositions_m(10)(0) = 0.5; localDofPositions_m(10)(1) = 1;
64 localDofPositions_m(10)(2) = 1;
65 localDofPositions_m(11)(0) = 1; localDofPositions_m(11)(1) = 0.5;
66 localDofPositions_m(11)(2) = 1;
67 }
68
69 // Initialize element indices Kokkos View
70 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
71 typename QuadratureType, typename FieldType>
72 void NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
73 ::initializeElementIndices(const Layout_t& layout) {
74
75 layout_m = layout;
76 const auto& ldom = layout.getLocalNDIndex();
77 int npoints = ldom.size();
78 auto first = ldom.first();
79 auto last = ldom.last();
80 ippl::Vector<double, Dim> bounds;
81
82 for (size_t d = 0; d < Dim; ++d) {
83 bounds[d] = this->nr_m[d] - 1;
84 }
85
86 int upperBoundaryPoints = -1;
87
88 Kokkos::View<size_t*> points("ComputeMapping", npoints);
89 Kokkos::parallel_reduce(
90 "ComputePoints", npoints,
91 KOKKOS_CLASS_LAMBDA(const int i, int& local) {
92 int idx = i;
93 indices_t val;
94 bool isBoundary = false;
95 for (unsigned int d = 0; d < Dim; ++d) {
96 int range = last[d] - first[d] + 1;
97 val[d] = first[d] + (idx % range);
98 idx /= range;
99 if (val[d] == bounds[d]) {
100 isBoundary = true;
101 }
102 }
103 points(i) = (!isBoundary) * (this->getElementIndex(val));
104 local += isBoundary;
105 },
106 Kokkos::Sum<int>(upperBoundaryPoints));
107 Kokkos::fence();
108
109 int elementsPerRank = npoints - upperBoundaryPoints;
110 elementIndices = Kokkos::View<size_t*>("i", elementsPerRank);
111 Kokkos::View<size_t> index("index");
112
113 Kokkos::parallel_for(
114 "RemoveNaNs", npoints, KOKKOS_CLASS_LAMBDA(const int i) {
115 if ((points(i) != 0) || (i == 0)) {
116 const size_t idx = Kokkos::atomic_fetch_add(&index(), 1);
117 elementIndices(idx) = points(i);
118 }
119 }
120 );
121 }
122
126
127 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
128 typename QuadratureType, typename FieldType>
129 KOKKOS_FUNCTION size_t NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
130 ::numGlobalDOFs() const {
131
132 size_t num_global_dofs = 0;
133
134 for (size_t d = 0; d < Dim; ++d) {
135 size_t accu = this->nr_m[d]-1;
136 for (size_t d2 = 0; d2 < Dim; ++d2) {
137 if (d == d2) continue;
138 accu *= this->nr_m[d2];
139 }
140 num_global_dofs += accu;
141 }
142
143 return num_global_dofs;
144 }
145
146 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
147 typename QuadratureType, typename FieldType>
148 KOKKOS_FUNCTION size_t NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
149 ::getLocalDOFIndex(const size_t& elementIndex,
150 const size_t& globalDOFIndex) const {
151
152 static_assert(Dim == 2 || Dim == 3, "Dim must be 2 or 3");
153
154 // Get all the global DOFs for the element
155 const Vector<size_t, numElementDOFs> global_dofs =
156 this->NedelecSpace::getGlobalDOFIndices(elementIndex);
157
158 ippl::Vector<size_t, numElementDOFs> dof_mapping;
159 if (Dim == 2) {
160 dof_mapping = {0, 1, 2, 3};
161 } else if (Dim == 3) {
162 dof_mapping = {0, 1, 2, 3, 4, 5, 7, 6, 8, 9, 10, 11};
163 }
164
165 // Find the global DOF in the vector and return the local DOF index
166 // TODO this can be done faster since the global DOFs are sorted
167 for (size_t i = 0; i < dof_mapping.dim; ++i) {
168 if (global_dofs[dof_mapping[i]] == globalDOFIndex) {
169 return dof_mapping[i];
170 }
171 }
172 // it would be good to throw an error in this case
173 // just like the comment in the LagrangeSpace::getLocalDOFIndex()
174 return 0;
175 }
176
177 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
178 typename QuadratureType, typename FieldType>
179 KOKKOS_FUNCTION size_t NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
180 ::getGlobalDOFIndex(const size_t& elementIndex,
181 const size_t& localDOFIndex) const {
182
183 const auto global_dofs = this->NedelecSpace::getGlobalDOFIndices(elementIndex);
184
185 return global_dofs[localDOFIndex];
186 }
187
188 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
189 typename QuadratureType, typename FieldType>
190 KOKKOS_FUNCTION Vector<size_t, NedelecSpace<T, Dim, Order, ElementType,
191 QuadratureType, FieldType>::numElementDOFs>
192 NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
193 ::getLocalDOFIndices() const {
194
195 Vector<size_t, numElementDOFs> localDOFs;
196
197 for (size_t dof = 0; dof < numElementDOFs; ++dof) {
198 localDOFs[dof] = dof;
199 }
200
201 return localDOFs;
202 }
203
204 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
205 typename QuadratureType, typename FieldType>
206 KOKKOS_FUNCTION Vector<size_t, NedelecSpace<T, Dim, Order, ElementType,
207 QuadratureType, FieldType>::numElementDOFs>
208 NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
209 ::getGlobalDOFIndices(const NedelecSpace<T, Dim, Order, ElementType,
210 QuadratureType, FieldType>::indices_t& elementIndex) const {
211
212
213 // These are simply some manual caclualtions that need to be done.
214
215 Vector<size_t, numElementDOFs> globalDOFs(0);
216
217 // Initialize a helper vector v
218 Vector<size_t, Dim> v(1);
219 if constexpr (Dim == 2) {
220 size_t nx = this->nr_m[0];
221 v(1) = 2*nx-1;
222 } else if constexpr (Dim == 3) {
223 size_t nx = this->nr_m[0];
224 size_t ny = this->nr_m[1];
225 v(1) = 2*nx -1;
226 v(2) = 3*nx*ny - nx - ny;
227 }
228
229 // For both 2D and 3D the first few DOF indices are the same
230 size_t nx = this->nr_m[0];
231 globalDOFs(0) = v.dot(elementIndex);
232 globalDOFs(1) = globalDOFs(0) + nx - 1;
233 globalDOFs(2) = globalDOFs(1) + nx;
234 globalDOFs(3) = globalDOFs(1) + 1;
235
236 if constexpr (Dim == 3) {
237 size_t ny = this->nr_m[1];
238
239 globalDOFs(4) = v(2)*elementIndex(2) + 2*nx*ny - nx - ny
240 + elementIndex(1)*nx + elementIndex(0);
241 globalDOFs(5) = globalDOFs(4) + 1;
242 globalDOFs(6) = globalDOFs(4) + nx + 1;
243 globalDOFs(7) = globalDOFs(4) + nx;
244 globalDOFs(8) = globalDOFs(0) + 3*nx*ny - nx - ny;
245 globalDOFs(9) = globalDOFs(8) + nx - 1;
246 globalDOFs(10) = globalDOFs(9) + nx;
247 globalDOFs(11) = globalDOFs(9) + 1;
248 }
249
250
251 return globalDOFs;
252 }
253
254 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
255 typename QuadratureType, typename FieldType>
256 KOKKOS_FUNCTION Vector<size_t, NedelecSpace<T, Dim, Order, ElementType,
257 QuadratureType, FieldType>::numElementDOFs>
258 NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
259 ::getGlobalDOFIndices(const size_t& elementIndex) const {
260
261 indices_t elementPos = this->getElementNDIndex(elementIndex);
262 return getGlobalDOFIndices(elementPos);
263 }
264
265
266
267 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
268 typename QuadratureType, typename FieldType>
269 KOKKOS_FUNCTION Vector<size_t, NedelecSpace<T, Dim, Order, ElementType,
270 QuadratureType, FieldType>::numElementDOFs>
271 NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
272 ::getFEMVectorDOFIndices(NedelecSpace<T, Dim, Order, ElementType,
273 QuadratureType, FieldType>::indices_t elementIndex,
274 NDIndex<Dim> ldom) const {
275
276
277 // This function here is pretty much the same as getGlobalDOFIndices()
278 // the only difference is the domain size and that we have an offset
279 // of the subdomain of the rank to the global one, else it is the same
280
281 Vector<size_t, numElementDOFs> FEMVectorDOFs(0);
282
283 // This corresponds to translating a global element position to one in
284 // the subdomain of the rank. For this we subtract the starting position
285 // the rank subdomain and add the "ghost" hyperplane.
286 elementIndex -= ldom.first();
287 elementIndex += 1;
288
289 indices_t dif(0);
290 dif = ldom.last() - ldom.first();
291 dif += 1 + 2; // plus 1 for last still being in +2 for ghosts.
292
293 // From here on out it is pretty much the same as the
294 // getGlobalDOFIndices() function.
295 Vector<size_t, Dim> v(1);
296 if constexpr (Dim == 2) {
297 size_t nx = dif[0];
298 v(1) = 2*nx-1;
299 } else if constexpr (Dim == 3) {
300 size_t nx = dif[0];
301 size_t ny = dif[1];
302 v(1) = 2*nx -1;
303 v(2) = 3*nx*ny - nx - ny;
304 }
305
306 size_t nx = dif[0];
307 FEMVectorDOFs(0) = v.dot(elementIndex);
308 FEMVectorDOFs(1) = FEMVectorDOFs(0) + nx - 1;
309 FEMVectorDOFs(2) = FEMVectorDOFs(1) + nx;
310 FEMVectorDOFs(3) = FEMVectorDOFs(1) + 1;
311
312 if constexpr (Dim == 3) {
313 size_t ny = dif[1];
314
315 FEMVectorDOFs(4) = v(2)*elementIndex(2) + 2*nx*ny - nx - ny
316 + elementIndex(1)*nx + elementIndex(0);
317 FEMVectorDOFs(5) = FEMVectorDOFs(4) + 1;
318 FEMVectorDOFs(6) = FEMVectorDOFs(4) + nx + 1;
319 FEMVectorDOFs(7) = FEMVectorDOFs(4) + nx;
320 FEMVectorDOFs(8) = FEMVectorDOFs(0) + 3*nx*ny - nx - ny;
321 FEMVectorDOFs(9) = FEMVectorDOFs(8) + nx - 1;
322 FEMVectorDOFs(10) = FEMVectorDOFs(9) + nx;
323 FEMVectorDOFs(11) = FEMVectorDOFs(9) + 1;
324 }
325
326
327 return FEMVectorDOFs;
328 }
329
330 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
331 typename QuadratureType, typename FieldType>
332 KOKKOS_FUNCTION Vector<size_t, NedelecSpace<T, Dim, Order, ElementType,
333 QuadratureType, FieldType>::numElementDOFs>
334 NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
335 ::getFEMVectorDOFIndices(const size_t& elementIndex, NDIndex<Dim> ldom) const {
336
337 // First get the global element position
338 indices_t elementPos = this->getElementNDIndex(elementIndex);
339 return getFEMVectorDOFIndices(elementPos, ldom);
340 }
341
342
343 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
344 typename QuadratureType, typename FieldType>
345 KOKKOS_FUNCTION typename NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>::point_t
346 NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
347 ::getLocalDOFPosition(size_t localDOFIndex) const {
348
349 // Hardcoded center of edges which are stored in the localDofPositions_m
350 // vector. If the DOF position of an edge element actually is the center
351 // of the edge is a different question...
352 return localDofPositions_m(localDOFIndex);
353 }
354
355
356
360 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
361 typename QuadratureType, typename FieldType>
362 template <typename F>
363 FEMVector<T> NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
364 ::evaluateAx(FEMVector<T>& x, F& evalFunction) const {
365 Inform m("");
366
367
368 IpplTimings::TimerRef timerAxInit = IpplTimings::getTimer("Ax init");
369 IpplTimings::startTimer(timerAxInit);
370
371 // create a new field for result, default initialized to zero thanks to
372 // the Kokkos::View
373 FEMVector<T> resultVector = x.template skeletonCopy<T>();
374
375 // List of quadrature weights
376 const Vector<T, QuadratureType::numElementNodes> w =
377 this->quadrature_m.getWeightsForRefElement();
378
379 // List of quadrature nodes
380 const Vector<point_t, QuadratureType::numElementNodes> q =
381 this->quadrature_m.getIntegrationNodesForRefElement();
382
383 // Get the values of the basis functions and their curl at the
384 // quadrature points.
385 Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> curl_b_q;
386 Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> val_b_q;
387 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
388 for (size_t i = 0; i < numElementDOFs; ++i) {
389 curl_b_q[k][i] = this->evaluateRefElementShapeFunctionCurl(i, q[k]);
390 val_b_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
391 }
392 }
393
394 // Get field data and atomic result data,
395 // since it will be added to during the kokkos loop
396 ViewType view = x.getView();
397 AtomicViewType resultView = resultVector.getView();
398
399
400 // Get domain information
401 auto ldom = layout_m.getLocalNDIndex();
402
403 using exec_space = typename Kokkos::View<const size_t*>::execution_space;
404 using policy_type = Kokkos::RangePolicy<exec_space>;
405
406
407 IpplTimings::stopTimer(timerAxInit);
408
409 // Here we assemble the local matrix of an element. In theory this would
410 // have to be done for each element individually, but because we have
411 // that in the case of IPPL all the elements have the same shape we can
412 // also just do it once and then use if all the time.
413 IpplTimings::TimerRef timerAxLocalMatrix = IpplTimings::getTimer("Ax local matrix");
414 IpplTimings::startTimer(timerAxLocalMatrix);
415 Vector<Vector<T, numElementDOFs>, numElementDOFs> A;
416 for (size_t i = 0; i < numElementDOFs; ++i) {
417 for (size_t j = 0; j < numElementDOFs; ++j) {
418 A[i][j] = 0.0;
419 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
420 A[i][j] += w[k] * evalFunction(i, j, curl_b_q[k], val_b_q[k]);
421 }
422 }
423 }
424 IpplTimings::stopTimer(timerAxLocalMatrix);
425
426
427 IpplTimings::TimerRef timerAxLoop = IpplTimings::getTimer("Ax Loop");
428 IpplTimings::startTimer(timerAxLoop);
429
430 // Loop over elements to compute contributions
431 Kokkos::parallel_for(
432 "Loop over elements", policy_type(0, elementIndices.extent(0)),
433 KOKKOS_CLASS_LAMBDA(const size_t index) {
434 const size_t elementIndex = elementIndices(index);
435
436 // Here we now retrieve the global DOF indices and their
437 // position inside of the FEMVector
438 const Vector<size_t, numElementDOFs> global_dofs =
439 this->NedelecSpace::getGlobalDOFIndices(elementIndex);
440
441 const Vector<size_t, numElementDOFs> vectorIndices =
442 this->getFEMVectorDOFIndices(elementIndex, ldom);
443
444
445 // local DOF indices
446 size_t i, j;
447
448 // global DOF n-dimensional indices (Vector of N indices
449 // representing indices in each dimension)
450 size_t I, J;
451
452 for (i = 0; i < numElementDOFs; ++i) {
453 I = global_dofs[i];
454
455 // Skip boundary DOFs (Zero Dirichlet BCs)
456 if (this->isDOFOnBoundary(I)) {
457 continue;
458 }
459
460 for (j = 0; j < numElementDOFs; ++j) {
461 J = global_dofs[j];
462
463 // Skip boundary DOFs (Zero Dirichlet BCs)
464 if (this->isDOFOnBoundary(J)) {
465 continue;
466 }
467
468 resultView(vectorIndices[i]) += A[i][j] * view(vectorIndices[j]);
469 }
470 }
471 }
472 );
473 IpplTimings::stopTimer(timerAxLoop);
474
475 return resultVector;
476
477 }
478
479
480 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
481 typename QuadratureType, typename FieldType>
482 FEMVector<T> NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
483 ::evaluateLoadVector(const FEMVector<NedelecSpace<T, Dim, Order, ElementType,
484 QuadratureType, FieldType>::point_t>& f) const {
485
486 // List of quadrature weights
487 const Vector<T, QuadratureType::numElementNodes> w =
488 this->quadrature_m.getWeightsForRefElement();
489
490 // List of quadrature nodes
491 const Vector<point_t, QuadratureType::numElementNodes> q =
492 this->quadrature_m.getIntegrationNodesForRefElement();
493
494 const indices_t zeroNdIndex = Vector<size_t, Dim>(0);
495
496
497
498 // Evaluate the basis functions for the DOF at the quadrature nodes
499 Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> basis_q;
500 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
501 for (size_t i = 0; i < numElementDOFs; ++i) {
502 basis_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
503 }
504 }
505
506
507 // Get the distance between the quadrature nodes and the DOFs
508 // we assume that the dofs are at the center of an edge, this is then
509 // going to be used to implement a very crude interpolation scheme.
510 Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes>
511 quadratureDOFDistances;
512 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
513 for (size_t i = 0; i < numElementDOFs; ++i) {
514 point_t dofPos = getLocalDOFPosition(i);
515 point_t d = dofPos - q[k];
516 quadratureDOFDistances[k][i] = Kokkos::sqrt(d.dot(d));
517 }
518 }
519
520 // Absolute value of det Phi_K
521 const T absDetDPhi = Kokkos::abs(this->ref_element_m.getDeterminantOfTransformationJacobian(
522 this->getElementMeshVertexPoints(zeroNdIndex)));
523
524 // Get domain information and ghost cells
525 auto ldom = layout_m.getLocalNDIndex();
526
527
528 // Get boundary conditions from field
529 FEMVector<T> resultVector = createFEMVector();
530
531 // Get field data and make it atomic,
532 // since it will be added to during the kokkos loop
533 AtomicViewType atomic_view = resultVector.getView();
534 typename detail::ViewType<point_t, 1>::view_type view = f.getView();
535
536 using exec_space = typename Kokkos::View<const size_t*>::execution_space;
537 using policy_type = Kokkos::RangePolicy<exec_space>;
538
539 // Loop over elements to compute contributions
540 Kokkos::parallel_for(
541 "Loop over elements", policy_type(0, elementIndices.extent(0)),
542 KOKKOS_CLASS_LAMBDA(size_t index) {
543 const size_t elementIndex = elementIndices(index);
544 const Vector<size_t, numElementDOFs> global_dofs =
545 this->NedelecSpace::getGlobalDOFIndices(elementIndex);
546
547 const Vector<size_t, numElementDOFs> vectorIndices =
548 this->getFEMVectorDOFIndices(elementIndex, ldom);
549
550 size_t i;
551
552 for (i = 0; i < numElementDOFs; ++i) {
553 size_t I = global_dofs[i];
554 if (this->isDOFOnBoundary(I)) {
555 continue;
556 }
557
558
559 // calculate the contribution of this element
560 T contrib = 0;
561 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
562 // We now have to interpolate the value of the field
563 // given at the DOF positions to the quadrature point.
564 point_t interpolatedVal(0);
565 T distSum = 0;
566 for (size_t j = 0; j < numElementDOFs; ++j) {
567 // get field index corresponding to this DOF
568
569 // the distance
570 T dist = quadratureDOFDistances[k][j];
571
572 // running variable used for normalization
573 distSum += 1/dist;
574
575 // get field value at DOF and interpolate to q_k
576 interpolatedVal += 1./dist * view(vectorIndices<:j:>);
577 }
578 // here we have to divide it by distSum in order to
579 // normalize it
580 interpolatedVal /= distSum;
581
582 // update contribution
583 contrib += w[k] * basis_q[k][i].dot(interpolatedVal) * absDetDPhi;
584 }
585
586 // add the contribution of the element to the field
587 atomic_view(vectorIndices<:i:>) += contrib;
588
589 }
590 });
591
592 return resultVector;
593 }
594
595 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
596 typename QuadratureType, typename FieldType>
597 template <typename F>
598 FEMVector<T> NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
599 ::evaluateLoadVectorFunctor(const F& f) const {
600
601 // List of quadrature weights
602 const Vector<T, QuadratureType::numElementNodes> w =
603 this->quadrature_m.getWeightsForRefElement();
604
605 // List of quadrature nodes
606 const Vector<point_t, QuadratureType::numElementNodes> q =
607 this->quadrature_m.getIntegrationNodesForRefElement();
608
609 const indices_t zeroNdIndex = Vector<size_t, Dim>(0);
610
611 // Evaluate the basis functions for the DOF at the quadrature nodes
612 Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> basis_q;
613 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
614 for (size_t i = 0; i < numElementDOFs; ++i) {
615 basis_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
616 }
617 }
618
619
620 // Absolute value of det Phi_K
621 const T absDetDPhi = Kokkos::abs(this->ref_element_m.getDeterminantOfTransformationJacobian(
622 this->getElementMeshVertexPoints(zeroNdIndex)));
623
624 // Get domain information and ghost cells
625 auto ldom = layout_m.getLocalNDIndex();
626
627
628 // Get boundary conditions from field
629 FEMVector<T> resultVector = createFEMVector();
630
631 // Get field data and make it atomic,
632 // since it will be added to during the kokkos loop
633 AtomicViewType atomic_view = resultVector.getView();
634
635 using exec_space = typename Kokkos::View<const size_t*>::execution_space;
636 using policy_type = Kokkos::RangePolicy<exec_space>;
637
638
639 // Loop over elements to compute contributions
640 Kokkos::parallel_for(
641 "Loop over elements", policy_type(0, elementIndices.extent(0)),
642 KOKKOS_CLASS_LAMBDA(size_t index) {
643 const size_t elementIndex = elementIndices(index);
644 const Vector<size_t, numElementDOFs> global_dofs =
645 this->NedelecSpace::getGlobalDOFIndices(elementIndex);
646
647 const Vector<size_t, numElementDOFs> vectorIndices =
648 this->getFEMVectorDOFIndices(elementIndex, ldom);
649
650
651 size_t i, I;
652
653 for (i = 0; i < numElementDOFs; ++i) {
654 I = global_dofs[i];
655
656
657 if (this->isDOFOnBoundary(I)) {
658 continue;
659 }
660
661 // calculate the contribution of this element
662 T contrib = 0;
663 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
664 // Get the global position of the quadrature point
665 point_t pos = this->ref_element_m.localToGlobal(
666 this->getElementMeshVertexPoints(this->getElementNDIndex(elementIndex)),
667 q[k]);
668
669 // evaluate the rhs function at this global position
670 point_t interpolatedVal = f(pos);
671
672 // update contribution
673 contrib += w[k] * basis_q[k][i].dot(interpolatedVal) * absDetDPhi;
674 }
675
676 // add the contribution of the element to the vector
677 atomic_view(vectorIndices[i]) += contrib;
678
679 }
680 });
681
682 return resultVector;
683 }
684
685
686
690
691 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
692 typename QuadratureType, typename FieldType>
693 KOKKOS_FUNCTION typename NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>::point_t
694 NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
695 ::evaluateRefElementShapeFunction(const size_t& localDOF,
696 const NedelecSpace<T, Dim, Order, ElementType,
697 QuadratureType, FieldType>::point_t& localPoint) const {
698
699 // Assert that the local vertex index is valid.
700 assert(localDOF < numElementDOFs && "The local vertex index is invalid");
701
702 assert(this->ref_element_m.isPointInRefElement(localPoint)
703 && "Point is not in reference element");
704
705
706
707 // Simply hardcoded
708 point_t result(0);
709 if constexpr (Dim == 2) {
710 T x = localPoint(0);
711 T y = localPoint(1);
712
713 switch (localDOF){
714 case 0: result(0) = 1 - y; break;
715 case 1: result(1) = 1 - x; break;
716 case 2: result(0) = y; break;
717 case 3: result(1) = x; break;
718 }
719 } else if constexpr (Dim == 3) {
720 T x = localPoint(0);
721 T y = localPoint(1);
722 T z = localPoint(2);
723
724 switch (localDOF){
725 case 0: result(0) = y*z - y - z + 1; break;
726 case 1: result(1) = x*z - x - z + 1; break;
727 case 2: result(0) = y*(1 - z); break;
728 case 3: result(1) = x*(1 - z); break;
729 case 4: result(2) = x*y - x - y + 1; break;
730 case 5: result(2) = x*(1 - y); break;
731 case 6: result(2) = x*y; break;
732 case 7: result(2) = y*(1 - x); break;
733 case 8: result(0) = z*(1 - y); break;
734 case 9: result(1) = z*(1 - x); break;
735 case 10: result(0) = y*z; break;
736 case 11: result(1) = x*z; break;
737 }
738 }
739
740
741 return result;
742 }
743
744
745 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
746 typename QuadratureType, typename FieldType>
747 KOKKOS_FUNCTION typename NedelecSpace<T, Dim, Order, ElementType,
748 QuadratureType, FieldType>::point_t
749 NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
750 ::evaluateRefElementShapeFunctionCurl(const size_t& localDOF,
751 const NedelecSpace<T, Dim, Order, ElementType,
752 QuadratureType, FieldType>::point_t& localPoint) const {
753
754 // Hard coded.
755 point_t result(0);
756
757 if constexpr (Dim == 2) {
758 // In case of 2d we would have that the curl would correspond to a
759 // scalar, but in order to keep the interface uniform across all the
760 // dimensions, we still use a 2d vector but only set the first entry
761 // of it. This should lead to no problems, as we later on only will
762 // take the dot product between two of them and therefore should not
763 // run into any problems
764
765 switch (localDOF) {
766 case 0: result(0) = 1; break;
767 case 1: result(0) = -1; break;
768 case 2: result(0) = -1; break;
769 case 3: result(0) = 1; break;
770 }
771 } else {
772 T x = localPoint(0);
773 T y = localPoint(1);
774 T z = localPoint(2);
775
776 switch (localDOF) {
777 case 0: result(0) = 0; result(1) = -1+y; result(2) = 1-z; break;
778 case 1: result(0) = 1-x; result(1) = 0; result(2) = -1+z; break;
779 case 2: result(0) = 0; result(1) = -y; result(2) = -1+z; break;
780 case 3: result(0) = x; result(1) = 0; result(2) = 1-z; break;
781 case 4: result(0) = -1+x; result(1) = 1-y; result(2) = 0; break;
782 case 5: result(0) = -x; result(1) = -1+y; result(2) = 0; break;
783 case 6: result(0) = x; result(1) = -y; result(2) = 0; break;
784 case 7: result(0) = 1-x; result(1) = y; result(2) = 0; break;
785 case 8: result(0) = 0; result(1) = 1-y; result(2) = z; break;
786 case 9: result(0) = -1+x; result(1) = 0; result(2) = -z; break;
787 case 10: result(0) = 0; result(1) = y; result(2) = -z; break;
788 case 11: result(0) = -x; result(1) = 0; result(2) = z; break;
789 }
790 }
791
792 return result;
793
794 }
795
796
800
801 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
802 typename QuadratureType, typename FieldType>
803 FEMVector<T> NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
804 ::createFEMVector() const {
805 // This function will simply call one of the other two depending on the
806 // dimension of the space
807 if constexpr (Dim == 2) {
808 return createFEMVector2d();
809 } else {
810 return createFEMVector3d();
811 }
812 }
813
814
815 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
816 typename QuadratureType, typename FieldType>
817 Kokkos::View<typename NedelecSpace<T, Dim, Order, ElementType, QuadratureType,
818 FieldType>::point_t*>
819 NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>::reconstructToPoints(
820 const Kokkos::View<typename NedelecSpace<T, Dim, Order, ElementType, QuadratureType,
821 FieldType>::point_t*>& positions, const FEMVector<T>& coef) const {
822
823 // The domain information of the subdomain of the MPI rank
824 auto ldom = layout_m.getLocalNDIndex();
825
826 // The domain information of the global domain
827 auto gdom = layout_m.getDomain();
828 indices_t gextent = gdom.last() - gdom.first();
829
830 // The size of the global domain.
831 point_t domainSize = (this->nr_m-1) * this->hr_m;
832
833
834 auto coefView = coef.getView();
835
836 Kokkos::View<point_t*> outView("reconstructed Func values at points", positions.extent(0));
837
838 Kokkos::parallel_for("reconstructToPoints", positions.extent(0),
839 KOKKOS_CLASS_LAMBDA(size_t i) {
840 // get the current position and for it figure out to which
841 // element it belongs
842 point_t pos = positions<:i:>;
843 indices_t elemIdx = ((pos - this->origin_m) / domainSize) * gextent;
844
845
846 // next up we have to handle the case of when a position that
847 // was provided to us lies on an edge at the upper bound of the
848 // local domain, because in this case we have that the above
849 // transformation gives back an element which is somewhat in the
850 // halo. In order to fix this we simply subtract one.
851 for (size_t d = 0; d < Dim; ++d) {
852 if (elemIdx<:d:> >= static_cast<size_t>(ldom.last()<:d:>)) {
853 elemIdx<:d:> -= 1;
854 }
855 }
856
857
858 // get correct indices
859 const Vector<size_t, numElementDOFs> vectorIndices =
860 this->getFEMVectorDOFIndices(elemIdx, ldom);
861
862
863 // figure out position inside of the reference element
864 point_t locPos = pos - (elemIdx * this->hr_m + this->origin_m);
865 locPos /= this->hr_m;
866
867 // because of numerical instabilities it might happen then when
868 // a point is on an edge this becomes marginally larger that 1
869 // or slightly negative which triggers an assertion. So this
870 // simply is to prevent this.
871 for (size_t d = 0; d < Dim; ++d) {
872 locPos<:d:> = Kokkos::min(T(1), locPos<:d:>);
873 locPos<:d:> = Kokkos::max(T(0), locPos<:d:>);
874 }
875
876
877 // interpolate the function value to the position, using the
878 // basis functions.
879 point_t val(0);
880 for (size_t j = 0; j < numElementDOFs; ++j) {
881 point_t funcVal = this->evaluateRefElementShapeFunction(j, locPos);
882 val += funcVal*coefView(vectorIndices<:j:>);
883 }
884 outView(i) = val;
885
886 }
887 );
888
889
890 return outView;
891 }
892
893
897
898 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
899 typename QuadratureType, typename FieldType>
900 template <typename F>
901 T NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>::computeError(
902 const FEMVector<T>& u_h, const F& u_sol) const {
903 if (this->quadrature_m.getOrder() < (2 * Order + 1)) {
904 // throw exception
905 throw IpplException( "NedelecSpace::computeError()",
906 "Order of quadrature rule for error computation should be > 2*p + 1");
907 }
908
909 // List of quadrature weights
910 const Vector<T, QuadratureType::numElementNodes> w =
911 this->quadrature_m.getWeightsForRefElement();
912
913 // List of quadrature nodes
914 const Vector<point_t, QuadratureType::numElementNodes> q =
915 this->quadrature_m.getIntegrationNodesForRefElement();
916
917 // Evaluate the basis functions for the DOF at the quadrature nodes
918 Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> basis_q;
919 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
920 for (size_t i = 0; i < numElementDOFs; ++i) {
921 basis_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
922 }
923 }
924
925 const indices_t zeroNdIndex = Vector<size_t, Dim>(0);
926
927 // Absolute value of det Phi_K
928 const T absDetDPhi = Kokkos::abs(this->ref_element_m
929 .getDeterminantOfTransformationJacobian(this->getElementMeshVertexPoints(zeroNdIndex)));
930
931 // Variable to sum the error to
932 T error = 0;
933
934 // Get domain information and ghost cells
935 auto ldom = layout_m.getLocalNDIndex();
936
937 using exec_space = typename Kokkos::View<const size_t*>::execution_space;
938 using policy_type = Kokkos::RangePolicy<exec_space>;
939
940 auto view = u_h.getView();
941
942
943 // Loop over elements to compute contributions
944 Kokkos::parallel_reduce("Compute error over elements",
945 policy_type(0, elementIndices.extent(0)),
946 KOKKOS_CLASS_LAMBDA(size_t index, double& local) {
947 const size_t elementIndex = elementIndices(index);
948 const Vector<size_t, numElementDOFs> global_dofs =
949 this->NedelecSpace::getGlobalDOFIndices(elementIndex);
950
951 const Vector<size_t, numElementDOFs> vectorIndices =
952 this->getFEMVectorDOFIndices(elementIndex, ldom);
953
954
955 // contribution of this element to the error
956 T contrib = 0;
957 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
958 // Evaluate the analystical solution at the global position
959 // of the quadrature point
960 point_t val_u_sol = u_sol(this->ref_element_m.localToGlobal(
961 this->getElementMeshVertexPoints(this->getElementNDIndex(elementIndex)),
962 q[k]));
963
964 // Here we now reconstruct the solution given the basis
965 // functions.
966 point_t val_u_h = 0;
967 for (size_t j = 0; j < numElementDOFs; ++j) {
968 // get field index corresponding to this DOF
969 val_u_h += basis_q[k][j] * view(vectorIndices[j]);
970 }
971
972 // calculate error and add to sum.
973 point_t dif = (val_u_sol - val_u_h);
974 T x = dif.dot(dif);
975 contrib += w[k] * x * absDetDPhi;
976 }
977 local += contrib;
978 },
979 Kokkos::Sum<double>(error)
980 );
981
982 // MPI reduce
983 T global_error = 0.0;
984 Comm->allreduce(error, global_error, 1, std::plus<T>());
985
986 return Kokkos::sqrt(global_error);
987 }
988
989 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
990 typename QuadratureType, typename FieldType>
991 KOKKOS_FUNCTION bool NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
992 ::isDOFOnBoundary(const size_t& dofIdx) const {
993
994 bool onBoundary = false;
995 if constexpr (Dim == 2) {
996 size_t nx = this->nr_m[0];
997 size_t ny = this->nr_m[1];
998 // South
999 bool sVal = (dofIdx < nx -1);
1000 onBoundary = onBoundary || sVal;
1001 // North
1002 onBoundary = onBoundary || (dofIdx > nx*(ny-1) + ny*(nx-1) - nx);
1003 // West
1004 onBoundary = onBoundary || ((dofIdx >= nx-1) && (dofIdx - (nx-1)) % (2*nx - 1) == 0);
1005 // East
1006 onBoundary = onBoundary || ((dofIdx >= 2*nx-2) && ((dofIdx - 2*nx + 2) % (2*nx - 1) == 0));
1007 }
1008
1009 if constexpr (Dim == 3) {
1010 size_t nx = this->nr_m[0];
1011 size_t ny = this->nr_m[1];
1012 size_t nz = this->nr_m[2];
1013
1014 size_t zOffset = dofIdx / (nx*(ny-1) + ny*(nx-1) + nx*ny);
1015
1016
1017 if (dofIdx - (nx*(ny-1) + ny*(nx-1) + nx*ny)*zOffset >= (nx*(ny-1) + ny*(nx-1))) {
1018 // we are parallel to z axis
1019 // therefore we have halve a cell offset and can never be on the ground or in
1020 // space
1021 size_t f = dofIdx - (nx*(ny-1) + ny*(nx-1) + nx*ny)*zOffset
1022 - (nx*(ny-1) + ny*(nx-1));
1023
1024 size_t yOffset = f / nx;
1025 // South
1026 onBoundary = onBoundary || yOffset == 0;
1027 // North
1028 onBoundary = onBoundary || yOffset == ny-1;
1029
1030 size_t xOffset = f % nx;
1031 // West
1032 onBoundary = onBoundary || xOffset == 0;
1033 // East
1034 onBoundary = onBoundary || xOffset == nx-1;
1035
1036 } else {
1037 // are parallel to one of the other axes
1038 // Ground
1039 onBoundary = onBoundary || zOffset == 0;
1040 // Space
1041 onBoundary = onBoundary || zOffset == nz-1;
1042
1043 size_t f = dofIdx - (nx*(ny-1) + ny*(nx-1) + nx*ny)*zOffset;
1044 size_t yOffset = f / (2*nx - 1);
1045 size_t xOffset = f - (2*nx - 1)*yOffset;
1046
1047 if (xOffset < (nx-1)) {
1048 // we are parallel to the x axis, therefore we cannot
1049 // be on an west or east boundary, but we still can
1050 // be on a north or south boundary
1051
1052 // South
1053 onBoundary = onBoundary || yOffset == 0;
1054 // North
1055 onBoundary = onBoundary || yOffset == ny-1;
1056
1057 } else {
1058 // we are parallel to the y axis, therefore we cannot be
1059 // on a south or north boundary, but we still can be on
1060 // a west or east boundary
1061 if (xOffset >= nx-1) {
1062 xOffset -= (nx-1);
1063 }
1064
1065 // West
1066 onBoundary = onBoundary || xOffset == 0;
1067 // East
1068 onBoundary = onBoundary || xOffset == nx-1;
1069 }
1070 }
1071 }
1072 return onBoundary;
1073 }
1074
1075
1076
1077 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
1078 typename QuadratureType, typename FieldType>
1079 KOKKOS_FUNCTION int NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
1080 ::getBoundarySide(const size_t& dofIdx) const {
1081
1082 if constexpr (Dim == 2) {
1083 size_t nx = this->nr_m[0];
1084 size_t ny = this->nr_m[1];
1085
1086 // South
1087 if (dofIdx < nx -1) return 0;
1088 // West
1089 if ((dofIdx - (nx-1)) % (2*nx - 1) == 0) return 1;
1090 // North
1091 if (dofIdx > nx*(ny-1) + ny*(nx-1) - nx) return 2;
1092 // East
1093 if ((dofIdx >= 2*nx-2) && (dofIdx - 2*nx + 2) % (2*nx - 1) == 0) return 3;
1094
1095 return -1;
1096 }
1097
1098 if constexpr (Dim == 3) {
1099 size_t nx = this->nr_m[0];
1100 size_t ny = this->nr_m[1];
1101 size_t nz = this->nr_m[2];
1102
1103 size_t zOffset = dofIdx / (nx*(ny-1) + ny*(nx-1) + nx*ny);
1104
1105
1106 if (dofIdx - (nx*(ny-1) + ny*(nx-1) + nx*ny)*zOffset >= (nx*(ny-1) + ny*(nx-1))) {
1107 // we are parallel to z axis
1108 // therefore we have halve a cell offset and can never be on the ground or in
1109 // space
1110 size_t f = dofIdx - (nx*(ny-1) + ny*(nx-1) + nx*ny)*zOffset
1111 - (nx*(ny-1) + ny*(nx-1));
1112
1113 size_t yOffset = f / nx;
1114 // South
1115 if (yOffset == 0) return 0;
1116 // North
1117 if (yOffset == ny-1) return 2;
1118
1119 size_t xOffset = f % nx;
1120 // West
1121 if (xOffset == 0) return 1;
1122 // East
1123 if (xOffset == nx-1) return 3;
1124
1125 } else {
1126 // are parallel to one of the other axes
1127 // Ground
1128 if (zOffset == 0) return 4;
1129 // Space
1130 if (zOffset == nz-1) return 5;
1131
1132 size_t f = dofIdx - (nx*(ny-1) + ny*(nx-1) + nx*ny)*zOffset;
1133 size_t yOffset = f / (2*nx - 1);
1134 size_t xOffset = f - (2*nx - 1)*yOffset;
1135
1136 if (xOffset < (nx-1)) {
1137 // we are parallel to the x axis, therefore we cannot
1138 // be on an west or east boundary, but we still can
1139 // be on a north or south boundary
1140
1141 // South
1142 if (yOffset == 0) return 0;
1143 // North
1144 if (yOffset == ny-1) return 2;
1145
1146 } else {
1147 // we are parallel to the y axis, therefore we cannot be
1148 // on a south or north boundary, but we still can be on
1149 // a west or east boundary
1150 if (xOffset >= nx-1) {
1151 xOffset -= (nx-1);
1152 }
1153
1154 // West
1155 if (xOffset == 0) return 1;
1156 // East
1157 if (xOffset == nx-1) return 3;
1158 }
1159 }
1160 return -1;
1161 }
1162
1163 }
1164
1165
1166 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
1167 typename QuadratureType, typename FieldType>
1168 FEMVector<T> NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
1169 ::createFEMVector2d() const{
1170
1171 // Here we will create an empty FEMVector for the case of the domain
1172 // being 2D.
1173 // The largest part of this is going to be the handling of the halo
1174 // cells, more specifically figuring out which entries of the vector are
1175 // part of the sendIdxs and which are part of the recvIdxs. For this we
1176 // loop through all the other domains and for each of them figure out if
1177 // we share a boundary with them. If we share a boundary we then have to
1178 // figure out which entries of the vector are part of this boundary. To
1179 // do this we loop though all the mesh elements which are on this
1180 // boundary and then for each of these elements we take the DOFs which
1181 // should be part of the sendIdxs and the ones which are part of the
1182 // recvIdxs. Currently in this step we have to manually select the
1183 // correct DOFs of the reference element corresponding to that specific
1184 // boundary type (north, south, west, east). This manual selection of
1185 // the correct DOFs might lead to an "out of the blue" feeling on why
1186 // exactly we selected those DOFs, but it should be easily verifyable
1187 // that those are the correct DOFs.
1188 // Also note that we are not exchanging any boundary information over
1189 // corners, test showed that this does not have any impact on the
1190 // correctness.
1191 // For more information regarding the domain decomposition refer to the
1192 // report available at: TODO add reference to report on AMAS website
1193
1194 auto ldom = layout_m.getLocalNDIndex();
1195 auto doms = layout_m.getHostLocalDomains();
1196
1197 // Create the temporaries and so on which will store the MPI
1198 // information.
1199 std::vector<size_t> neighbors;
1200 std::vector< Kokkos::View<size_t*> > sendIdxs;
1201 std::vector< Kokkos::View<size_t*> > recvIdxs;
1202 std::vector< std::vector<size_t> > sendIdxsTemp;
1203 std::vector< std::vector<size_t> > recvIdxsTemp;
1204
1205 // Here we loop thought all the domains to figure out how we are related
1206 // to them and if we have to do any kind of exchange.
1207 size_t myRank = Comm->rank();
1208 for (size_t i = 0; i < doms.extent(0); ++i) {
1209 if (i == myRank) {
1210 // We are looking at ourself
1211 continue;
1212 }
1213 auto odom = doms(i);
1214
1215 // East boundary
1216 if (ldom.last()[0] == odom.first()[0]-1 &&
1217 !(odom.last()[1] < ldom.first()[1] || odom.first()[1] > ldom.last()[1])) {
1218 // Extract the range of the boundary.
1219 int begin = std::max(odom.first()[1], ldom.first()[1]);
1220 int end = std::min(odom.last()[1], ldom.last()[1]);
1221 int pos = ldom.last()[0];
1222
1223 // Add this to the neighbour list.
1224 neighbors.push_back(i);
1225 sendIdxsTemp.push_back(std::vector<size_t>());
1226 recvIdxsTemp.push_back(std::vector<size_t>());
1227 size_t idx = neighbors.size() - 1;
1228
1229 // Add all the halo
1230 indices_t elementPosHalo(0);
1231 elementPosHalo(0) = pos;
1232 indices_t elementPosSend(0);
1233 elementPosSend(0) = pos;
1234 for (int k = begin; k <= end; ++k) {
1235 elementPosHalo(1) = k;
1236 elementPosSend(1) = k;
1237
1238 auto dofIndicesHalo = getFEMVectorDOFIndices(elementPosHalo, ldom);
1239 recvIdxsTemp[idx].push_back(dofIndicesHalo[3]);
1240
1241 auto dofIndicesSend = getFEMVectorDOFIndices(elementPosSend, ldom);
1242 sendIdxsTemp[idx].push_back(dofIndicesSend[0]);
1243 sendIdxsTemp[idx].push_back(dofIndicesSend[1]);
1244 }
1245 // Check if on very north
1246 if (end == layout_m.getDomain().last()[1] || ldom.last()[1] > odom.last()[1]) {
1247 elementPosSend(1) = end;
1248 auto dofIndicesSend = getFEMVectorDOFIndices(elementPosSend, ldom);
1249 // also have to add dof 2
1250 sendIdxsTemp[idx].push_back(dofIndicesSend[2]);
1251 }
1252 }
1253
1254 // West boundary
1255 if (ldom.first()[0] == odom.last()[0]+1 &&
1256 !(odom.last()[1] < ldom.first()[1] || odom.first()[1] > ldom.last()[1])) {
1257 // Extract the range of the boundary.
1258 int begin = std::max(odom.first()[1], ldom.first()[1]);
1259 int end = std::min(odom.last()[1], ldom.last()[1]);
1260 int pos = ldom.first()[0];
1261
1262 // Add this to the neighbour list.
1263 neighbors.push_back(i);
1264 sendIdxsTemp.push_back(std::vector<size_t>());
1265 recvIdxsTemp.push_back(std::vector<size_t>());
1266 size_t idx = neighbors.size() - 1;
1267
1268 // Add all the halo
1269 indices_t elementPosHalo(0);
1270 elementPosHalo(0) = pos-1;
1271 indices_t elementPosSend(0);
1272 elementPosSend(0) = pos;
1273 for (int k = begin; k <= end; ++k) {
1274 elementPosHalo(1) = k;
1275 elementPosSend(1) = k;
1276
1277 auto dofIndicesHalo = getFEMVectorDOFIndices(elementPosHalo, ldom);
1278 recvIdxsTemp[idx].push_back(dofIndicesHalo[0]);
1279 recvIdxsTemp[idx].push_back(dofIndicesHalo[1]);
1280
1281 auto dofIndicesSend = getFEMVectorDOFIndices(elementPosSend, ldom);
1282 sendIdxsTemp[idx].push_back(dofIndicesSend[1]);
1283 }
1284 // Check if on very north
1285 if (end == layout_m.getDomain().last()[1] || odom.last()[1] > ldom.last()[1]) {
1286 elementPosHalo(1) = end;
1287 auto dofIndicesHalo = getFEMVectorDOFIndices(elementPosHalo, ldom);
1288 // also have to add dof 2
1289 recvIdxsTemp[idx].push_back(dofIndicesHalo[2]);
1290 }
1291 }
1292
1293 // North boundary
1294 if (ldom.last()[1] == odom.first()[1]-1 &&
1295 !(odom.last()[0] < ldom.first()[0] || odom.first()[0] > ldom.last()[0])) {
1296 // Extract the range of the boundary.
1297 int begin = std::max(odom.first()[0], ldom.first()[0]);
1298 int end = std::min(odom.last()[0], ldom.last()[0]);
1299 int pos = ldom.last()[1];
1300
1301 // Add this to the neighbour list.
1302 neighbors.push_back(i);
1303 sendIdxsTemp.push_back(std::vector<size_t>());
1304 recvIdxsTemp.push_back(std::vector<size_t>());
1305 size_t idx = neighbors.size() - 1;
1306
1307 // Add all the halo
1308 indices_t elementPosHalo(0);
1309 elementPosHalo(1) = pos;
1310 indices_t elementPosSend(0);
1311 elementPosSend(1) = pos;
1312 for (int k = begin; k <= end; ++k) {
1313 elementPosHalo(0) = k;
1314 elementPosSend(0) = k;
1315
1316 auto dofIndicesHalo = getFEMVectorDOFIndices(elementPosHalo, ldom);
1317 recvIdxsTemp[idx].push_back(dofIndicesHalo[2]);
1318
1319 auto dofIndicesSend = getFEMVectorDOFIndices(elementPosSend, ldom);
1320 sendIdxsTemp[idx].push_back(dofIndicesSend[0]);
1321 sendIdxsTemp[idx].push_back(dofIndicesSend[1]);
1322 }
1323 // Check if on very east
1324 if (end == layout_m.getDomain().last()[0] || ldom.last()[0] > odom.last()[0]) {
1325 elementPosSend(0) = end;
1326 auto dofIndicesSend = getFEMVectorDOFIndices(elementPosSend, ldom);
1327 // also have to add dof 3
1328 sendIdxsTemp[idx].push_back(dofIndicesSend[3]);
1329 }
1330 }
1331
1332 // South boundary
1333 if (ldom.first()[1] == odom.last()[1]+1 &&
1334 !(odom.last()[0] < ldom.first()[0] || odom.first()[0] > ldom.last()[0])) {
1335 // Extract the range of the boundary.
1336 int begin = std::max(odom.first()[0], ldom.first()[0]);
1337 int end = std::min(odom.last()[0], ldom.last()[0]);
1338 int pos = ldom.first()[1];
1339
1340 // Add this to the neighbour list.
1341 neighbors.push_back(i);
1342 sendIdxsTemp.push_back(std::vector<size_t>());
1343 recvIdxsTemp.push_back(std::vector<size_t>());
1344 size_t idx = neighbors.size() - 1;
1345
1346 // Add all the halo
1347 indices_t elementPosHalo(0);
1348 elementPosHalo(1) = pos-1;
1349 indices_t elementPosSend(0);
1350 elementPosSend(1) = pos;
1351 for (int k = begin; k <= end; ++k) {
1352 elementPosHalo(0) = k;
1353 elementPosSend(0) = k;
1354
1355 auto dofIndicesHalo = getFEMVectorDOFIndices(elementPosHalo, ldom);
1356 recvIdxsTemp[idx].push_back(dofIndicesHalo[0]);
1357 recvIdxsTemp[idx].push_back(dofIndicesHalo[1]);
1358
1359 auto dofIndicesSend = getFEMVectorDOFIndices(elementPosSend, ldom);
1360 sendIdxsTemp[idx].push_back(dofIndicesSend[0]);
1361 }
1362 // Check if on very east
1363 if (end == layout_m.getDomain().last()[0] || odom.last()[0] > ldom.last()[0]) {
1364 elementPosHalo(0) = end;
1365 auto dofIndicesHalo = getFEMVectorDOFIndices(elementPosHalo, ldom);
1366 // also have to add dof 3
1367 recvIdxsTemp[idx].push_back(dofIndicesHalo[3]);
1368 }
1369 }
1370 }
1371
1372
1373
1374 // Here we now have to translate the sendIdxsTemp and recvIdxsTemp which
1375 // are std::vectors<std::vector> into the correct list type which
1376 // is std::vector<Kokkos::View>
1377 for (size_t i = 0; i < neighbors.size(); ++i) {
1378 sendIdxs.push_back(Kokkos::View<size_t*>("FEMvector::sendIdxs[" + std::to_string(i) +
1379 "]", sendIdxsTemp[i].size()));
1380 recvIdxs.push_back(Kokkos::View<size_t*>("FEMvector::recvIdxs[" + std::to_string(i) +
1381 "]", recvIdxsTemp[i].size()));
1382 auto sendView = sendIdxs[i];
1383 auto recvView = recvIdxs[i];
1384 auto hSendView = Kokkos::create_mirror_view(sendView);
1385 auto hRecvView = Kokkos::create_mirror_view(recvView);
1386
1387 for (size_t j = 0; j < sendIdxsTemp[i].size(); ++j) {
1388 hSendView(j) = sendIdxsTemp[i][j];
1389 }
1390
1391 for (size_t j = 0; j < recvIdxsTemp[i].size(); ++j) {
1392 hRecvView(j) = recvIdxsTemp[i][j];
1393 }
1394
1395 Kokkos::deep_copy(sendView, hSendView);
1396 Kokkos::deep_copy(recvView, hRecvView);
1397 }
1398
1399
1400
1401 // Now finaly create the FEMVector
1402 indices_t extents(0);
1403 extents = (ldom.last() - ldom.first()) + 3;
1404 size_t nx = extents(0);
1405 size_t ny = extents(1);
1406 size_t n = nx*(ny-1) + ny*(nx-1);
1407 FEMVector<T> vec(n, neighbors, sendIdxs, recvIdxs);
1408
1409 return vec;
1410 }
1411
1412
1413
1414 template <typename T, unsigned Dim, unsigned Order, typename ElementType,
1415 typename QuadratureType, typename FieldType>
1416 FEMVector<T> NedelecSpace<T, Dim, Order, ElementType, QuadratureType, FieldType>
1417 ::createFEMVector3d() const{
1418
1419
1420
1421
1422 // Here we will create an empty FEMVector for the case of the domain
1423 // being 3D.
1424 // It follows the same principle as the 2D case (check comment there).
1425 // The major difference now is that we have more types of boundaries.
1426 // Namely we have 6 directions: west, east, south, north, ground, and
1427 // space. Where west-east is on the x-axis, south-north on the y-axis,
1428 // and ground-space on the z-axis. For this we have 3 major types of
1429 // boundaries, namely "flat" boundaries which are along the coordinate
1430 // axes (your standard west-east, south-north, and ground-space
1431 // exchanges), then we have two different types of diagonal exchanges,
1432 // which we will call "positive" and "negative", we have two types as
1433 // a diagonal exchange always happens over an edge and this edges is
1434 // shared by 4 different ranks, so we have two different diagonal
1435 // exchanges per edge, which we differentiate with "positive" and
1436 // "negative".
1437 // If we now look at these types independently we have that the code
1438 // needed to perform one such exchange is large going to be independent
1439 // of the direction of the exchange (i.e. is the "flat" exchange
1440 // happening over a west-est or a ground-space boundary) the major
1441 // difference is the DOF indices we have to chose for the elements on
1442 // the boundary (check the 2D case for more info).
1443 // We therefore create 3 lambdas for these different types of boundaries
1444 // which we then call with appropriate arguments for the direction of the
1445 // exchange.
1446 // Note that like with the 2D case we do not consider any exchanges over
1447 // corners.
1448 // For more information regarding the domain decomposition refer to the
1449 // report available at: TODO add reference to report on AMAS website
1450
1451 using indices_t = Vector<int, Dim>;
1452
1453 auto ldom = layout_m.getLocalNDIndex();
1454 auto doms = layout_m.getHostLocalDomains();
1455
1456 // Create the temporaries and so on which will store the MPI
1457 // information.
1458 std::vector<size_t> neighbors;
1459 std::vector< Kokkos::View<size_t*> > sendIdxs;
1460 std::vector< Kokkos::View<size_t*> > recvIdxs;
1461 std::vector< std::vector<size_t> > sendIdxsTemp;
1462 std::vector< std::vector<size_t> > recvIdxsTemp;
1463
1464 // The parameters are:
1465 // i: The index of the other dom we are looking at (according to doms).
1466 // a: Along which axis the exchange happens, 0 = x-axis, 1 = y-axis,
1467 // 2 = z-axis.
1468 // f,s: While the exchange happens over the axis "a" we have that
1469 // elements which are part of the boundary are on a plane spanned
1470 // by the other two axes, these other two axes are then given by
1471 // these two variables "f" and "s" (they then also define the
1472 // order in which we go though these axes).
1473 // va,vb: These are the placeholders for the sendIdxs and recvIdxs
1474 // arrays, note that depending on if an exchange happens from,
1475 // e.g. west to east or from east to west the role of which of
1476 // these placeholders stores the sendIdxs and which one stores
1477 // the recvIdxs changes.
1478 // posA, posB: The "a"-axis coordinate of the elements which are part of
1479 // the boundary, we have two of them as the coordinate
1480 // can be different depending on if we are looking at the
1481 // sendIdxs or the recvIdxs.
1482 // idxsA, idxsB: These are going to be the local DOF indices of the
1483 // elements which are part of the boundary and which
1484 // need to be exchanged. Again we have two as the
1485 // indices are going to depend on if we are looking at
1486 // the sendIdxs or the recvIdxs
1487 // adom, bdom: The domain extents of the two domains which are part of
1488 // this exchange.
1489 auto flatBoundaryExchange = [this, &neighbors, &ldom](
1490 size_t i, size_t a, size_t f, size_t s,
1491 std::vector<std::vector<size_t> >& va, std::vector<std::vector<size_t> >& vb,
1492 int posA, int posB,
1493 const std::vector<size_t>& idxsA, const std::vector<size_t>& idxsB,
1494 NDIndex<3>& adom, NDIndex<3>& bdom) {
1495
1496 int beginF = std::max(bdom.first()[f], adom.first()[f]);
1497 int endF = std::min(bdom.last()[f], adom.last()[f]);
1498 int beginS = std::max(bdom.first()[s], adom.first()[s]);
1499 int endS = std::min(bdom.last()[s], adom.last()[s]);
1500
1501 neighbors.push_back(i);
1502 va.push_back(std::vector<size_t>());
1503 vb.push_back(std::vector<size_t>());
1504 size_t idx = neighbors.size() - 1;
1505
1506
1507 indices_t elementPosA(0);
1508 elementPosA(a) = posA;
1509 indices_t elementPosB(0);
1510 elementPosB(a) = posB;
1511
1512 // Here we now have the double loop that goes though all the
1513 // elements spanned by the plane given by the "f" and "s" axis.
1514 for (int k = beginF; k <= endF; ++k) {
1515 elementPosA(f) = k;
1516 elementPosB(f) = k;
1517 for (int l = beginS; l <= endS; ++l) {
1518 elementPosA(s) = l;
1519 elementPosB(s) = l;
1520
1521 auto dofIndicesA = this->getFEMVectorDOFIndices(elementPosA, ldom);
1522 va[idx].push_back(dofIndicesA[idxsA[0]]);
1523 va[idx].push_back(dofIndicesA[idxsA[1]]);
1524
1525 auto dofIndicesB = this->getFEMVectorDOFIndices(elementPosB, ldom);
1526 vb[idx].push_back(dofIndicesB[idxsB[0]]);
1527 vb[idx].push_back(dofIndicesB[idxsB[1]]);
1528 vb[idx].push_back(dofIndicesB[idxsB[2]]);
1529
1530
1531 // We now have reached the end of the first axis and have to
1532 // figure out if we need to add any additional DOFs. If we
1533 // need to add DOFs depends on if we are at the mesh
1534 // boundary and if one of the two domains does not end here
1535 // and "overlaps" the other one.
1536 if (k == endF) {
1537 if (endF == layout_m.getDomain().last()[f] ||
1538 bdom.last()[f] > adom.last()[f]) {
1539 va[idx].push_back(dofIndicesA[idxsA[2]]);
1540 }
1541
1542 if (endF == layout_m.getDomain().last()[f] ||
1543 adom.last()[f] > bdom.last()[f]) {
1544 vb[idx].push_back(dofIndicesB[idxsB[3]]);
1545 vb[idx].push_back(dofIndicesB[idxsB[4]]);
1546 }
1547
1548 // This is a modification to the beginning of the f axis
1549 // we still put it on here as we are guaranteed that
1550 // like this it only gets called once
1551 // call this last, as modifies elementPosA(s)
1552 if (bdom.first()[f] < adom.first()[f]) {
1553 indices_t tmpPos = elementPosA;
1554 tmpPos(f) = beginF-1;
1555 auto dofIndicestmp = this->getFEMVectorDOFIndices(tmpPos, ldom);
1556 va[idx].push_back(dofIndicestmp[idxsA[0]]);
1557 va[idx].push_back(dofIndicestmp[idxsA[1]]);
1558 }
1559 }
1560 }
1561
1562 // We now have reached the end of the second axis and have to
1563 // figure out if we need to add any additional DOFs. If we need
1564 // to add DOFs depends on if we are at the mesh boundary and if
1565 // one of the two domains does not end here and "overlaps" the
1566 // other one.
1567
1568 if (endS == layout_m.getDomain().last()[s] || bdom.last()[s] > adom.last()[s]) {
1569 elementPosA(s) = endS;
1570 auto dofIndicesA = this->getFEMVectorDOFIndices(elementPosA, ldom);
1571 va[idx].push_back(dofIndicesA[idxsA[3]]);
1572 }
1573
1574 if (endS == layout_m.getDomain().last()[s] || adom.last()[s] > bdom.last()[s]) {
1575 elementPosB(s) = endS;
1576 auto dofIndicesB = this->getFEMVectorDOFIndices(elementPosB, ldom);
1577 vb[idx].push_back(dofIndicesB[idxsB[5]]);
1578 vb[idx].push_back(dofIndicesB[idxsB[6]]);
1579 }
1580
1581 // This is a modification to the beginning of the s axis
1582 // we still put it on here as we are guaranteed that
1583 // like this it only gets called once
1584 // call this last, as modifies elementPosA(f);
1585 if (bdom.first()[f] < adom.first()[f]) {
1586 indices_t tmpPos = elementPosA;
1587 tmpPos(s) = beginS-1;
1588 auto dofIndicestmp = this->getFEMVectorDOFIndices(tmpPos, ldom);
1589 va[idx].push_back(dofIndicestmp[idxsA[0]]);
1590 va[idx].push_back(dofIndicestmp[idxsA[1]]);
1591 }
1592 }
1593 // At this point we have reached the end of both axes f and s and
1594 // therefore now have to make one final check.
1595 if ((endF == layout_m.getDomain().last()[f] || adom.last()[f] > bdom.last()[f]) &&
1596 (endS == layout_m.getDomain().last()[s] || adom.last()[s] > bdom.last()[s])) {
1597 elementPosB(f) = endF;
1598 elementPosB(s) = endS;
1599 auto dofIndicesB = this->getFEMVectorDOFIndices(elementPosB, ldom);
1600 vb[idx].push_back(dofIndicesB[idxsB[7]]);
1601 }
1602 };
1603
1604 // The parameters are:
1605 // i: The index of the other dom we are looking at (according to doms).
1606 // a: Along which axis the edge is over which the exchange happens,
1607 // 0 = x-axis, 1 = y-axis, 2 = z-axis.
1608 // f,s: While the exchange happens over the edge along the axis "a" we
1609 // have store in those two the other two axes.
1610 // ao,bo: These are offset variables as certain exchanges require
1611 // offsets to certain values.
1612 // va,vb: These are the placeholders for the sendIdxs and recvIdxs
1613 // arrays, note that depending on if an exchange happens from,
1614 // e.g. west to east or from east to west the role of which of
1615 // these placeholders stores the sendIdxs and which one stores
1616 // the recvIdxs changes.
1617 // idxsA, idxsB: These are going to be the local DOF indices of the
1618 // elements which are part of the boundary and which
1619 // need to be exchanged. Again we have two as the
1620 // indices are going to depend on if we are looking at
1621 // the sendIdxs or the recvIdxs
1622 // odom: The other domain we are exchanging to.
1623 auto negativeDiagonalExchange = [this, &neighbors, &ldom](
1624 size_t i, size_t a, size_t f, size_t s, int ao, int bo,
1625 std::vector<std::vector<size_t> >& va, std::vector<std::vector<size_t> >& vb,
1626 const std::vector<size_t>& idxsA, const std::vector<size_t>& idxsB,
1627 NDIndex<3>& odom) {
1628
1629 neighbors.push_back(i);
1630 va.push_back(std::vector<size_t>());
1631 vb.push_back(std::vector<size_t>());
1632 size_t idx = neighbors.size() - 1;
1633
1634 indices_t elementPosA(0);
1635 elementPosA(f) = ldom.last()[f];
1636 elementPosA(s) = ldom.first()[s] + ao;
1637
1638 indices_t elementPosB(0);
1639 elementPosB(f) = ldom.last()[f];
1640 elementPosB(s) = ldom.first()[s] + bo;
1641
1642 int begin = std::max(odom.first()[a], ldom.first()[a]);
1643 int end = std::min(odom.last()[a], ldom.last()[a]);
1644 // Loop through all the elements along the edge.
1645 for (int k = begin; k <= end; ++k) {
1646 elementPosA(a) = k;
1647 elementPosB(a) = k;
1648
1649 auto dofIndicesA = this->getFEMVectorDOFIndices(elementPosA, ldom);
1650 va[idx].push_back(dofIndicesA[idxsA[0]]);
1651 va[idx].push_back(dofIndicesA[idxsA[1]]);
1652
1653 auto dofIndicesB = this->getFEMVectorDOFIndices(elementPosB, ldom);
1654 vb[idx].push_back(dofIndicesB[idxsB[0]]);
1655 vb[idx].push_back(dofIndicesB[idxsB[1]]);
1656 }
1657 };
1658
1659
1660 // The parameters are:
1661 // i: The index of the other dom we are looking at (according to doms).
1662 // a: Along which axis the edge is over which the exchange happens,
1663 // 0 = x-axis, 1 = y-axis, 2 = z-axis.
1664 // f,s: While the exchange happens over the edge along the axis "a" we
1665 // have store in those two the other two axes.
1666 // va,vb: These are the placeholders for the sendIdxs and recvIdxs
1667 // arrays, note that depending on if an exchange happens from,
1668 // e.g. west to east or from east to west the role of which of
1669 // these placeholders stores the sendIdxs and which one stores
1670 // the recvIdxs changes.
1671 // idxsA, idxsB: These are going to be the local DOF indices of the
1672 // elements which are part of the boundary and which
1673 // need to be exchanged. Again we have two as the
1674 // indices are going to depend on if we are looking at
1675 // the sendIdxs or the recvIdxs
1676 // odom: The other domain we are exchanging to.
1677 auto positiveDiagonalExchange = [this, &neighbors, &ldom](
1678 size_t i, size_t a, size_t f, size_t s,
1679 indices_t posA, indices_t posB,
1680 std::vector<std::vector<size_t> >& va, std::vector<std::vector<size_t> >& vb,
1681 const std::vector<size_t>& idxsA, const std::vector<size_t>& idxsB,
1682 NDIndex<3>& odom) {
1683
1684 neighbors.push_back(i);
1685 va.push_back(std::vector<size_t>());
1686 vb.push_back(std::vector<size_t>());
1687 size_t idx = neighbors.size() - 1;
1688
1689 indices_t elementPosA(0);
1690 elementPosA(f) = posA(f);
1691 elementPosA(s) = posA(s);
1692
1693 indices_t elementPosB(0);
1694 elementPosB(f) = posB(f);
1695 elementPosB(s) = posB(s);
1696
1697 int begin = std::max(odom.first()[a], ldom.first()[a]);
1698 int end = std::min(odom.last()[a], ldom.last()[a]);
1699
1700 for (int k = begin; k <= end; ++k) {
1701 elementPosA(a) = k;
1702 elementPosB(a) = k;
1703
1704 auto dofIndicesA = this->getFEMVectorDOFIndices(elementPosA, ldom);
1705 va[idx].push_back(dofIndicesA[idxsA[0]]);
1706 va[idx].push_back(dofIndicesA[idxsA[1]]);
1707 va[idx].push_back(dofIndicesA[idxsA[2]]);
1708
1709 auto dofIndicesB = this->getFEMVectorDOFIndices(elementPosB, ldom);
1710 vb[idx].push_back(dofIndicesB[idxsB[0]]);
1711 }
1712 };
1713
1714 // After we now have defined the code required for each of the exchanges
1715 // we can look at all the exchanges we need to make and call the lambdas
1716 // with the appropriate parameters.
1717
1718
1719 // Here we loop through all the domains to figure out how we are related
1720 // to them and if we have to do any kind of exchange.
1721 size_t myRank = Comm->rank();
1722 for (size_t i = 0; i < doms.extent(0); ++i) {
1723 if (i == myRank) {
1724 // We are looking at ourself
1725 continue;
1726 }
1727 auto odom = doms(i);
1728
1729 // East boundary
1730 if (ldom.last()[0] == odom.first()[0]-1 &&
1731 !(odom.last()[1] < ldom.first()[1] || odom.first()[1] > ldom.last()[1]) &&
1732 !(odom.last()[2] < ldom.first()[2] || odom.first()[2] > ldom.last()[2])) {
1733
1734 int pos = ldom.last()[0];
1735 flatBoundaryExchange(
1736 i, 0, 1, 2,
1737 recvIdxsTemp, sendIdxsTemp,
1738 pos, pos,
1739 {3,5,6,11}, {0,1,4,2,7,8,9,10},
1740 ldom, odom
1741 );
1742 }
1743
1744 // West boundary
1745 if (ldom.first()[0] == odom.last()[0]+1 &&
1746 !(odom.last()[1] < ldom.first()[1] || odom.first()[1] > ldom.last()[1]) &&
1747 !(odom.last()[2] < ldom.first()[2] || odom.first()[2] > ldom.last()[2])) {
1748
1749 int pos = ldom.first()[0];
1750 flatBoundaryExchange(
1751 i, 0, 1, 2,
1752 sendIdxsTemp, recvIdxsTemp,
1753 pos, pos-1,
1754 {1,4,7,9}, {0,1,4,2,7,8,9,10},
1755 odom, ldom
1756 );
1757 }
1758
1759 // North boundary
1760 if (ldom.last()[1] == odom.first()[1]-1 &&
1761 !(odom.last()[0] < ldom.first()[0] || odom.first()[0] > ldom.last()[0]) &&
1762 !(odom.last()[2] < ldom.first()[2] || odom.first()[2] > ldom.last()[2])) {
1763
1764 int pos = ldom.last()[1];
1765 flatBoundaryExchange(
1766 i, 1, 0, 2,
1767 recvIdxsTemp, sendIdxsTemp,
1768 pos, pos,
1769 {2,7,6,10}, {0,1,4,3,5,8,9,11},
1770 ldom, odom
1771 );
1772 }
1773
1774 // South boundary
1775 if (ldom.first()[1] == odom.last()[1]+1 &&
1776 !(odom.last()[0] < ldom.first()[0] || odom.first()[0] > ldom.last()[0]) &&
1777 !(odom.last()[2] < ldom.first()[2] || odom.first()[2] > ldom.last()[2])) {
1778
1779 int pos = ldom.first()[1];
1780 flatBoundaryExchange(
1781 i, 1, 0, 2,
1782 sendIdxsTemp, recvIdxsTemp,
1783 pos, pos-1,
1784 {0,4,5,8}, {0,1,4,3,5,8,9,11},
1785 odom, ldom
1786 );
1787
1788 }
1789
1790 // Space boundary
1791 if (ldom.last()[2] == odom.first()[2]-1 &&
1792 !(odom.last()[0] < ldom.first()[0] || odom.first()[0] > ldom.last()[0]) &&
1793 !(odom.last()[1] < ldom.first()[1] || odom.first()[1] > ldom.last()[1])) {
1794
1795 int pos = ldom.last()[2];
1796 flatBoundaryExchange(
1797 i, 2, 0, 1,
1798 recvIdxsTemp, sendIdxsTemp,
1799 pos, pos,
1800 {8,9,11,10}, {0,1,4,3,5,2,7,6},
1801 ldom, odom
1802 );
1803 }
1804
1805 // Ground boundary
1806 if (ldom.first()[2] == odom.last()[2]+1 &&
1807 !(odom.last()[0] < ldom.first()[0] || odom.first()[0] > ldom.last()[0]) &&
1808 !(odom.last()[1] < ldom.first()[1] || odom.first()[1] > ldom.last()[1])) {
1809
1810 int pos = ldom.first()[2];
1811 flatBoundaryExchange(
1812 i, 2, 0, 1,
1813 sendIdxsTemp, recvIdxsTemp,
1814 pos, pos-1,
1815 {0,1,3,2}, {0,1,4,3,5,2,7,6},
1816 odom, ldom
1817 );
1818 }
1819
1820
1821
1822 // Next up we handle all the annoying diagonals.
1823 // The negative ones:
1824 // Parallel to y from space to ground, west to east
1825 if (ldom.last()[0] == odom.first()[0]-1 && ldom.first()[2] == odom.last()[2]+1 &&
1826 !(odom.last()[1] < ldom.first()[1] || odom.first()[1] > ldom.last()[1])) {
1827
1828 negativeDiagonalExchange(
1829 i, 1, 0, 2, 0, -1,
1830 sendIdxsTemp, recvIdxsTemp,
1831 {0,1}, {3,5},
1832 odom
1833 );
1834 }
1835
1836 // Parallel to y from ground to space, east to west
1837 if (ldom.first()[0] == odom.last()[0]+1 && ldom.last()[2] == odom.first()[2]-1 &&
1838 !(odom.last()[1] < ldom.first()[1] || odom.first()[1] > ldom.last()[1])) {
1839
1840 negativeDiagonalExchange(
1841 i, 1, 2, 0, -1, 0,
1842 recvIdxsTemp, sendIdxsTemp,
1843 {8,9}, {1,4},
1844 odom
1845 );
1846 }
1847
1848
1849 // Parallel to x from space to ground, south to north
1850 if (ldom.last()[1] == odom.first()[1]-1 && ldom.first()[2] == odom.last()[2]+1 &&
1851 !(odom.last()[0] < ldom.first()[0] || odom.first()[0] > ldom.last()[0])) {
1852 negativeDiagonalExchange(
1853 i, 0, 1, 2, 0, -1,
1854 sendIdxsTemp, recvIdxsTemp,
1855 {0,1}, {2,7},
1856 odom
1857 );
1858 }
1859
1860 // Parallel to x from ground to space, north to south
1861 if (ldom.first()[1] == odom.last()[1]+1 && ldom.last()[2] == odom.first()[2]-1 &&
1862 !(odom.last()[0] < ldom.first()[0] || odom.first()[0] > ldom.last()[0])) {
1863 negativeDiagonalExchange(
1864 i, 0, 2, 1, -1, 0,
1865 recvIdxsTemp, sendIdxsTemp,
1866 {8,9}, {0,4},
1867 odom
1868 );
1869 }
1870
1871
1872 // Parallel to z from west to east, north to south
1873 if (ldom.last()[0] == odom.first()[0]-1 && ldom.first()[1] == odom.last()[1]+1 &&
1874 !(odom.last()[2] < ldom.first()[2] || odom.first()[2] > ldom.last()[2])) {
1875 negativeDiagonalExchange(
1876 i, 2, 0, 1, 0, -1,
1877 sendIdxsTemp, recvIdxsTemp,
1878 {0,4}, {3,5},
1879 odom
1880 );
1881 }
1882
1883 // Parallel to z from east to west, south to north
1884 if (ldom.first()[0] == odom.last()[0]+1 && ldom.last()[1] == odom.first()[1]-1 &&
1885 !(odom.last()[2] < ldom.first()[2] || odom.first()[2] > ldom.last()[2])) {
1886 negativeDiagonalExchange(
1887 i, 2, 1, 0, -1, 0,
1888 recvIdxsTemp, sendIdxsTemp,
1889 {2,7}, {1,4},
1890 odom
1891 );
1892 }
1893
1894
1895
1896 // The positive ones
1897 // Parallel to y from ground to space, west to east
1898 if (ldom.last()[0] == odom.first()[0]-1 && ldom.last()[2] == odom.first()[2]-1 &&
1899 !(odom.last()[1] < ldom.first()[1] || odom.first()[1] > ldom.last()[1])) {
1900 positiveDiagonalExchange(
1901 i, 1, 0, 2,
1902 ldom.last(), ldom.last(),
1903 sendIdxsTemp, recvIdxsTemp,
1904 {0,1,4}, {11},
1905 odom
1906 );
1907 }
1908
1909 // Parallel to y from space to ground, east to west
1910 if (ldom.first()[0] == odom.last()[0]+1 && ldom.first()[2] == odom.last()[2]+1 &&
1911 !(odom.last()[1] < ldom.first()[1] || odom.first()[1] > ldom.last()[1])) {
1912 positiveDiagonalExchange(
1913 i, 1, 0, 2,
1914 ldom.first()-1, ldom.first(),
1915 recvIdxsTemp, sendIdxsTemp,
1916 {0,1,4}, {1},
1917 odom
1918 );
1919 }
1920
1921
1922 // Parallel to x from ground to space, south to north
1923 if (ldom.last()[1] == odom.first()[1]-1 && ldom.last()[2] == odom.first()[2]-1 &&
1924 !(odom.last()[0] < ldom.first()[0] || odom.first()[0] > ldom.last()[0])) {
1925 positiveDiagonalExchange(
1926 i, 0, 1, 2,
1927 ldom.last(), ldom.last(),
1928 sendIdxsTemp, recvIdxsTemp,
1929 {0,1,4}, {10},
1930 odom
1931 );
1932 }
1933
1934 // Parallel to x from space to ground, north to south
1935 if (ldom.first()[1] == odom.last()[1]+1 && ldom.first()[2] == odom.last()[2]+1 &&
1936 !(odom.last()[0] < ldom.first()[0] || odom.first()[0] > ldom.last()[0])) {
1937 positiveDiagonalExchange(
1938 i, 0, 1, 2,
1939 ldom.first()-1, ldom.first(),
1940 recvIdxsTemp, sendIdxsTemp,
1941 {0,1,4}, {0},
1942 odom
1943 );
1944 }
1945
1946
1947 // Parallel to z from west to east, south to north
1948 if (ldom.last()[0] == odom.first()[0]-1 && ldom.last()[1] == odom.first()[1]-1 &&
1949 !(odom.last()[2] < ldom.first()[2] || odom.first()[2] > ldom.last()[2])) {
1950 positiveDiagonalExchange(
1951 i, 2, 0, 1,
1952 ldom.last(), ldom.last(),
1953 sendIdxsTemp, recvIdxsTemp,
1954 {0,1,4}, {6},
1955 odom
1956 );
1957 }
1958
1959 // Parallel to z from east to west, north to south
1960 if (ldom.first()[0] == odom.last()[0]+1 && ldom.first()[1] == odom.last()[1]+1 &&
1961 !(odom.last()[2] < ldom.first()[2] || odom.first()[2] > ldom.last()[2])) {
1962 positiveDiagonalExchange(
1963 i, 2, 0, 1,
1964 ldom.first()-1, ldom.first(),
1965 recvIdxsTemp, sendIdxsTemp,
1966 {0,1,4}, {4},
1967 odom
1968 );
1969 }
1970
1971 }
1972
1973
1974
1975
1976 // Here we now have to translate the sendIdxsTemp and recvIdxsTemp which
1977 // are std::vectors<std::vector> into the correct list type which
1978 // is std::vector<Kokkos::View>
1979 for (size_t i = 0; i < neighbors.size(); ++i) {
1980 sendIdxs.push_back(Kokkos::View<size_t*>("FEMvector::sendIdxs[" + std::to_string(i) +
1981 "]", sendIdxsTemp[i].size()));
1982 recvIdxs.push_back(Kokkos::View<size_t*>("FEMvector::recvIdxs[" + std::to_string(i) +
1983 "]", recvIdxsTemp[i].size()));
1984 auto sendView = sendIdxs[i];
1985 auto recvView = recvIdxs[i];
1986 auto hSendView = Kokkos::create_mirror_view(sendView);
1987 auto hRecvView = Kokkos::create_mirror_view(recvView);
1988
1989 for (size_t j = 0; j < sendIdxsTemp[i].size(); ++j) {
1990 hSendView(j) = sendIdxsTemp[i][j];
1991 }
1992
1993 for (size_t j = 0; j < recvIdxsTemp[i].size(); ++j) {
1994 hRecvView(j) = recvIdxsTemp[i][j];
1995 }
1996
1997 Kokkos::deep_copy(sendView, hSendView);
1998 Kokkos::deep_copy(recvView, hRecvView);
1999 }
2000
2001
2002
2003 // Now finaly create the FEMVector
2004 indices_t extents(0);
2005 extents = (ldom.last() - ldom.first()) + 3;
2006 size_t nx = extents(0);
2007 size_t ny = extents(1);
2008 size_t nz = extents(2);
2009 size_t n = (nz-1)*(nx*(ny-1) + ny*(nx-1) + nx*ny) + nx*(ny-1) + ny*(nx-1);
2010 FEMVector<T> vec(n, neighbors, sendIdxs, recvIdxs);
2011
2012 return vec;
2013 }
2014
2015
2016
2017} // namespace ippl
Dim
constexpr unsigned Dim
Definition BumponTailInstability.cpp:22
getNedelecNumElementDOFs
constexpr unsigned getNedelecNumElementDOFs(unsigned Dim, unsigned Order)
Definition NedelecSpace.h:13
first
constexpr KOKKOS_INLINE_FUNCTION auto first()
Definition AbsorbingBC.h:10
ippl
Definition Archive.h:20
ippl::Comm
std::unique_ptr< mpi::Communicator > Comm
Definition Ippl.h:22
ippl::FEMVector
1D vector used in the context of FEM.
Definition FEMVector.h:32
ippl::FEMVector::getView
const Kokkos::View< T * > & getView() const
Get underlying data view.
Definition FEMVector.hpp:223
ippl::FiniteElementSpace< T, Dim, getNedelecNumElementDOFs(Dim, Order), ElementType, QuadratureType, FEMVector< T >, FEMVector< T > >::FiniteElementSpace
FiniteElementSpace(UniformCartesian< T, Dim > &mesh, ElementType &ref_element, const QuadratureType &quadrature)
Definition FiniteElementSpace.hpp:6
ippl::NedelecSpace
A class representing a Nedelec space for finite element methods on a structured, rectilinear grid.
Definition NedelecSpace.h:35
ippl::NedelecSpace::ViewType
detail::ViewType< T, 1 >::view_type ViewType
Definition NedelecSpace.h:66
ippl::NedelecSpace::getFEMVectorDOFIndices
KOKKOS_FUNCTION Vector< size_t, numElementDOFs > getFEMVectorDOFIndices(const size_t &elementIndex, NDIndex< Dim > ldom) const
Get the DOF indices (vector of indices) corresponding to the position inside the FEMVector of an elem...
Definition NedelecSpace.hpp:335
ippl::NedelecSpace::computeError
T computeError(const FEMVector< T > &u_h, const F &u_sol) const
Error norm computations ///////////////////////////////////////////.
Definition NedelecSpace.hpp:901
ippl::NedelecSpace::elementIndices
Kokkos::View< size_t * > elementIndices
Stores which elements (squares or cubes) belong to the current MPI rank.
Definition NedelecSpace.h:407
ippl::NedelecSpace::initializeElementIndices
void initializeElementIndices(const Layout_t &layout)
Initialize a Kokkos view containing the element indices.
Definition NedelecSpace.hpp:73
ippl::NedelecSpace::Layout_t
FieldLayout< Dim > Layout_t
Definition NedelecSpace.h:63
ippl::NedelecSpace::point_t
FiniteElementSpace< T, Dim, numElementDOFs, ElementType, QuadratureType, FEMVector< T >, FEMVector< T > >::point_t point_t
Definition NedelecSpace.h:57
ippl::NedelecSpace::getLocalDOFPosition
KOKKOS_FUNCTION point_t getLocalDOFPosition(size_t localDOFIndex) const
Get the cartesion position of a local DOF in the reference element.
Definition NedelecSpace.hpp:347
ippl::NedelecSpace::evaluateRefElementShapeFunctionCurl
KOKKOS_FUNCTION point_t evaluateRefElementShapeFunctionCurl(const size_t &localDOF, const point_t &localPoint) const
Evaluate the curl of the shape function of a local degree of of freedom at ta given point in the refe...
Definition NedelecSpace.hpp:750
ippl::NedelecSpace::createFEMVector3d
FEMVector< T > createFEMVector3d() const
Implementation of the NedelecSpace::createFEMVector function for 3d.
Definition NedelecSpace.hpp:1417
ippl::NedelecSpace::createFEMVector
FEMVector< T > createFEMVector() const
FEMVector conversion and creation//////////////////////////////////.
Definition NedelecSpace.hpp:804
ippl::NedelecSpace::indices_t
FiniteElementSpace< T, Dim, numElementDOFs, ElementType, QuadratureType, FEMVector< T >, FEMVector< T > >::indices_t indices_t
Definition NedelecSpace.h:53
ippl::NedelecSpace::layout_m
Layout_t layout_m
The layout of the MPI ranks over the mesh.
Definition NedelecSpace.h:425
ippl::NedelecSpace::isDOFOnBoundary
KOKKOS_FUNCTION bool isDOFOnBoundary(const size_t &dofIdx) const
Check if a DOF is on the boundary of the mesh.
Definition NedelecSpace.hpp:992
ippl::NedelecSpace::evaluateRefElementShapeFunction
KOKKOS_FUNCTION point_t evaluateRefElementShapeFunction(const size_t &localDOF, const point_t &localPoint) const
Basis functions and gradients /////////////////////////////////////.
Definition NedelecSpace.hpp:695
ippl::NedelecSpace::getGlobalDOFIndices
KOKKOS_FUNCTION Vector< size_t, numElementDOFs > getGlobalDOFIndices(const size_t &elementIndex) const override
Get the global DOF indices (vector of global DOF indices) of an element.
Definition NedelecSpace.hpp:259
ippl::NedelecSpace::NedelecSpace
NedelecSpace(UniformCartesian< T, Dim > &mesh, ElementType &ref_element, const QuadratureType &quadrature, const Layout_t &layout)
Construct a new NedelecSpace object.
Definition NedelecSpace.hpp:9
ippl::NedelecSpace::AtomicViewType
detail::ViewType< T, 1, Kokkos::MemoryTraits< Kokkos::Atomic > >::view_type AtomicViewType
Definition NedelecSpace.h:68
ippl::NedelecSpace::numElementDOFs
static constexpr unsigned numElementDOFs
Definition NedelecSpace.h:38
ippl::NedelecSpace::reconstructToPoints
Kokkos::View< point_t * > reconstructToPoints(const Kokkos::View< point_t * > &positions, const FEMVector< T > &coef) const
Reconstructs function values at arbitrary points in the mesh given the Nedelec DOF coefficients.
Definition NedelecSpace.hpp:819
ippl::NedelecSpace::createFEMVector2d
FEMVector< T > createFEMVector2d() const
Implementation of the NedelecSpace::createFEMVector function for 2d.
Definition NedelecSpace.hpp:1169
ippl::NedelecSpace::localDofPositions_m
Vector< point_t, 12 > localDofPositions_m
Stores the positions of the local Degrees of Freedoms on the reference elements.
Definition NedelecSpace.h:416
ippl::FieldLayout::getLocalNDIndex
const NDIndex_t & getLocalNDIndex() const
Definition FieldLayout.hpp:116
ippl::NDIndex::last
KOKKOS_INLINE_FUNCTION Vector< int, Dim > last() const
Definition NDIndex.hpp:178
ippl::NDIndex::size
KOKKOS_INLINE_FUNCTION unsigned size() const noexcept
Definition NDIndex.hpp:43
ippl::NDIndex::first
KOKKOS_INLINE_FUNCTION Vector< int, Dim > first() const
Definition NDIndex.hpp:170
ippl::Vector::dim
static constexpr unsigned dim
Definition Vector.h:26
ippl::Vector::dot
KOKKOS_INLINE_FUNCTION T dot(const Vector< T, Dim > &rhs) const
Definition Vector.hpp:182
ippl::detail::ViewType::view_type
Kokkos::View< typename NPtr< T, Dim >::type, Properties... > view_type
Definition ViewTypes.h:45
IpplTimings::TimerRef
Timing::TimerRef TimerRef
Definition IpplTimings.h:144
IpplTimings::getTimer
static TimerRef getTimer(const char *nm)
Definition IpplTimings.h:150
IpplTimings::stopTimer
static void stopTimer(TimerRef t)
Definition IpplTimings.h:156
IpplTimings::startTimer
static void startTimer(TimerRef t)
Definition IpplTimings.h:153