Skip to content

Arbitrary order 1d Lagrange FEs #4360

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

Draft
nmnobre wants to merge 5 commits into
base: devel
Choose a base branch
from
+340 −96
Draft

Arbitrary order 1d Lagrange FEs #4360

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
33329d4
Add support to generate 1d Lagrange shape functions of arbitrary order
nmnobre Dec 30, 2025
609c028
Add unit test for generated 1d Lagrange shape functions
nmnobre Dec 30, 2025
c1fed60
Re-bootstrap
nmnobre Dec 30, 2025
7c76f9c
Generalise EDGE3 for arbitrary order Lagrange
nmnobre Jan 2, 2026
741dbfa
Expand unit test to check for symmetries in the shape functions
nmnobre Jan 6, 2026
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
138 changes: 122 additions & 16 deletions include/fe/fe_lagrange_shape_1D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,13 +94,36 @@ Real fe_lagrange_1D_cubic_shape(const unsigned int i,



inline
Real fe_lagrange_1D_any_shape(const Order order,
const unsigned int i,
const Real xi)
{
libmesh_assert_greater_equal (order, 0);
libmesh_assert_less (i, order + 1);

// So the interval endpoints are at i = 0 and i = 1
const unsigned ii = i == 0 ? 0 : i == 1 ? int(order) : i - 1;

auto root = [&](const unsigned int k) { return 1. * k / order; };
const unsigned int nzeros = order + 1;
Real x = (xi + 1.) / 2.;
Real val = 1.;

for (const auto z : make_range(nzeros))
if (z != ii)
val *= (x - root(z)) / (root(ii) - root(z));

return val;
}



inline
Real fe_lagrange_1D_shape(const Order order,
const unsigned int i,
const Real xi)
{
libmesh_assert_less_equal(order, THIRD);

switch (order)
{
// Lagrange linears
Expand All @@ -112,9 +135,12 @@ Real fe_lagrange_1D_shape(const Order order,
return fe_lagrange_1D_quadratic_shape(i, xi);

// Lagrange cubics
// case THIRD
default:
case THIRD:
return fe_lagrange_1D_cubic_shape(i, xi);

// Lagrange arbitrary shape functions
default:
return fe_lagrange_1D_any_shape(order, i, xi);
}
}

Expand Down Expand Up @@ -196,14 +222,51 @@ Real fe_lagrange_1D_cubic_shape_deriv(const unsigned int i,



inline
Real fe_lagrange_1D_any_shape_deriv(const Order order,
const unsigned int i,
const unsigned int libmesh_dbg_var(j),
const Real xi)
{
// only d()/dxi in 1D!
libmesh_assert_equal_to (j, 0);

libmesh_assert_greater_equal (order, 0);
libmesh_assert_less (i, order + 1);

// So the interval endpoints are at i = 0 and i = 1
const unsigned ii = i == 0 ? 0 : i == 1 ? int(order) : i - 1;

auto root = [&](const unsigned int k) { return 1. * k / order; };
const unsigned int nzeros = order + 1;
Real x = (xi + 1.) / 2.;
Real val = 0.;

for (const auto z : make_range(nzeros))
if (z != ii)
{
Real prod = .5;
for (const auto zz : make_range(nzeros))
if (zz != z && zz != ii)
prod *= (x - root(zz));
val += prod;
}

for (const auto z : make_range(nzeros))
if (z != ii)
val /= (root(ii) - root(z));

return val;
}



inline
Real fe_lagrange_1D_shape_deriv(const Order order,
const unsigned int i,
const unsigned int j,
const Real xi)
{
libmesh_assert_less_equal(order, THIRD);

switch (order)
{
case FIRST:
Expand All @@ -212,9 +275,11 @@ Real fe_lagrange_1D_shape_deriv(const Order order,
case SECOND:
return fe_lagrange_1D_quadratic_shape_deriv(i, j, xi);

// case THIRD
default:
case THIRD:
return fe_lagrange_1D_cubic_shape_deriv(i, j, xi);

default:
return fe_lagrange_1D_any_shape_deriv(order, i, j, xi);
}
}

Expand All @@ -229,9 +294,9 @@ Real fe_lagrange_1D_quadratic_shape_second_deriv(const unsigned int i,
const unsigned int libmesh_dbg_var(j),
const Real)
{
// Don't need to switch on j. 1D shape functions
// depend on xi only!
// Don't need to switch on j. 1D shape functions depend on xi only!
libmesh_assert_equal_to (j, 0);

libmesh_assert_less(i, 3);

switch (i)
Expand All @@ -254,9 +319,9 @@ Real fe_lagrange_1D_cubic_shape_second_deriv(const unsigned int i,
const unsigned int libmesh_dbg_var(j),
const Real xi)
{
// Don't need to switch on j. 1D shape functions
// depend on xi only!
// Don't need to switch on j. 1D shape functions depend on xi only!
libmesh_assert_equal_to (j, 0);

libmesh_assert_less(i, 4);

switch (i)
Expand All @@ -278,14 +343,53 @@ Real fe_lagrange_1D_cubic_shape_second_deriv(const unsigned int i,



inline
Real fe_lagrange_1D_any_shape_second_deriv(const Order order,
const unsigned int i,
const unsigned int libmesh_dbg_var(j),
const Real xi)
{
// Don't need to switch on j. 1D shape functions depend on xi only!
libmesh_assert_equal_to (j, 0);

libmesh_assert_greater_equal (order, 0);
libmesh_assert_less (i, order + 1);

// So the interval endpoints are at i = 0 and i = 1
const unsigned ii = i == 0 ? 0 : i == 1 ? int(order) : i - 1;

auto root = [&](const unsigned int k) { return 1. * k / order; };
const unsigned int nzeros = order + 1;
Real x = (xi + 1.) / 2.;
Real val = 0.;

for (const auto z : make_range(nzeros))
if (z != ii)
for (const auto zz : make_range(nzeros))
if (zz != z && zz != ii)
{
Real prod = .25;
for (const auto zzz : make_range(nzeros))
if (zzz != zz && zzz != z && zzz != ii)
prod *= (x - root(zzz));
val += prod;
}

for (const auto z : make_range(nzeros))
if (z != ii)
val /= (root(ii) - root(z));

return val;
}



inline
Real fe_lagrange_1D_shape_second_deriv(const Order order,
const unsigned int i,
const unsigned int j,
const Real xi)
{
libmesh_assert_less_equal(order, THIRD);

switch (order)
{
// All second derivatives of linears are zero....
Expand All @@ -295,9 +399,11 @@ Real fe_lagrange_1D_shape_second_deriv(const Order order,
case SECOND:
return fe_lagrange_1D_quadratic_shape_second_deriv(i, j, xi);

// case THIRD
default:
case THIRD:
return fe_lagrange_1D_cubic_shape_second_deriv(i, j, xi);

default:
return fe_lagrange_1D_any_shape_second_deriv(order, i, j, xi);
} // end switch (order)
}

Expand Down
3 changes: 2 additions & 1 deletion src/fe/fe_interface.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2149,9 +2149,10 @@ unsigned int FEInterface::max_order(const FEType & fe_t,
switch (el_t)
{
case EDGE2:
case EDGE3:
case EDGE4:
return 3;
case EDGE3:
return unlimited;
case TRI3:
case TRISHELL3:
case C0POLYGON:
Expand Down
Loading