Skip to content

Pass through failure of compressible material failure. #35

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
hverhelst merged 2 commits into from
Sep 22, 2026
Merged

Pass through failure of compressible material failure. #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0f60e40
Signal compressible plane-stress material failure as data, and fix ps…
hverhelst Sep 7, 2026
de2e45f
copilot suggestion approved
hverhelst Sep 22, 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
44 changes: 42 additions & 2 deletions src/gsMaterialMatrixBaseDim.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1012,7 +1012,27 @@ void gsMaterialMatrixBaseDim<dim,T>::_getMetric(index_t k, T z, const gsMatrix<T
else
ratio = det_def / det_ori;

GISMO_ENSURE(ratio >= 0, "Jacobian determinant is negative! det(Gcov_def) = "<<det_def<<"; det(Gcov_ori) = "<<det_ori);
if (!(ratio >= 0))
{
// An inverted (ratio < 0) or undefined (ratio NaN) element has no admissible squared
// thickness stretch J0^2. The failure is signalled as data - a quiet NaN that propagates
// through the material tensors into the assembled system - and not as an exception: this
// routine is evaluated inside the expression assembler's OpenMP region, where a throw
// cannot be caught by any handler above the parallel construct and terminates the process
// instead. The test is written negated so that a NaN ratio also enters the branch.
bool warn = false;
# pragma omp critical (gsMaterialMatrixBaseDim_getMetric_negativeJacobian_C)
{
static bool warnedNegativeJacobian = false;
warn = !warnedNegativeJacobian;
warnedNegativeJacobian = true;
}
if (warn)
gsWarn<<"_getMetric (Cmat): J0^2 = det(Gcov_def)/det(Gcov_ori) is negative or undefined "
<<"(inverted or degenerate element), m_J0_sq set to NaN (det(Gcov_def)="<<det_def
<<", det(Gcov_ori)="<<det_ori<<")\n";
ratio = std::numeric_limits<T>::quiet_NaN();
}
m_data.mine().m_J0_sq = ratio;
}

Expand All @@ -1034,7 +1054,27 @@ void gsMaterialMatrixBaseDim<dim,T>::_getMetric(const index_t k, const T z) cons
else
ratio = det_def / det_ori;

GISMO_ENSURE(ratio >= 0, "Jacobian determinant is negative! det(Gcov_def) = "<<det_def<<"; det(Gcov_ori) = "<<det_ori<<"\nGcov_def = "<<m_data.mine().m_Gcov_def<<"\n"<<"Acov_def = "<<m_data.mine().m_Acov_def<<"\nBcov_def = "<<m_data.mine().m_Bcov_def);
if (!(ratio >= 0))
{
// An inverted (ratio < 0) or undefined (ratio NaN) element has no admissible squared
// thickness stretch J0^2. The failure is signalled as data - a quiet NaN that propagates
// through the material tensors into the assembled system - and not as an exception: this
// routine is evaluated inside the expression assembler's OpenMP region, where a throw
// cannot be caught by any handler above the parallel construct and terminates the process
// instead. The test is written negated so that a NaN ratio also enters the branch.
bool warn = false;
# pragma omp critical (gsMaterialMatrixBaseDim_getMetric_negativeJacobian)
{
static bool warnedNegativeJacobian = false;
warn = !warnedNegativeJacobian;
warnedNegativeJacobian = true;
}
if (warn)
gsWarn<<"_getMetric: J0^2 = det(Gcov_def)/det(Gcov_ori) is negative or undefined "
<<"(inverted or degenerate element), m_J0_sq set to NaN (det(Gcov_def)="<<det_def
<<", det(Gcov_ori)="<<det_ori<<")\n";
ratio = std::numeric_limits<T>::quiet_NaN();
}
m_data.mine().m_J0_sq = ratio;
}

Expand Down
174 changes: 165 additions & 9 deletions src/gsMaterialMatrixNonlinear.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -591,7 +591,7 @@ gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_pstretch_impl(const
colIdx = j*u.cols()+k;
this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)

C33 = C33s(0,k);
C33 = C33s(0,colIdx);

// Compute c
c.setZero();
Expand Down Expand Up @@ -1872,6 +1872,11 @@ gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compres
gsMatrix<T> result(1, u.cols() * z.rows());
result.setZero();
index_t colIdx;
bool warnedNonFinite = false;
bool warnedDamped = false;
bool warnedExhausted = false;
bool warnedNotConverged = false;
bool warnedNonFiniteGuess = false;
for (index_t k=0; k!=u.cols(); k++)
{
// Evaluate material properties on the quadrature point
Expand Down Expand Up @@ -1899,6 +1904,26 @@ gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compres
c.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
c(2,2) = math::pow(m_data.mine().m_J0_sq,-1.0); // c33
// c(2,2) = 1.0; // c33
// The Newton initial guess c33 = 1/J0^2 is not finite when J0^2 is zero (a
// degenerate element) or NaN (an inverted or undefined one, already flagged
// upstream in gsMaterialMatrixBaseDim::_getMetric). No admissible squared
// thickness stretch exists at such a point and the iteration cannot recover
// from a non-finite start, so the failure is signalled as data - a quiet NaN,
// as on the itmax path below - and not as an exception: this routine is
// evaluated inside the expression assembler's OpenMP region, where a throw
// cannot be caught by any handler above the parallel construct.
if (!math::isfinite(c(2,2)))
{
result(0,colIdx) = std::numeric_limits<T>::quiet_NaN();
if (!warnedNonFiniteGuess)
{
warnedNonFiniteGuess = true;
gsWarn<<"_eval3D_Compressible_C33: non-finite Newton initial guess c33 = 1/J0^2, "
<<"c33 set to NaN without iterating (c(2,2)="<<c(2,2)
<<", m_J0_sq="<<m_data.mine().m_J0_sq<<")\n";
}
continue;
}
cinv.setZero();
cinv.block(0,0,2,2) = m_data.mine().m_Gcon_def.block(0,0,2,2);
cinv(2,2) = 1.0/c(2,2);
Expand All @@ -1911,9 +1936,44 @@ gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compres
dc33 = -2. * S33 / C3333;
for (index_t it = 0; it < itmax; it++)
{
c(2,2) += dc33;
// c(2,2) is a squared through-thickness stretch and is admissible only for
// c(2,2) > 0: the compressible stress evaluates pow(m_J_sq,-1/3) with
// m_J_sq = m_J0_sq*c(2,2), which is NaN for a negative argument, and the NaN
// then propagates through every remaining iterate. The Newton increment is
// therefore backtracked until the trial iterate stays admissible. A non-finite
// increment (C3333 -> 0 at loss of ellipticity, or S33 already NaN) cannot be
// repaired by halving and is left to the convergence check below.
if (math::isfinite(dc33))
{
const T dc33_0 = dc33;
index_t nhalve = 0;
while (!(c(2,2) + dc33 > 0) && nhalve < 50) { dc33 /= 2; ++nhalve; }
if (!(c(2,2) + dc33 > 0))
{
dc33 = dc33_0;
if (!warnedExhausted)
{
warnedExhausted = true;
gsWarn<<"_eval3D_Compressible_C33: backtracking exhausted after 50 halvings, dc33 left unchanged (it="<<it<<", c(2,2)="<<c(2,2)<<", dc33="<<dc33<<", m_J0_sq="<<m_data.mine().m_J0_sq<<")\n";
}
}
else if (nhalve > 0 && !warnedDamped)
{
warnedDamped = true;
gsWarn<<"_eval3D_Compressible_C33: Newton increment dc33 backtracked to keep c(2,2) admissible (it="<<it<<", c(2,2)="<<c(2,2)<<", dc33="<<dc33_0<<" -> "<<dc33<<", m_J0_sq="<<m_data.mine().m_J0_sq<<")\n";
}
}
else
{
if (!warnedNonFinite)
{
warnedNonFinite = true;
gsWarn<<"_eval3D_Compressible_C33: non-finite Newton increment dc33, iterate held (it="<<it<<", c(2,2)="<<c(2,2)<<", dc33="<<dc33<<", m_J0_sq="<<m_data.mine().m_J0_sq<<")\n";
}
dc33 = 0.0;
}

//GISMO_ENSURE(c(2,2)>= 0,"ERROR in iteration "<<it<<"; c(2,2) = " << c(2,2) << " C3333=" << C3333 <<" S33=" << S33<<" dc33 = "<<dc33);
c(2,2) += dc33;
cinv(2,2) = 1.0/c(2,2);

m_data.mine().m_J_sq = m_data.mine().m_J0_sq * c(2,2) ;
Expand All @@ -1922,12 +1982,30 @@ gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compres
C3333 = _Cijkl3D(2,2,2,2,c,cinv); // or _Cijkl???

dc33 = -2. * S33 / C3333;
if (math::lessthan(math::abs(dc33),tol))
if (math::lessthan(math::abs(dc33),tol) && c(2,2) > 0)
{
result(0,colIdx) = c(2,2);
break;
}
GISMO_ENSURE(it != itmax-1,"Error: Method did not converge, S33 = "<<S33<<", dc33 = "<<dc33<<" and tolerance = "<<tol<<"\n");
// Newton on S33(c33)=0 exhausted itmax iterations: no admissible squared
// thickness stretch was found at this quadrature point. The failure is
// signalled as data — a quiet NaN that propagates through the material
// tensors into the assembled system — and not as an exception: this routine
// is evaluated inside the expression assembler's OpenMP region, where a
// throw cannot be caught by any handler above the parallel construct and
// terminates the process instead.
if (it == itmax-1)
{
result(0,colIdx) = std::numeric_limits<T>::quiet_NaN();
if (!warnedNotConverged)
{
warnedNotConverged = true;
gsWarn<<"_eval3D_Compressible_C33: Newton did not converge in "<<itmax
<<" iterations, c33 set to NaN (S33="<<S33<<", dc33="<<dc33
<<", tol="<<tol<<", c(2,2)="<<c(2,2)
<<", m_J0_sq="<<m_data.mine().m_J0_sq<<")\n";
}
}
}
}
}
Expand All @@ -1946,6 +2024,11 @@ gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compres
gsMatrix<T> result(1, u.cols() * z.rows());
result.setZero();
index_t colIdx;
bool warnedNonFinite = false;
bool warnedDamped = false;
bool warnedExhausted = false;
bool warnedNotConverged = false;
bool warnedNonFiniteGuess = false;
for (index_t k=0; k!=u.cols(); k++)
{
// Evaluate material properties on the quadrature point
Expand Down Expand Up @@ -1973,6 +2056,26 @@ gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compres
c.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
c(2,2) = math::pow(m_data.mine().m_J0_sq,-1.0); // c33
// c(2,2) = 1.0; // c33
// The Newton initial guess c33 = 1/J0^2 is not finite when J0^2 is zero (a
// degenerate element) or NaN (an inverted or undefined one, already flagged
// upstream in gsMaterialMatrixBaseDim::_getMetric). No admissible squared
// thickness stretch exists at such a point and the iteration cannot recover
// from a non-finite start, so the failure is signalled as data - a quiet NaN,
// as on the itmax path below - and not as an exception: this routine is
// evaluated inside the expression assembler's OpenMP region, where a throw
// cannot be caught by any handler above the parallel construct.
if (!math::isfinite(c(2,2)))
{
result(0,colIdx) = std::numeric_limits<T>::quiet_NaN();
if (!warnedNonFiniteGuess)
{
warnedNonFiniteGuess = true;
gsWarn<<"_eval3D_Compressible_C33 (Cmat): non-finite Newton initial guess c33 = 1/J0^2, "
<<"c33 set to NaN without iterating (c(2,2)="<<c(2,2)
<<", m_J0_sq="<<m_data.mine().m_J0_sq<<")\n";
}
continue;
}
cinv.setZero();
cinv.block(0,0,2,2) = m_data.mine().m_Gcon_def.block(0,0,2,2);
cinv(2,2) = 1.0/c(2,2);
Expand All @@ -1985,9 +2088,44 @@ gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compres
dc33 = -2. * S33 / C3333;
for (index_t it = 0; it < itmax; it++)
{
c(2,2) += dc33;
// c(2,2) is a squared through-thickness stretch and is admissible only for
// c(2,2) > 0: the compressible stress evaluates pow(m_J_sq,-1/3) with
// m_J_sq = m_J0_sq*c(2,2), which is NaN for a negative argument, and the NaN
// then propagates through every remaining iterate. The Newton increment is
// therefore backtracked until the trial iterate stays admissible. A non-finite
// increment (C3333 -> 0 at loss of ellipticity, or S33 already NaN) cannot be
// repaired by halving and is left to the convergence check below.
if (math::isfinite(dc33))
{
const T dc33_0 = dc33;
index_t nhalve = 0;
while (!(c(2,2) + dc33 > 0) && nhalve < 50) { dc33 /= 2; ++nhalve; }
if (!(c(2,2) + dc33 > 0))
{
dc33 = dc33_0;
if (!warnedExhausted)
{
warnedExhausted = true;
gsWarn<<"_eval3D_Compressible_C33 (Cmat): backtracking exhausted after 50 halvings, dc33 left unchanged (it="<<it<<", c(2,2)="<<c(2,2)<<", dc33="<<dc33<<", m_J0_sq="<<m_data.mine().m_J0_sq<<")\n";
}
}
else if (nhalve > 0 && !warnedDamped)
{
warnedDamped = true;
gsWarn<<"_eval3D_Compressible_C33 (Cmat): Newton increment dc33 backtracked to keep c(2,2) admissible (it="<<it<<", c(2,2)="<<c(2,2)<<", dc33="<<dc33_0<<" -> "<<dc33<<", m_J0_sq="<<m_data.mine().m_J0_sq<<")\n";
}
}
else
{
if (!warnedNonFinite)
{
warnedNonFinite = true;
gsWarn<<"_eval3D_Compressible_C33 (Cmat): non-finite Newton increment dc33, iterate held (it="<<it<<", c(2,2)="<<c(2,2)<<", dc33="<<dc33<<", m_J0_sq="<<m_data.mine().m_J0_sq<<")\n";
}
dc33 = 0.0;
}

//GISMO_ENSURE(c(2,2)>= 0,"ERROR in iteration "<<it<<"; c(2,2) = " << c(2,2) << " C3333=" << C3333 <<" S33=" << S33<<" dc33 = "<<dc33);
c(2,2) += dc33;
cinv(2,2) = 1.0/c(2,2);

m_data.mine().m_J_sq = m_data.mine().m_J0_sq * c(2,2) ;
Expand All @@ -1996,12 +2134,30 @@ gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compres
C3333 = _Cijkl3D(2,2,2,2,c,cinv); // or _Cijkl???

dc33 = -2. * S33 / C3333;
if (math::lessthan(math::abs(dc33),tol))
if (math::lessthan(math::abs(dc33),tol) && c(2,2) > 0)
{
result(0,colIdx) = c(2,2);
break;
}
GISMO_ENSURE(it != itmax-1,"Error: Method did not converge, S33 = "<<S33<<", dc33 = "<<dc33<<" and tolerance = "<<tol<<"\n");
// Newton on S33(c33)=0 exhausted itmax iterations: no admissible squared
// thickness stretch was found at this quadrature point. The failure is
// signalled as data — a quiet NaN that propagates through the material
// tensors into the assembled system — and not as an exception: this routine
// is evaluated inside the expression assembler's OpenMP region, where a
// throw cannot be caught by any handler above the parallel construct and
// terminates the process instead.
if (it == itmax-1)
{
result(0,colIdx) = std::numeric_limits<T>::quiet_NaN();
if (!warnedNotConverged)
{
warnedNotConverged = true;
gsWarn<<"_eval3D_Compressible_C33 (Cmat): Newton did not converge in "<<itmax
<<" iterations, c33 set to NaN (S33="<<S33<<", dc33="<<dc33
<<", tol="<<tol<<", c(2,2)="<<c(2,2)
<<", m_J0_sq="<<m_data.mine().m_J0_sq<<")\n";
}
}
}
}
}
Expand Down
5 changes: 5 additions & 0 deletions src/gsThinShellAssembler.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -533,6 +533,11 @@ class gsThinShellAssembler : public gsThinShellAssemblerBase<T>
bool _isInPlane(const boundaryInterface & ifc, const T tol = 1e-2);

private:
/// True if the assembled matrix's stored values are all finite (no NaN/Inf).
bool _isSystemFinite(const gsSparseMatrix<T> & K) const;
/// True if the assembled vector's entries are all finite (no NaN/Inf).
bool _isSystemFinite(const gsMatrix<T> & v) const;

template<short_t _d>
typename std::enable_if<(_d==3), void>::type
_assembleNeumann_impl();
Expand Down
Loading
Loading