Skip to content

DegenerateMapping exception for bad Jacobians and/or Newton convergence in FEMap #4286

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

Open
wants to merge 12 commits into
base: devel
Choose a base branch
from
Open

DegenerateMapping exception for bad Jacobians and/or Newton convergence in FEMap #4286

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
e6b6de6
Add DegenerateMapping exception class
roystgnr Oct 21, 2025
ef0983a
Add libmesh_degenerate_mapping macros
roystgnr Oct 21, 2025
1a05f44
Use libmesh_degenerate_mapping in FEMap
roystgnr Oct 21, 2025
67fb1ff
Refactor degenerate map checks
roystgnr Oct 21, 2025
5a7c22c
Put elem info in exception, not stderr
roystgnr Oct 21, 2025
cc3899a
Use libmesh_catch in reduced_basis
roystgnr Oct 21, 2025
50aae1e
More expansive DegenerateMapping message
roystgnr Oct 21, 2025
66ff4bc
inverse_map non-convergence isn't a LogicError
roystgnr Oct 21, 2025
e0b97e1
DegenerateMap is a better name
roystgnr Nov 4, 2025
1f8c0fd
Use DegenerateMap for degenerate linear maps
roystgnr Nov 4, 2025
c1ad436
Use the proper Jac for degenerate map output
roystgnr Nov 4, 2025
1f8de15
Newline for cleaner DegenerateMap output
roystgnr Nov 4, 2025
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
12 changes: 12 additions & 0 deletions include/base/libmesh_common.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -398,6 +398,9 @@ struct casting_compare {
//
// The libmesh_convergence_failure() macro
// throws a ConvergenceFailure exception
//
// The libmesh_degenerate_mapping() macro prints a message into and
// throws a DegenerateMap exception
#define libmesh_error_msg(msg) \
do { \
std::stringstream message_stream; \
Expand Down Expand Up @@ -451,6 +454,15 @@ struct casting_compare {
LIBMESH_THROW(libMesh::ConvergenceFailure()); \
} while (0)

#define libmesh_degenerate_mapping_msg(msg) \
do { \
std::stringstream message_stream; \
message_stream << msg << '\n'; \
LIBMESH_THROW(libMesh::DegenerateMap(message_stream.str())); \
} while (0)

#define libmesh_degenerate_mapping(filename) libmesh_degenerate_mapping_msg("")

// The libmesh_example_requires() macro prints a message and calls
// "return 77;" if the condition specified by the macro is not true. This
// macro is used in the example executables, which should run when the
Expand Down
24 changes: 24 additions & 0 deletions include/base/libmesh_exceptions.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,30 @@ class FileError : public std::runtime_error
};


/**
* A class representing the detection of an unexpected degeneracy,
* e.g. a negative-determinant Jacobian in a map expected to be
* positive, or a non-trivial kernel in a map expected to be a
* bijection (such as a singular matrix).
*
* libMesh::FEMap throws this if it encounters a point xi in an
* element's "master space" at which the mapping to physical space has
* a too-small (negative, or zero, or nearly zero) Jacobian
* determinant, where "too-small" is determined by a particular
* library method's assigned tolerance.
*
* libMesh::DenseMatrix throws this if it is asked to solve a system
* with a singular matrix and a method (such as lu_solve()) that
* cannot handle singularities.
*/
class DegenerateMap : public std::runtime_error
{
public:
DegenerateMap(std::string msg="") :
std::runtime_error( "Degenerate map, e.g. negative Jacobian or singular matrix.\n" + msg ) {}
};


/**
* A class representing a solver's failure to converge, to be thrown
* by "libmesh_convergence_failure();" This should be a last
Expand Down
3 changes: 2 additions & 1 deletion include/numerics/dense_matrix_impl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -672,7 +672,8 @@ void DenseMatrix<T>::_lu_decompose ()
}

// If the max abs entry found is zero, the matrix is singular
libmesh_error_msg_if(A(i,i) == libMesh::zero, "Matrix A is singular!");
if (A(i,i) == libMesh::zero)
libmesh_degenerate_mapping_msg ("Matrix A is singular!");

// Scale upper triangle entries of row i by the diagonal entry
// Note: don't scale the diagonal entry itself!
Expand Down
251 changes: 68 additions & 183 deletions src/fe/fe_map.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -464,6 +464,47 @@ void FEMap::compute_single_point_map(const unsigned int dim,
if (calculate_xyz)
libmesh_assert_equal_to(phi_map.size(), elem_nodes.size());

auto check_for_degenerate_map =
[this, elem, p]
(Real det_J) {
if (det_J <= jacobian_tolerance)
{
// Don't call get_info() recursively if we're already
// failing. get_info() calls Elem::volume() which may
// call FE::reinit() and trigger the same failure again.
static bool failing = false;
if (!failing)
{
failing = true;
std::string elem_info = elem->get_info();
failing = false;
if (calculate_xyz)
{
libmesh_degenerate_mapping_msg
("Jacobian " << det_J << " at or under tolerance " <<
jacobian_tolerance << " at point " << xyz[p] <<
" in element:\n" << elem_info);
}
else
{
// In this case xyz[p] is not defined, so don't
// try to print it out.
libmesh_degenerate_mapping_msg
("Jacobian " << det_J << " at or under tolerance " <<
jacobian_tolerance << " at point index " << p <<
" in element:\n" << elem_info);
}
}
else
{
// We were already failing when we called this, so just
// stop the current computation and return with
// incomplete results.
return;
}
}
};

switch (dim)
{
//--------------------------------------------------------------------
Expand Down Expand Up @@ -539,46 +580,7 @@ void FEMap::compute_single_point_map(const unsigned int dim,
{
jac[p] = dxyzdxi_map[p].norm();

if (jac[p] <= jacobian_tolerance)
{
// Don't call print_info() recursively if we're already
// failing. print_info() calls Elem::volume() which may
// call FE::reinit() and trigger the same failure again.
static bool failing = false;
if (!failing)
{
failing = true;
elem->print_info(libMesh::err);
failing = false;
if (calculate_xyz)
{
libmesh_error_msg("ERROR: negative Jacobian " \
<< jac[p] \
<< " at point " \
<< xyz[p] \
<< " in element " \
<< elem->id());
}
else
{
// In this case xyz[p] is not defined, so don't
// try to print it out.
libmesh_error_msg("ERROR: negative Jacobian " \
<< jac[p] \
<< " at point index " \
<< p \
<< " in element " \
<< elem->id());
}
}
else
{
// We were already failing when we called this, so just
// stop the current computation and return with
// incomplete results.
return;
}
}
check_for_degenerate_map(jac[p]);

// The inverse Jacobian entries also come from the
// generalized inverse of T (see also the 2D element
Expand Down Expand Up @@ -787,46 +789,7 @@ void FEMap::compute_single_point_map(const unsigned int dim,
// jac = dx/dxi*dy/deta - dx/deta*dy/dxi
jac[p] = (dx_dxi*dy_deta - dx_deta*dy_dxi);

if (jac[p] <= jacobian_tolerance)
{
// Don't call print_info() recursively if we're already
// failing. print_info() calls Elem::volume() which may
// call FE::reinit() and trigger the same failure again.
static bool failing = false;
if (!failing)
{
failing = true;
elem->print_info(libMesh::err);
failing = false;
if (calculate_xyz)
{
libmesh_error_msg("ERROR: negative Jacobian " \
<< jac[p] \
<< " at point " \
<< xyz[p] \
<< " in element " \
<< elem->id());
}
else
{
// In this case xyz[p] is not defined, so don't
// try to print it out.
libmesh_error_msg("ERROR: negative Jacobian " \
<< jac[p] \
<< " at point index " \
<< p \
<< " in element " \
<< elem->id());
}
}
else
{
// We were already failing when we called this, so just
// stop the current computation and return with
// incomplete results.
return;
}
}
check_for_degenerate_map(jac[p]);

JxW[p] = jac[p]*qw[p];

Expand Down Expand Up @@ -888,47 +851,7 @@ void FEMap::compute_single_point_map(const unsigned int dim,

const Real det = (g11*g22 - g12*g21);

if (det <= jacobian_tolerance)
{
// Don't call print_info() recursively if we're already
// failing. print_info() calls Elem::volume() which may
// call FE::reinit() and trigger the same failure again.
thread_local bool failing = false;
if (!failing)
{
failing = true;
Threads::spin_mutex::scoped_lock lock(_point_inv_err_mutex);
elem->print_info(libMesh::err);
failing = false;
if (calculate_xyz)
{
libmesh_error_msg("ERROR: negative Jacobian " \
<< det \
<< " at point " \
<< xyz[p] \
<< " in element " \
<< elem->id());
}
else
{
// In this case xyz[p] is not defined, so don't
// try to print it out.
libmesh_error_msg("ERROR: negative Jacobian " \
<< det \
<< " at point index " \
<< p \
<< " in element " \
<< elem->id());
}
}
else
{
// We were already failing when we called this, so just
// stop the current computation and return with
// incomplete results.
return;
}
}
check_for_degenerate_map(det);

const Real inv_det = 1./det;
jac[p] = std::sqrt(det);
Expand Down Expand Up @@ -1124,46 +1047,7 @@ void FEMap::compute_single_point_map(const unsigned int dim,
dy_dxi*(dz_deta*dx_dzeta - dx_deta*dz_dzeta) +
dz_dxi*(dx_deta*dy_dzeta - dy_deta*dx_dzeta));

if (jac[p] <= jacobian_tolerance)
{
// Don't call print_info() recursively if we're already
// failing. print_info() calls Elem::volume() which may
// call FE::reinit() and trigger the same failure again.
static bool failing = false;
if (!failing)
{
failing = true;
elem->print_info(libMesh::err);
failing = false;
if (calculate_xyz)
{
libmesh_error_msg("ERROR: negative Jacobian " \
<< jac[p] \
<< " at point " \
<< xyz[p] \
<< " in element " \
<< elem->id());
}
else
{
// In this case xyz[p] is not defined, so don't
// try to print it out.
libmesh_error_msg("ERROR: negative Jacobian " \
<< jac[p] \
<< " at point index " \
<< p \
<< " in element " \
<< elem->id());
}
}
else
{
// We were already failing when we called this, so just
// stop the current computation and return with
// incomplete results.
return;
}
}
check_for_degenerate_map(jac[p]);

JxW[p] = jac[p]*qw[p];

Expand Down Expand Up @@ -1730,8 +1614,11 @@ Point FEMap::inverse_map (const unsigned int dim,
// G = [J]^T [J]
const Real G = dxi*dxi;

if (secure)
libmesh_assert_greater (G, 0.);
if (secure && G <= 0)
libmesh_degenerate_mapping_msg
("inverse_map found a singular Jacobian " <<
" at master point " << p << " in element " <<
elem->id());

const Real Ginv = 1./G;

Expand Down Expand Up @@ -1783,8 +1670,11 @@ Point FEMap::inverse_map (const unsigned int dim,

const Real det = (G11*G22 - G12*G21);

if (secure)
libmesh_assert_not_equal_to (det, 0.);
if (secure && det == 0)
libmesh_degenerate_mapping_msg
("inverse_map found a singular Jacobian " <<
" at master point " << p << " in element " <<
elem->id());

const Real inv_det = 1./det;

Expand Down Expand Up @@ -1867,13 +1757,10 @@ Point FEMap::inverse_map (const unsigned int dim,
// case we can just return a far away point.
if (secure)
{
libMesh::err << "ERROR: Newton scheme encountered a singular Jacobian in element: "
<< elem->id()
<< std::endl;

elem->print_info(libMesh::err);

libmesh_error_msg("Exiting...");
libmesh_degenerate_mapping_msg(
"inverse_map found a singular Jacobian" <<
" at master point " << p << " in element " <<
elem->id());
}
else
{
Expand Down Expand Up @@ -1957,17 +1844,15 @@ Point FEMap::inverse_map (const unsigned int dim,

if (cnt > 2*max_cnt)
{
libMesh::err << "ERROR: Newton scheme FAILED to converge in "
<< cnt
<< " iterations in element "
<< elem->id()
<< " for physical point = "
<< physical_point
<< std::endl;

elem->print_info(libMesh::err);

libmesh_error_msg("Exiting...");
// Whether or not this is a degenerate map in the
// "Jacobian becomes singular or negative" sense,
// it's at least degenerate in the "straightforward
// Newton is failing to invert it" sense.
libmesh_degenerate_mapping_msg(
"inverse_map Newton FAILED to converge in " <<
cnt << " iterations in element " << elem->id() <<
" for physical point = " << physical_point <<
std::endl << elem->get_info());
}
}
// Return a far off point when secure is false - this
Expand Down
Loading