DOLFINx 0.9.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
CoordinateElement.h
1// Copyright (C) 2018-2020 Garth N. Wells and Chris N. Richardson
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include "ElementDofLayout.h"
10#include <algorithm>
11#include <array>
12#include <basix/element-families.h>
13#include <basix/mdspan.hpp>
14#include <cmath>
15#include <concepts>
16#include <cstdint>
17#include <dolfinx/common/math.h>
18#include <dolfinx/mesh/cell_types.h>
19#include <limits>
20#include <memory>
21#include <span>
22
23namespace basix
24{
25template <std::floating_point T>
26class FiniteElement;
27}
28
29namespace dolfinx::fem
30{
36template <std::floating_point T>
37class CoordinateElement
38{
39public:
42 explicit CoordinateElement(
43 std::shared_ptr<const basix::FiniteElement<T>> element);
44
51 basix::element::lagrange_variant type
52 = basix::element::lagrange_variant::unset);
53
55 virtual ~CoordinateElement() = default;
56
60
63 int degree() const;
64
71 int dim() const;
72
74 basix::element::lagrange_variant variant() const;
75
82 std::array<std::size_t, 4> tabulate_shape(std::size_t nd,
83 std::size_t num_points) const;
84
93 void tabulate(int nd, std::span<const T> X, std::array<std::size_t, 2> shape,
94 std::span<T> basis) const;
95
112 void permute_subentity_closure(std::span<std::int32_t> d,
113 std::uint32_t cell_info,
114 mesh::CellType entity_type,
115 int entity_index) const;
116
125 template <typename U, typename V, typename W>
126 static void compute_jacobian(const U& dphi, const V& cell_geometry, W&& J)
127 {
128 math::dot(cell_geometry, dphi, J, true);
129 }
130
134 template <typename U, typename V>
135 static void compute_jacobian_inverse(const U& J, V&& K)
136 {
137 int gdim = J.extent(0);
138 int tdim = K.extent(0);
139 if (gdim == tdim)
140 math::inv(J, K);
141 else
142 math::pinv(J, K);
143 }
144
150 template <typename U>
151 static double
152 compute_jacobian_determinant(const U& J, std::span<typename U::value_type> w)
153 {
154 static_assert(U::rank() == 2, "Must be rank 2");
155 if (J.extent(0) == J.extent(1))
156 return math::det(J);
157 else
158 {
159 assert(w.size() >= 2 * J.extent(0) * J.extent(1));
160 using X = typename U::element_type;
161 using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
162 X, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
163 mdspan2_t B(w.data(), J.extent(1), J.extent(0));
164 mdspan2_t BA(w.data() + J.extent(0) * J.extent(1), B.extent(0),
165 J.extent(1));
166 for (std::size_t i = 0; i < B.extent(0); ++i)
167 for (std::size_t j = 0; j < B.extent(1); ++j)
168 B(i, j) = J(j, i);
169
170 // Zero working memory of BA
171 std::fill_n(BA.data_handle(), BA.size(), 0);
172 math::dot(B, J, BA);
173 return std::sqrt(math::det(BA));
174 }
175 }
176
179
187 template <typename U, typename V, typename W>
188 static void push_forward(U&& x, const V& cell_geometry, const W& phi)
189 {
190 for (std::size_t i = 0; i < x.extent(0); ++i)
191 for (std::size_t j = 0; j < x.extent(1); ++j)
192 x(i, j) = 0;
193
194 // Compute x = phi * cell_geometry;
195 math::dot(phi, cell_geometry, x);
196 }
197
208 template <typename U, typename V, typename W>
209 static void pull_back_affine(U&& X, const V& K, const std::array<T, 3>& x0,
210 const W& x)
211 {
212 assert(X.extent(0) == x.extent(0));
213 assert(X.extent(1) == K.extent(0));
214 assert(x.extent(1) == K.extent(1));
215 for (std::size_t i = 0; i < X.extent(0); ++i)
216 for (std::size_t j = 0; j < X.extent(1); ++j)
217 X(i, j) = 0;
218
219 // Calculate X for each point
220 for (std::size_t p = 0; p < x.extent(0); ++p)
221 for (std::size_t i = 0; i < K.extent(0); ++i)
222 for (std::size_t j = 0; j < K.extent(1); ++j)
223 X(p, i) += K(i, j) * (x(p, j) - x0[j]);
224 }
225
227 template <typename X>
228 using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
229 X, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
230
243 mdspan2_t<const T> cell_geometry,
244 double tol = 1.0e-6, int maxit = 15) const;
245
247 void permute(std::span<std::int32_t> dofs, std::uint32_t cell_perm) const;
248
250 void permute_inv(std::span<std::int32_t> dofs, std::uint32_t cell_perm) const;
251
257 bool needs_dof_permutations() const;
258
261 bool is_affine() const noexcept { return _is_affine; }
262
263private:
264 // Flag denoting affine map
265 bool _is_affine;
266
267 // Basix Element
268 std::shared_ptr<const basix::FiniteElement<T>> _element;
269};
270} // namespace dolfinx::fem
Definition utils.h:42
ElementDofLayout create_dof_layout() const
Compute and return the dof layout.
Definition CoordinateElement.cpp:75
void permute_subentity_closure(std::span< std::int32_t > d, std::uint32_t cell_info, mesh::CellType entity_type, int entity_index) const
Given the closure DOFs of a cell sub-entity in reference ordering, this function computes the permut...
Definition CoordinateElement.cpp:64
static void compute_jacobian(const U &dphi, const V &cell_geometry, W &&J)
Definition CoordinateElement.h:126
static void pull_back_affine(U &&X, const V &K, const std::array< T, 3 > &x0, const W &x)
Compute reference coordinates X for physical coordinates x for an affine map. For the affine case,...
Definition CoordinateElement.h:209
void tabulate(int nd, std::span< const T > X, std::array< std::size_t, 2 > shape, std::span< T > basis) const
Evaluate basis values and derivatives at set of points.
Definition CoordinateElement.cpp:55
CoordinateElement(std::shared_ptr< const basix::FiniteElement< T > > element)
Create a coordinate element from a Basix element.
Definition CoordinateElement.cpp:19
void permute_inv(std::span< std::int32_t > dofs, std::uint32_t cell_perm) const
Reverses a DOF permutation.
Definition CoordinateElement.cpp:183
basix::element::lagrange_variant variant() const
Variant of the element.
Definition CoordinateElement.cpp:213
mesh::CellType cell_shape() const
Cell shape.
Definition CoordinateElement.cpp:41
static void compute_jacobian_inverse(const U &J, V &&K)
Compute the inverse of the Jacobian.
Definition CoordinateElement.h:135
MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< X, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 > > mdspan2_t
mdspan typedef
Definition CoordinateElement.h:228
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Shape of array to fill when calling tabulate.
Definition CoordinateElement.cpp:48
int degree() const
The polynomial degree of the element.
Definition CoordinateElement.cpp:199
void permute(std::span< std::int32_t > dofs, std::uint32_t cell_perm) const
Permute a list of DOF numbers on a cell.
Definition CoordinateElement.cpp:175
int dim() const
The dimension of the coordinate element space.
Definition CoordinateElement.cpp:206
void pull_back_nonaffine(mdspan2_t< T > X, mdspan2_t< const T > x, mdspan2_t< const T > cell_geometry, double tol=1.0e-6, int maxit=15) const
Compute reference coordinates X for physical coordinates x for a non-affine map.
Definition CoordinateElement.cpp:83
virtual ~CoordinateElement()=default
Destructor.
bool is_affine() const noexcept
Definition CoordinateElement.h:261
static void push_forward(U &&x, const V &cell_geometry, const W &phi)
Compute physical coordinates x for points X in the reference configuration.
Definition CoordinateElement.h:188
static double compute_jacobian_determinant(const U &J, std::span< typename U::value_type > w)
Compute the determinant of the Jacobian.
Definition CoordinateElement.h:152
bool needs_dof_permutations() const
Indicates whether the geometry DOF numbers on each cell need permuting.
Definition CoordinateElement.cpp:191
Definition ElementDofLayout.h:30
Finite element method functionality.
Definition assemble_matrix_impl.h:26
CellType
Cell type identifier.
Definition cell_types.h:22