Line data Source code
1 : // Class LagrangeSpace
2 : // This is the LagrangeSpace class. It is a class representing a Lagrange space
3 : // for finite element methods on a structured grid.
4 :
5 : #ifndef IPPL_LAGRANGESPACE_H
6 : #define IPPL_LAGRANGESPACE_H
7 :
8 : #include <cmath>
9 :
10 : #include "FEM/FiniteElementSpace.h"
11 :
12 : constexpr unsigned getLagrangeNumElementDOFs(unsigned Dim, unsigned Order) {
13 : // needs to be constexpr pow function to work at compile time. Kokkos::pow doesn't work.
14 : return static_cast<unsigned>(power(static_cast<int>(Order + 1), static_cast<int>(Dim)));
15 : }
16 :
17 : namespace ippl {
18 :
19 : /**
20 : * @brief A class representing a Lagrange space for finite element methods on a structured,
21 : * rectilinear grid.
22 : *
23 : * @tparam T The floating point number type of the field values
24 : * @tparam Dim The dimension of the mesh
25 : * @tparam Order The order of the Lagrange space
26 : * @tparam QuadratureType The type of the quadrature rule
27 : * @tparam FieldLHS The type of the left hand side field
28 : * @tparam FieldRHS The type of the right hand side field
29 : */
30 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
31 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
32 : // requires IsQuadrature<QuadratureType>
33 : class LagrangeSpace
34 : : public FiniteElementSpace<T, Dim, getLagrangeNumElementDOFs(Dim, Order), ElementType,
35 : QuadratureType, FieldLHS, FieldRHS> {
36 : public:
37 : // The number of degrees of freedom per element
38 : static constexpr unsigned numElementDOFs = getLagrangeNumElementDOFs(Dim, Order);
39 :
40 : // The dimension of the mesh
41 : static constexpr unsigned dim = FiniteElementSpace<T, Dim, numElementDOFs, ElementType,
42 : QuadratureType, FieldLHS, FieldRHS>::dim;
43 :
44 : // The order of the Lagrange space
45 : static constexpr unsigned order = Order;
46 :
47 : // The number of mesh vertices per element
48 : static constexpr unsigned numElementVertices =
49 : FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType, FieldLHS,
50 : FieldRHS>::numElementVertices;
51 :
52 : // A vector with the position of the element in the mesh in each dimension
53 : typedef typename FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType,
54 : FieldLHS, FieldRHS>::indices_t indices_t;
55 :
56 : // A point in the global coordinate system
57 : typedef typename FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType,
58 : FieldLHS, FieldRHS>::point_t point_t;
59 :
60 : typedef typename FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType,
61 : FieldLHS, FieldRHS>::vertex_points_t vertex_points_t;
62 :
63 : // Field layout type for domain decomposition info
64 : typedef FieldLayout<Dim> Layout_t;
65 :
66 : // View types
67 : typedef typename detail::ViewType<T, Dim>::view_type ViewType;
68 : typedef typename detail::ViewType<T, Dim, Kokkos::MemoryTraits<Kokkos::Atomic>>::view_type
69 : AtomicViewType;
70 :
71 : ///////////////////////////////////////////////////////////////////////
72 : // Constructors ///////////////////////////////////////////////////////
73 : ///////////////////////////////////////////////////////////////////////
74 :
75 : /**
76 : * @brief Construct a new LagrangeSpace object
77 : *
78 : * @param mesh Reference to the mesh
79 : * @param ref_element Reference to the reference element
80 : * @param quadrature Reference to the quadrature rule
81 : * @param layout Reference to the field layout
82 : */
83 : LagrangeSpace(UniformCartesian<T, Dim>& mesh, ElementType& ref_element,
84 : const QuadratureType& quadrature, const Layout_t& layout);
85 :
86 : /**
87 : * @brief Construct a new LagrangeSpace object (without layout)
88 : * This constructor is made to work with the default constructor in
89 : * FEMPoissonSolver.h such that it is compatible with alpine.
90 : *
91 : * @param mesh Reference to the mesh
92 : * @param ref_element Reference to the reference element
93 : * @param quadrature Reference to the quadrature rule
94 : */
95 : LagrangeSpace(UniformCartesian<T, Dim>& mesh, ElementType& ref_element,
96 : const QuadratureType& quadrature);
97 :
98 : /**
99 : * @brief Initialize a LagrangeSpace object created with the default constructor
100 : *
101 : * @param mesh Reference to the mesh
102 : * @param layout Reference to the field layout
103 : */
104 : void initialize(UniformCartesian<T, Dim>& mesh, const Layout_t& layout);
105 :
106 : ///////////////////////////////////////////////////////////////////////
107 : /**
108 : * @brief Initialize a Kokkos view containing the element indices.
109 : * This distributes the elements among MPI ranks.
110 : */
111 : void initializeElementIndices(const Layout_t& layout);
112 :
113 : /// Degree of Freedom operations //////////////////////////////////////
114 : ///////////////////////////////////////////////////////////////////////
115 :
116 : /**
117 : * @brief Get the number of global degrees of freedom in the space
118 : *
119 : * @return size_t - unsigned integer number of global degrees of freedom
120 : */
121 : KOKKOS_FUNCTION size_t numGlobalDOFs() const override;
122 :
123 : /**
124 : * @brief Get the elements local DOF from the element index and global DOF
125 : * index
126 : *
127 : * @param elementIndex size_t - The index of the element
128 : * @param globalDOFIndex size_t - The global DOF index
129 : *
130 : * @return size_t - The local DOF index
131 : */
132 : KOKKOS_FUNCTION size_t getLocalDOFIndex(const size_t& elementIndex,
133 : const size_t& globalDOFIndex) const override;
134 :
135 : /**
136 : * @brief Get the global DOF index from the element index and local DOF
137 : *
138 : * @param elementIndex size_t - The index of the element
139 : * @param localDOFIndex size_t - The local DOF index
140 : *
141 : * @return size_t - The global DOF index
142 : */
143 : KOKKOS_FUNCTION size_t getGlobalDOFIndex(const size_t& elementIndex,
144 : const size_t& localDOFIndex) const override;
145 :
146 : /**
147 : * @brief Get the local DOF indices (vector of local DOF indices)
148 : * They are independent of the specific element because it only depends on
149 : * the reference element type
150 : *
151 : * @return Vector<size_t, NumElementDOFs> - The local DOF indices
152 : */
153 : KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getLocalDOFIndices() const override;
154 :
155 : /**
156 : * @brief Get the global DOF indices (vector of global DOF indices) of an element
157 : *
158 : * @param elementIndex size_t - The index of the element
159 : *
160 : * @return Vector<size_t, NumElementDOFs> - The global DOF indices
161 : */
162 : KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getGlobalDOFIndices(
163 : const size_t& element_index) const override;
164 :
165 : ///////////////////////////////////////////////////////////////////////
166 : /// Basis functions and gradients /////////////////////////////////////
167 : ///////////////////////////////////////////////////////////////////////
168 :
169 : /**
170 : * @brief Evaluate the shape function of a local degree of freedom at a given point in the
171 : * reference element
172 : *
173 : * @param localDOF size_t - The local degree of freedom index
174 : * @param localPoint point_t (Vector<T, Dim>) - The point in the reference element
175 : *
176 : * @return T - The value of the shape function at the given point
177 : */
178 : KOKKOS_FUNCTION T evaluateRefElementShapeFunction(const size_t& localDOF,
179 : const point_t& localPoint) const override;
180 :
181 : /**
182 : * @brief Evaluate the gradient of the shape function of a local degree of freedom at a
183 : * given point in the reference element
184 : *
185 : * @param localDOF size_t - The local degree of freedom index
186 : * @param localPoint point_t (Vector<T, Dim>) - The point in the reference element
187 : *
188 : * @return point_t (Vector<T, Dim>) - The gradient of the shape function at the given
189 : * point
190 : */
191 : KOKKOS_FUNCTION point_t evaluateRefElementShapeFunctionGradient(
192 : const size_t& localDOF, const point_t& localPoint) const override;
193 :
194 : ///////////////////////////////////////////////////////////////////////
195 : /// Assembly operations ///////////////////////////////////////////////
196 : ///////////////////////////////////////////////////////////////////////
197 :
198 : /**
199 : * @brief Assemble the left stiffness matrix A of the system Ax = b
200 : *
201 : * @param field The field to assemble the matrix for
202 : * @param evalFunction The lambda telling us the form which A takes
203 : *
204 : * @return FieldLHS - The LHS field containing A*x
205 : */
206 : template <typename F>
207 : FieldLHS evaluateAx(FieldLHS& field, F& evalFunction) const;
208 :
209 : /**
210 : * @brief Assemble the left stiffness matrix A of the system
211 : * but only for the boundary values, so that they can be
212 : * subtracted from the RHS for treatment of Dirichlet BCs
213 : *
214 : * @param field The field to assemble the matrix for
215 : * @param evalFunction The lambda telling us the form which A takes
216 : *
217 : * @return FieldLHS - The LHS field containing A*x
218 : */
219 : template <typename F>
220 : FieldLHS evaluateAx_lift(FieldLHS& field, F& evalFunction) const;
221 :
222 : /**
223 : * @brief Assemble the load vector b of the system Ax = b
224 : *
225 : * @param rhs_field The field to set with the load vector
226 : *
227 : * @return FieldRHS - The RHS field containing b
228 : */
229 : void evaluateLoadVector(FieldRHS& field) const override;
230 :
231 : ///////////////////////////////////////////////////////////////////////
232 : /// Error norm computations ///////////////////////////////////////////
233 : ///////////////////////////////////////////////////////////////////////
234 :
235 : /**
236 : * @brief Given two fields, compute the L2 norm error
237 : *
238 : * @param u_h The numerical solution found using FEM
239 : * @param u_sol The analytical solution (functor)
240 : *
241 : * @return error - The error ||u_h - u_sol||_L2
242 : */
243 : template <typename F>
244 : T computeErrorL2(const FieldLHS& u_h, const F& u_sol) const;
245 :
246 : /**
247 : * @brief Given a field, compute the average
248 : *
249 : * @param u_h The numerical solution found using FEM
250 : *
251 : * @return avg The average of the field on the domain
252 : */
253 : T computeAvg(const FieldLHS& u_h) const;
254 :
255 : private:
256 : /**
257 : * @brief Check if a DOF is on the boundary of the mesh
258 : *
259 : * @param ndindex The NDIndex of the global DOF
260 : *
261 : * @return true - If the DOF is on the boundary
262 : * @return false - If the DOF is not on the boundary
263 : */
264 216 : KOKKOS_FUNCTION bool isDOFOnBoundary(const indices_t& ndindex) const {
265 388 : for (size_t d = 0; d < Dim; ++d) {
266 216 : if (ndindex[d] <= 0 || ndindex[d] >= this->nr_m[d] - 1) {
267 44 : return true;
268 : }
269 : }
270 172 : return false;
271 : }
272 :
273 : Kokkos::View<size_t*> elementIndices;
274 : };
275 :
276 : } // namespace ippl
277 :
278 : #include "FEM/LagrangeSpace.hpp"
279 :
280 : #endif
|
