DOLFINx
DOLFINx C++ interface
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 <basix/element-families.h>
11#include <cstdint>
12#include <dolfinx/common/math.h>
13#include <dolfinx/mesh/cell_types.h>
14#include <memory>
15#include <xtensor/xtensor.hpp>
16#include <xtl/xspan.hpp>
17
18namespace basix
19{
20class FiniteElement;
21}
22
23namespace dolfinx::fem
24{
25
31{
32public:
35 explicit CoordinateElement(
36 std::shared_ptr<const basix::FiniteElement> element);
37
44 basix::element::lagrange_variant type
45 = basix::element::lagrange_variant::equispaced);
46
48 virtual ~CoordinateElement() = default;
49
53
55 int degree() const;
56
63 int dim() const;
64
66 basix::element::lagrange_variant variant() const;
67
74 std::array<std::size_t, 4> tabulate_shape(std::size_t nd,
75 std::size_t num_points) const;
76
84 xt::xtensor<double, 4> tabulate(int nd,
85 const xt::xtensor<double, 2>& X) const;
86
95 void tabulate(int nd, const xt::xtensor<double, 2>& X,
96 xt::xtensor<double, 4>& basis) const;
97
106 template <typename U, typename V, typename W>
107 static void compute_jacobian(const U& dphi, const V& cell_geometry, W&& J)
108 {
109 math::dot(cell_geometry, dphi, J, true);
110 }
111
115 template <typename U, typename V>
116 static void compute_jacobian_inverse(const U& J, V&& K)
117 {
118 const int gdim = J.shape(1);
119 const int tdim = K.shape(1);
120 if (gdim == tdim)
121 math::inv(J, K);
122 else
123 math::pinv(J, K);
124 }
125
129 template <typename U>
130 static double compute_jacobian_determinant(const U& J)
131 {
132 if (J.shape(0) == J.shape(1))
133 return math::det(J);
134 else
135 {
136 using T = typename U::value_type;
137 auto B = xt::transpose(J);
138 xt::xtensor<T, 2> BA = xt::zeros<T>({B.shape(0), J.shape(1)});
139 math::dot(B, J, BA);
140 return std::sqrt(math::det(BA));
141 }
142 }
143
146
152 static void push_forward(xt::xtensor<double, 2>& x,
153 const xt::xtensor<double, 2>& cell_geometry,
154 const xt::xtensor<double, 2>& phi);
155
162 static std::array<double, 3> x0(const xt::xtensor<double, 2>& cell_geometry);
163
174 static void pull_back_affine(xt::xtensor<double, 2>& X,
175 const xt::xtensor<double, 2>& K,
176 const std::array<double, 3>& x0,
177 const xt::xtensor<double, 2>& x);
178
190 void pull_back_nonaffine(xt::xtensor<double, 2>& X,
191 const xt::xtensor<double, 2>& x,
192 const xt::xtensor<double, 2>& cell_geometry,
193 double tol = 1.0e-8, int maxit = 10) const;
194
196 void permute_dofs(const xtl::span<std::int32_t>& dofs,
197 std::uint32_t cell_perm) const;
198
200 void unpermute_dofs(const xtl::span<std::int32_t>& dofs,
201 std::uint32_t cell_perm) const;
202
208 bool needs_dof_permutations() const;
209
212 bool is_affine() const noexcept { return _is_affine; }
213
214private:
215 // Flag denoting affine map
216 bool _is_affine;
217
218 // Basix Element
219 std::shared_ptr<const basix::FiniteElement> _element;
220};
221} // namespace dolfinx::fem
Definition: CoordinateElement.h:31
ElementDofLayout create_dof_layout() const
Compute and return the dof layout.
Definition: CoordinateElement.cpp:61
static void compute_jacobian(const U &dphi, const V &cell_geometry, W &&J)
Compute Jacobian for a cell with given geometry using the basis functions and first order derivatives...
Definition: CoordinateElement.h:107
static double compute_jacobian_determinant(const U &J)
Compute the determinant of the Jacobian.
Definition: CoordinateElement.h:130
static void push_forward(xt::xtensor< double, 2 > &x, const xt::xtensor< double, 2 > &cell_geometry, const xt::xtensor< double, 2 > &phi)
Compute physical coordinates x for points X in the reference configuration.
Definition: CoordinateElement.cpp:68
static void pull_back_affine(xt::xtensor< double, 2 > &X, const xt::xtensor< double, 2 > &K, const std::array< double, 3 > &x0, const xt::xtensor< double, 2 > &x)
Compute reference coordinates X for physical coordinates x for an affine map. For the affine case,...
Definition: CoordinateElement.cpp:89
basix::element::lagrange_variant variant() const
The variant of the element.
Definition: CoordinateElement.cpp:213
mesh::CellType cell_shape() const
Cell shape.
Definition: CoordinateElement.cpp:38
static void compute_jacobian_inverse(const U &J, V &&K)
Compute the inverse of the Jacobian.
Definition: CoordinateElement.h:116
void pull_back_nonaffine(xt::xtensor< double, 2 > &X, const xt::xtensor< double, 2 > &x, const xt::xtensor< double, 2 > &cell_geometry, double tol=1.0e-8, int maxit=10) const
Compute reference coordinates X for physical coordinates x for a non-affine map.
Definition: CoordinateElement.cpp:106
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 FiniteElement::tabulate
Definition: CoordinateElement.cpp:44
static std::array< double, 3 > x0(const xt::xtensor< double, 2 > &cell_geometry)
Compute the physical coordinate of the reference point X=(0 , 0, 0)
Definition: CoordinateElement.cpp:81
int degree() const
The polynomial degree of the element.
Definition: CoordinateElement.cpp:201
int dim() const
The dimension of the geometry element space.
Definition: CoordinateElement.cpp:207
void unpermute_dofs(const xtl::span< std::int32_t > &dofs, std::uint32_t cell_perm) const
Reverses a DOF permutation.
Definition: CoordinateElement.cpp:187
virtual ~CoordinateElement()=default
Destructor.
bool is_affine() const noexcept
Check is geometry map is affine.
Definition: CoordinateElement.h:212
CoordinateElement(std::shared_ptr< const basix::FiniteElement > element)
Create a coordinate element from a Basix element.
Definition: CoordinateElement.cpp:19
xt::xtensor< double, 4 > tabulate(int nd, const xt::xtensor< double, 2 > &X) const
Evaluate basis values and derivatives at set of points.
Definition: CoordinateElement.cpp:50
void permute_dofs(const xtl::span< std::int32_t > &dofs, std::uint32_t cell_perm) const
Permutes a list of DOF numbers on a cell.
Definition: CoordinateElement.cpp:180
bool needs_dof_permutations() const
Indicates whether the geometry DOF numbers on each cell need permuting.
Definition: CoordinateElement.cpp:194
The class represents the degree-of-freedom (dofs) for an element. Dofs are associated with a mesh ent...
Definition: ElementDofLayout.h:31
Finite element method functionality.
Definition: assemble_matrix_impl.h:24
CellType
Cell type identifier.
Definition: cell_types.h:22