Skip to content

Test for out-of-mesh-value on bad MeshFunction variable indices #4391

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

Merged
roystgnr merged 8 commits into from
Feb 11, 2026
Merged

Test for out-of-mesh-value on bad MeshFunction variable indices #4391

Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
79cdfad
New unit test from ChatGPT 5.2 Codex
roystgnr Feb 10, 2026
f0d4398
Fix test assertion
roystgnr Feb 10, 2026
6de0698
Add an analogous test for MeshFunction hessian
roystgnr Feb 10, 2026
391d628
Better names for new tests
roystgnr Feb 10, 2026
02ca102
Weaken FP tolerance for new hessian tests
roystgnr Feb 10, 2026
aa04683
Fix new tests for parallel runs
roystgnr Feb 10, 2026
eb8ae09
Fix new test for --enable-complex builds
roystgnr Feb 10, 2026
5244b48
No, really fix new tests in parallel
roystgnr Feb 10, 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
161 changes: 161 additions & 0 deletions tests/mesh/mesh_function.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,25 @@ Number trilinear_function (const Point & p,
}


Number bilinear_function (const Point & p,
const Parameters &,
const std::string &,
const std::string &)
{
return 8*p(0) + 80*p(1);
}


Number biquadratic_function (const Point & p,
const Parameters &,
const std::string &,
const std::string &)
{
return 7.5*p(0)*p(0) + 75*p(1)*p(1) + 0.75*p(1)*p(0);
}



class MeshFunctionTest : public CppUnit::TestCase
{
/**
Expand All @@ -44,6 +63,10 @@ public:

#if LIBMESH_DIM > 1
CPPUNIT_TEST( test_subdomain_id_sets );
CPPUNIT_TEST( test_bad_gradient_var_with_out_of_mesh_value );
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
CPPUNIT_TEST( test_bad_hessian_var_with_out_of_mesh_value );
#endif
#ifdef LIBMESH_HAVE_PETSC
CPPUNIT_TEST( vectorMeshFunctionLagrange );
CPPUNIT_TEST( vectorMeshFunctionNedelec );
Expand Down Expand Up @@ -154,6 +177,144 @@ public:
}
}

void test_bad_gradient_var_with_out_of_mesh_value()
{
LOG_UNIT_TEST;

ReplicatedMesh mesh(*TestCommWorld);

MeshTools::Generation::build_square(mesh,
4, 4,
0., 1.,
0., 1.,
QUAD4);

EquationSystems es(mesh);
System & sys = es.add_system<System>("SimpleSystem");
unsigned int u_var = sys.add_variable("u", FIRST, LAGRANGE);

es.init();
sys.project_solution(bilinear_function, nullptr, es.parameters);

std::vector<unsigned int> variables {u_var, libMesh::invalid_uint};
MeshFunction mesh_function(es,
*sys.current_local_solution,
sys.get_dof_map(),
variables);
mesh_function.init();

DenseVector<Number> out_of_mesh_value(2);
out_of_mesh_value(0) = -99.0;
out_of_mesh_value(1) = 5.25;
mesh_function.enable_out_of_mesh_mode(out_of_mesh_value);

std::vector<Gradient> gradients;
const Point p(0.4, 0.6, 0.0);

// On a distributed mesh not every processor may have every
// element
const Elem * elem = (*mesh.sub_point_locator())(p);
bool someone_found_elem = elem;
mesh.comm().max(someone_found_elem);
CPPUNIT_ASSERT(someone_found_elem);

mesh_function.gradient(p, 0.0, gradients);

// Let's only test our evaluation where we know we can evaluate, in parallel
if (!elem || elem->processor_id() != mesh.processor_id())
return;

CPPUNIT_ASSERT_EQUAL(std::size_t(2), gradients.size());

LIBMESH_ASSERT_NUMBERS_EQUAL(8.0, gradients[0](0), TOLERANCE * TOLERANCE);
LIBMESH_ASSERT_NUMBERS_EQUAL(80.0, gradients[0](1), TOLERANCE * TOLERANCE);
if (LIBMESH_DIM > 2)
LIBMESH_ASSERT_NUMBERS_EQUAL(0.0, gradients[0](2), TOLERANCE * TOLERANCE);

LIBMESH_ASSERT_NUMBERS_EQUAL(out_of_mesh_value(1),
gradients[1](0),
TOLERANCE * TOLERANCE);
for (unsigned int d = 1; d < LIBMESH_DIM; ++d)
LIBMESH_ASSERT_NUMBERS_EQUAL(0, gradients[1](d),
TOLERANCE * TOLERANCE);
}

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
void test_bad_hessian_var_with_out_of_mesh_value()
{
LOG_UNIT_TEST;

// Hessian calculations are a little weaker in FP on some systems
const Real tol = TOLERANCE * std::sqrt(TOLERANCE);

ReplicatedMesh mesh(*TestCommWorld);

MeshTools::Generation::build_square(mesh,
4, 4,
0., 1.,
0., 1.,
QUAD9);

EquationSystems es(mesh);
System & sys = es.add_system<System>("SimpleSystem");
unsigned int u_var = sys.add_variable("u", SECOND, LAGRANGE);

es.init();
sys.project_solution(biquadratic_function, nullptr, es.parameters);

std::vector<unsigned int> variables {u_var, libMesh::invalid_uint};
MeshFunction mesh_function(es,
*sys.current_local_solution,
sys.get_dof_map(),
variables);
mesh_function.init();

DenseVector<Number> out_of_mesh_value(2);
out_of_mesh_value(0) = -99.0;
out_of_mesh_value(1) = 5.25;
mesh_function.enable_out_of_mesh_mode(out_of_mesh_value);

std::vector<Tensor> hessians;
const Point p(0.4, 0.6, 0.0);

// On a distributed mesh not every processor may have every
// element
const Elem * elem = (*mesh.sub_point_locator())(p);
bool someone_found_elem = elem;
mesh.comm().max(someone_found_elem);
CPPUNIT_ASSERT(someone_found_elem);

mesh_function.hessian(p, 0.0, hessians);

// Let's only test our evaluation where we know we can evaluate, in parallel
if (!elem || elem->processor_id() != mesh.processor_id())
return;

CPPUNIT_ASSERT_EQUAL(std::size_t(2), hessians.size());

LIBMESH_ASSERT_NUMBERS_EQUAL(15.0, hessians[0](0,0), tol);
LIBMESH_ASSERT_NUMBERS_EQUAL(0.75, hessians[0](0,1), tol);
LIBMESH_ASSERT_NUMBERS_EQUAL(0.75, hessians[0](1,0), tol);
LIBMESH_ASSERT_NUMBERS_EQUAL(150.0, hessians[0](1,1), tol);

if (LIBMESH_DIM > 2)
for (unsigned int d = 0; d < LIBMESH_DIM; ++d)
for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2)
if (d>1 || d2>1)
LIBMESH_ASSERT_NUMBERS_EQUAL(0, hessians[1](d,d2),
tol);

LIBMESH_ASSERT_NUMBERS_EQUAL(out_of_mesh_value(1),
hessians[1](0,0),
tol);
for (unsigned int d = 0; d < LIBMESH_DIM; ++d)
for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2)
if (d || d2)
LIBMESH_ASSERT_NUMBERS_EQUAL(0, hessians[1](d,d2),
tol);
}
#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES

// test that mesh function works correctly with non-zero
// Elem::p_level() values.
#ifdef LIBMESH_ENABLE_AMR
Expand Down