Branch data 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
|
