Skip to content

Commit

Permalink
serendipity 2d but this time it's correct
Browse files Browse the repository at this point in the history
  • Loading branch information
mscroggs committed Sep 7, 2023
1 parent f6481a5 commit 217549c
Showing 1 changed file with 2 additions and 29 deletions.
31 changes: 2 additions & 29 deletions cpp/basix/e-serendipity.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -42,37 +42,10 @@ impl::mdarray_t<T, 2> make_serendipity_space_2d(int degree)
impl::mdarray_t<T, 2> wcoeffs(ndofs, psize);
int row_n = 0;
for (int i = 0; i <= degree; ++i)
for (int j = 0; j <= degree - i; ++j)
for (int j = 0; j <= (i < 2 ? degree : (i == degree ? 1 : degree - i)); ++j)
wcoeffs(row_n++, i * (degree + 1) + j) = 1;

if (degree == 1)
{
for (std::size_t i = 0; i < psize; ++i)
{
wcoeffs(row_n, i) = 0.0;
for (std::size_t j = 0; j < wts.size(); ++j)
wcoeffs(row_n, i) += wts[j] * pts(j, 0) * pts(j, 1) * Pq(0, i, j);
}
return wcoeffs;
}

std::vector<T> integrand(wts.size());
for (std::size_t k = 0; k < psize; ++k)
{
for (std::size_t a = 0; a < 2; ++a)
{
for (std::size_t i = 0; i < integrand.size(); ++i)
integrand[i] = wts[i] * pts(i, 0) * pts(i, 1) * Pq(0, k, i);

for (int i = 1; i < degree; ++i)
for (std::size_t j = 0; j < integrand.size(); ++j)
integrand[j] *= pts(j, a);

wcoeffs(row_n + a, k) = 0;
for (std::size_t i = 0; i < integrand.size(); ++i)
wcoeffs(row_n + a, k) += integrand[i];
}
}
assert(std::size_t(row_n) == ndofs);

return wcoeffs;
}
Expand Down

0 comments on commit 217549c

Please sign in to comment.