Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix tolerances for determining DOF transformation type and add test #875

Draft
wants to merge 8 commits into
base: main
Choose a base branch
from
+173 −122
Draft

Fix tolerances for determining DOF transformation type and add test #875

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
6bb7f06
Fix tolerances and add test
garth-wells Nov 15, 2024
1868dd3
Make eps consistent
garth-wells Nov 15, 2024
eec2489
Major speed improvement
garth-wells Nov 15, 2024
ff1dfaa
Tidy up
garth-wells Nov 15, 2024
4de32c0
Tidy up
garth-wells Nov 16, 2024
9daa18b
Remove debug include
garth-wells Nov 16, 2024
669c6a0
Small update
garth-wells Nov 16, 2024
1a39ff4
Move include
garth-wells Nov 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
146 changes: 80 additions & 66 deletions cpp/basix/dof-transformations.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@
#include <span>
#include <tuple>

#include <iostream>
garth-wells marked this conversation as resolved.
Show resolved Hide resolved

using namespace basix;

namespace stdex
Expand All @@ -39,8 +41,7 @@ int find_first_subentity(cell::type cell_type, cell::type entity_type)
{
const int edim = cell::topological_dimension(entity_type);
std::vector<cell::type> entities = cell::subentity_types(cell_type)[edim];
if (auto it = std::ranges::find(entities, entity_type);
it != entities.end())
if (auto it = std::ranges::find(entities, entity_type); it != entities.end())
{
return std::distance(entities.begin(), it);
}
Expand Down Expand Up @@ -93,9 +94,7 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
{
mapinfo_t<T> mapinfo;
auto& data = mapinfo.try_emplace(cell::type::interval).first->second;
auto map = [](auto pt) -> std::array<T, 3> {
return {pt[1], pt[0], 0.0};
};
auto map = [](auto pt) -> std::array<T, 3> { return {pt[1], pt[0], 0.0}; };
stdex::mdarray<T, stdex::extents<std::size_t, 2, 2>> J(
stdex::extents<std::size_t, 2, 2>{}, {0., 1., 1., 0.});

Expand All @@ -109,9 +108,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
{
mapinfo_t<T> mapinfo;
auto& data = mapinfo.try_emplace(cell::type::interval).first->second;
auto map = [](auto pt) -> std::array<T, 3> {
return {1 - pt[0], pt[1], 0};
};
auto map
= [](auto pt) -> std::array<T, 3> { return {1 - pt[0], pt[1], 0}; };
stdex::mdarray<T, stdex::extents<std::size_t, 2, 2>> J(
stdex::extents<std::size_t, 2, 2>{}, {-1., 0., 0., 1.});

Expand All @@ -126,9 +124,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
mapinfo_t<T> mapinfo;
{
auto& data = mapinfo.try_emplace(cell::type::interval).first->second;
auto map = [](auto pt) -> std::array<T, 3> {
return {pt[0], pt[2], pt[1]};
};
auto map
= [](auto pt) -> std::array<T, 3> { return {pt[0], pt[2], pt[1]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{1., 0., 0., 0., 0., 1., 0., 1., 0.});
Expand All @@ -143,9 +140,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
{
auto& data = mapinfo.try_emplace(cell::type::triangle).first->second;
{
auto map = [](auto pt) -> std::array<T, 3> {
return {pt[2], pt[0], pt[1]};
};
auto map
= [](auto pt) -> std::array<T, 3> { return {pt[2], pt[0], pt[1]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{0., 1., 0., 0., 0., 1., 1., 0., 0.});
Expand All @@ -157,9 +153,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
data.push_back(std::tuple(map, J, detJ, K));
}
{
auto map = [](auto pt) -> std::array<T, 3> {
return {pt[0], pt[2], pt[1]};
};
auto map
= [](auto pt) -> std::array<T, 3> { return {pt[0], pt[2], pt[1]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{1., 0., 0., 0., 0., 1., 0., 1., 0.});
Expand All @@ -179,9 +174,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
mapinfo_t<T> mapinfo;
{
auto& data = mapinfo.try_emplace(cell::type::interval).first->second;
auto map = [](auto pt) -> std::array<T, 3> {
return {1 - pt[0], pt[1], pt[2]};
};
auto map = [](auto pt) -> std::array<T, 3>
{ return {1 - pt[0], pt[1], pt[2]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{-1., 0., 0., 0., 1., 0., 0., 0., 1.});
Expand All @@ -196,9 +190,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
{
auto& data = mapinfo.try_emplace(cell::type::quadrilateral).first->second;
{
auto map = [](auto pt) -> std::array<T, 3> {
return {1 - pt[1], pt[0], pt[2]};
};
auto map = [](auto pt) -> std::array<T, 3>
{ return {1 - pt[1], pt[0], pt[2]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{0., 1., 0., -1., 0., 0., 0., 0., 1.});
Expand All @@ -210,9 +203,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
data.push_back(std::tuple(map, J, detJ, K));
}
{
auto map = [](auto pt) -> std::array<T, 3> {
return {pt[1], pt[0], pt[2]};
};
auto map
= [](auto pt) -> std::array<T, 3> { return {pt[1], pt[0], pt[2]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{0., 1., 0., 1., 0., 0., 0., 0., 1.});
Expand All @@ -231,9 +223,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
mapinfo_t<T> mapinfo;
{
auto& data = mapinfo.try_emplace(cell::type::interval).first->second;
auto map = [](auto pt) -> std::array<T, 3> {
return {1 - pt[0], pt[1], pt[2]};
};
auto map = [](auto pt) -> std::array<T, 3>
{ return {1 - pt[0], pt[1], pt[2]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{-1., 0., 0., 0., 1., 0., 0., 0., 1.});
Expand All @@ -246,9 +237,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
{
auto& data = mapinfo.try_emplace(cell::type::triangle).first->second;
{
auto map = [](auto pt) -> std::array<T, 3> {
return {1 - pt[1] - pt[0], pt[0], pt[2]};
};
auto map = [](auto pt) -> std::array<T, 3>
{ return {1 - pt[1] - pt[0], pt[0], pt[2]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{0., 1., 0., -1., -1., 0., 0., 0., 1.});
Expand All @@ -259,9 +249,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
data.push_back(std::tuple(map, J, detJ, K));
}
{
auto map = [](auto pt) -> std::array<T, 3> {
return {pt[1], pt[0], pt[2]};
};
auto map
= [](auto pt) -> std::array<T, 3> { return {pt[1], pt[0], pt[2]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{0., 1., 0., 1., 0., 0., 0., 0., 1.});
Expand All @@ -275,9 +264,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
{
auto& data = mapinfo.try_emplace(cell::type::quadrilateral).first->second;
{
auto map = [](auto pt) -> std::array<T, 3> {
return {1 - pt[2], pt[1], pt[0]};
};
auto map = [](auto pt) -> std::array<T, 3>
{ return {1 - pt[2], pt[1], pt[0]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{0., 0., 1., 0., 1., 0., -1., 0., 0.});
Expand All @@ -288,9 +276,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
data.push_back(std::tuple(map, J, detJ, K));
}
{ // scope
auto map = [](auto pt) -> std::array<T, 3> {
return {pt[2], pt[1], pt[0]};
};
auto map
= [](auto pt) -> std::array<T, 3> { return {pt[2], pt[1], pt[0]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{0., 0., 1., 0., 1., 0., 1., 0., 0.});
Expand All @@ -309,9 +296,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
mapinfo_t<T> mapinfo;
{
auto& data = mapinfo.try_emplace(cell::type::interval).first->second;
auto map = [](auto pt) -> std::array<T, 3> {
return {1 - pt[0], pt[1], pt[2]};
};
auto map = [](auto pt) -> std::array<T, 3>
{ return {1 - pt[0], pt[1], pt[2]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{-1., 0., 0., 0., 1., 0., 0., 0., 1.});
Expand All @@ -325,9 +311,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
{
auto& data = mapinfo.try_emplace(cell::type::quadrilateral).first->second;
{
auto map = [](auto pt) -> std::array<T, 3> {
return {1 - pt[1], pt[0], pt[2]};
};
auto map = [](auto pt) -> std::array<T, 3>
{ return {1 - pt[1], pt[0], pt[2]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{0., 1., 0., -1., 0., 0., 0., 0., 1.});
Expand All @@ -338,9 +323,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
data.push_back(std::tuple(map, J, detJ, K));
}
{
auto map = [](auto pt) -> std::array<T, 3> {
return {pt[1], pt[0], pt[2]};
};
auto map
= [](auto pt) -> std::array<T, 3> { return {pt[1], pt[0], pt[2]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{0., 1., 0., 1., 0., 0., 0., 0., 1.});
Expand All @@ -355,9 +339,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
{
auto& data = mapinfo.try_emplace(cell::type::triangle).first->second;
{
auto map = [](auto pt) -> std::array<T, 3> {
return {1 - pt[2] - pt[0], pt[1], pt[0]};
};
auto map = [](auto pt) -> std::array<T, 3>
{ return {1 - pt[2] - pt[0], pt[1], pt[0]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{0., 0., 1., 0., 1., 0., -1., 0., -1.});
Expand All @@ -368,9 +351,8 @@ mapinfo_t<T> get_mapinfo(cell::type cell_type)
data.push_back(std::tuple(map, J, detJ, K));
}
{
auto map = [](auto pt) -> std::array<T, 3> {
return {pt[2], pt[1], pt[0]};
};
auto map
= [](auto pt) -> std::array<T, 3> { return {pt[2], pt[1], pt[0]}; };
stdex::mdarray<T, stdex::extents<std::size_t, 3, 3>> J(
stdex::extents<std::size_t, 3, 3>{},
{0., 0., 1., 0., 1., 0., 1., 0., 0.});
Expand Down Expand Up @@ -426,7 +408,7 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_transformation(
mdarray_t<T, 2> mapped_pts(pts.extents());
for (std::size_t p = 0; p < mapped_pts.extent(0); ++p)
{
auto mp = map_point(
std::array<T, 3> mp = map_point(
std::span(pts.data_handle() + p * pts.extent(1), pts.extent(1)));
for (std::size_t k = 0; k < mapped_pts.extent(1); ++k)
mapped_pts(p, k) = mp[k];
Expand All @@ -439,18 +421,44 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_transformation(
mdspan_t<const T, 2> polyset_vals(polyset_vals_b.data(), polyset_shape[1],
polyset_shape[2]);

mdarray_t<T, 3> tabulated_data(npts, total_ndofs, vs);
std::vector<T> tabulated_data_b(npts * total_ndofs * vs);
mdspan_t<T, 3> tabulated_data(tabulated_data_b.data(), npts, total_ndofs, vs);
std::vector<T> result_b(polyset_vals.extent(1) * coeffs.extent(0));
mdspan_t<T, 2> result(result_b.data(), coeffs.extent(0),
polyset_vals.extent(1));
std::vector<T> coeffs_b(coeffs.extent(0) * polyset_vals.extent(0));
mdspan_t<T, 2> _coeffs(coeffs_b.data(), coeffs.extent(0),
polyset_vals.extent(0));
for (std::size_t j = 0; j < vs; ++j)
{
mdarray_t<T, 2> result(polyset_vals.extent(1), coeffs.extent(0));
// std::fill(result_b.begin(), result_b.end(), 0);

for (std::size_t k0 = 0; k0 < coeffs.extent(0); ++k0)
for (std::size_t k1 = 0; k1 < polyset_vals.extent(1); ++k1)
for (std::size_t k2 = 0; k2 < polyset_vals.extent(0); ++k2)
result(k1, k0) += coeffs(k0, k2 + psize * j) * polyset_vals(k2, k1);
for (std::size_t k2 = 0; k2 < polyset_vals.extent(0); ++k2)
_coeffs(k0, k2) = coeffs(k0, k2 + psize * j);

// R^T = coeffs * polyset_vals
// for (std::size_t k0 = 0; k0 < coeffs.extent(0); ++k0) // big
// for (std::size_t k1 = 0; k1 < polyset_vals.extent(1); ++k1)
// for (std::size_t k2 = 0; k2 < polyset_vals.extent(0); ++k2) // big
// result(k1, k0) += _coeffs(k0, k2) * polyset_vals(k2, k1);
// result(k1, k0) += coeffs(k0, k2 + psize * j) * polyset_vals(k2, k1);

// r^t: coeffs.extent(0) x polyset_vals.extent(1) [k0, k1]
// c: coeffs.extent(1) x polyset_vals.extent(0) [k0, k2]
// p: polyset_vals.extent(0) x polyset_vals.extent(1) [k2, k1]

math::dot(_coeffs, polyset_vals, result);

for (std::size_t k0 = 0; k0 < result.extent(0); ++k0)
for (std::size_t k1 = 0; k1 < result.extent(1); ++k1)
tabulated_data(k0, k1, j) = result(k0, k1);
// std::cout << "4: transform" << std::endl;

for (std::size_t k0 = 0; k0 < result.extent(1); ++k0)
for (std::size_t k1 = 0; k1 < result.extent(0); ++k1)
tabulated_data(k0, k1, j) = result(k1, k0);

// for (std::size_t k0 = 0; k0 < result.extent(0); ++k0)
// for (std::size_t k1 = 0; k1 < result.extent(1); ++k1)
// tabulated_data(k0, k1, j) = result(k0, k1);
}

// push forward
Expand All @@ -460,7 +468,7 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_transformation(
for (std::size_t i = 0; i < npts; ++i)
{
mdspan_t<const T, 2> tab(
tabulated_data.data()
tabulated_data_b.data()
+ i * tabulated_data.extent(1) * tabulated_data.extent(2),
tabulated_data.extent(1), tabulated_data.extent(2));

Expand All @@ -482,10 +490,16 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_transformation(
for (std::size_t i = 0; i < vs; ++i)
{
for (std::size_t k0 = 0; k0 < transform.extent(1); ++k0)
{
for (std::size_t k1 = 0; k1 < transform.extent(0); ++k1)
{
for (std::size_t k2 = 0; k2 < imat.extent(2); ++k2)
{
transform(k1, k0)
+= imat(k0, i, k2, d) * pushed_data(k2, k1 + dofstart, i);
}
}
}
}
}

Expand Down
5 changes: 3 additions & 2 deletions cpp/basix/finite-element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
#include <concepts>
#include <limits>
#include <numeric>

#define str_macro(X) #X
#define str(X) str_macro(X)

Expand Down Expand Up @@ -965,7 +966,7 @@ FiniteElement<F>::FiniteElement(
rtot += r;
}

constexpr F eps = 10.0 * std::numeric_limits<float>::epsilon();
const F eps = 10 * degree * degree * std::numeric_limits<F>::epsilon();
if ((trans.extent(2) != 1 and std::abs(rmin) > eps)
or std::abs(rmax - 1.0) > eps or std::abs(rtot - 1.0) > eps)
{
Expand Down Expand Up @@ -1436,7 +1437,7 @@ FiniteElement<F>::FiniteElement(
for (std::size_t col = 0; col < matM.extent(1); ++col)
{
F v = col == row ? 1.0 : 0.0;
constexpr F eps = 100 * std::numeric_limits<F>::epsilon();
const F eps = 10 * degree * degree * std::numeric_limits<F>::epsilon();
if (std::abs(matM(row, col) - v) > eps)
{
_interpolation_is_identity = false;
Expand Down
1 change: 0 additions & 1 deletion cpp/basix/finite-element.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@
/// Basix: FEniCS runtime basis evaluation library
namespace basix
{

namespace impl
{
template <typename T, std::size_t d>
Expand Down
6 changes: 6 additions & 0 deletions cpp/basix/math.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -330,6 +330,12 @@ void dot(const U& A, const V& B, W&& C)
}
else
{
static_assert(std::is_same_v<typename U::layout_type, std::layout_right>);
static_assert(std::is_same_v<typename V::layout_type, std::layout_right>);
// static_assert(
// std::is_same_v<typename std::remove_cv<typename W::layout_type>,
// std::layout_right>);

using T = typename std::decay_t<U>::value_type;
impl::dot_blas<T>(
std::span(A.data_handle(), A.size()), {A.extent(0), A.extent(1)},
Expand Down
Loading