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