%~Mouliné par MaN_auto v.0.40.4 (9e9f6e09) 2026-06-23 15:45:30
\documentclass[CRMECA,Unicode,XML,screen,biblatex,published]{cedram}

\COI{The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.}
\addbibresource{CRMECA_MOUMNI_20250960.bib}




\newcommand{\jump}[1]{[\![#1]\!]}
\newcommand{\SBV}{\mathbf{SBV}}
\newcommand{\BV}{\mathbf{BV}}
\newcommand{\Haus}{\mathcal{H}}
\newcommand{\clL}{\mathcal{L}}

\newcommand{\bfx}{\mathbf{x}}
\newcommand{\bfn}{\mathbf{n}}
\newcommand{\bfK}{\mathbf{K}}
\newcommand{\bfU}{\mathbf{U}}
\newcommand{\bfJ}{\mathbf{J}}
\newcommand{\bfF}{\mathbf{F}}
\newcommand{\bfZO}{\mathbf{0}}


\newcommand{\bbR}{\mathbb{R}}

\DeclareMathOperator{\Span}{span}
%\rmd
%\newcommand{\rmd}{\,\mathrm{d}}
\newcommand{\rmd}{\,d}
%A proposer aux auteurs.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Better widebar

\makeatletter
\let\save@mathaccent\mathaccent
\newcommand*\if@single[3]{%
  \setbox0\hbox{${\mathaccent"0362{#1}}^H$}%
  \setbox2\hbox{${\mathaccent"0362{\kern0pt#1}}^H$}%
  \ifdim\ht0=\ht2 #3\else #2\fi
  }
%The bar will be moved to the right by a half of \macc@kerna, which is computed by amsmath:
\newcommand*\rel@kern[1]{\kern#1\dimexpr\macc@kerna}
%If there's a superscript following the bar, then no negative kern may follow the bar;
%an additional {} makes sure that the superscript is high enough in this case:
\newcommand*\widebar[1]{\@ifnextchar^{{\wide@bar{#1}{0}}}{\wide@bar{#1}{1}}}
%Use a separate algorithm for single symbols:
\newcommand*\wide@bar[2]{\if@single{#1}{\wide@bar@{#1}{#2}{1}}{\wide@bar@{#1}{#2}{2}}}
\newcommand*\wide@bar@[3]{%
  \begingroup
  \def\mathaccent##1##2{%
%Enable nesting of accents:
    \let\mathaccent\save@mathaccent
%If there's more than a single symbol, use the first character instead (see below):
    \if#32 \let\macc@nucleus\first@char \fi
%Determine the italic correction:
    \setbox\z@\hbox{$\macc@style{\macc@nucleus}_{}$}%
    \setbox\tw@\hbox{$\macc@style{\macc@nucleus}{}_{}$}%
    \dimen@\wd\tw@
    \advance\dimen@-\wd\z@
%Now \dimen@ is the italic correction of the symbol.
    \divide\dimen@ 3
    \@tempdima\wd\tw@
    \advance\@tempdima-\scriptspace
%Now \@tempdima is the width of the symbol.
    \divide\@tempdima 10
    \advance\dimen@-\@tempdima
%Now \dimen@ = (italic correction / 3) - (Breite / 10)
    \ifdim\dimen@>\z@ \dimen@0pt\fi
%The bar will be shortened in the case \dimen@<0 !
    \rel@kern{0.6}\kern-\dimen@
    \if#31
      \overline{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@}%
      \advance\dimen@0.4\dimexpr\macc@kerna
%Place the combined final kern (-\dimen@) if it is >0 or if a superscript follows:
      \let\final@kern#2%
      \ifdim\dimen@<\z@ \let\final@kern1\fi
      \if\final@kern1 \kern-\dimen@\fi
    \else
      \overline{\rel@kern{-0.6}\kern\dimen@#1}%
    \fi
  }%
  \macc@depth\@ne
  \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar
  \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax}%
  \macc@set@skewchar\relax
  \let\mathaccentV\macc@nested@a
%The following initialises \macc@kerna and calls \mathaccent:
  \if#31
    \macc@nested@a\relax111{#1}%
  \else
%If the argument consists of more than one symbol, and if the first token is
%a letter, use that letter for the computations:
    \def\gobble@till@marker##1\endmarker{}%
    \futurelet\first@char\gobble@till@marker#1\endmarker
    \ifcat\noexpand\first@char A\else
      \def\first@char{}%
    \fi
    \macc@nested@a\relax111{\first@char}%
  \fi
  \endgroup
}
\makeatother

\let\oldbar\bar
\renewcommand*{\bar}[1]{{\mathchoice{\widebar{#1}}{\widebar{#1}}{\widebar{#1}}{\oldbar{#1}}}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}
\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}
\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{A pseudo-continuous media approach}
\alttitle{Mécanique des milieux pseudo-continus}

\keywords{Continuous media, discontinuities, distributions, variational formulation}
\altkeywords{Milieux continus, discontinuités, distributions, formulation variationnelle}

\author{\firstname{Hua} \lastname{Ding} \CDRorcid{0009-0005-5107-3773}\IsCorresp}
\address{Shenzhen TenFeng Technology Co., Xueyuan Avenue, Nanshan District, Shenzhen, China}
\email[H.~Ding]{dingh@tenfong.cn}

\author{\firstname{Ziad} \lastname{Moumni}\CDRorcid{0000-0002-2644-1133}\IsCorresp}
\address{LMI, UME, ENSTA, Institue Polytechnique de Paris, 828 Boulevard des maréchaux, 91120 Palaiseau, Cedex, France}
\email[Z.~Moumni]{ziad.moumni@ensta.fr}

\begin{abstract}
This paper introduces a novel mathematical framework, termed the pseudo-continuous media approach, for treating strong discontinuities (such as cracks, interfaces, and contact surfaces) within continuum mechanics. Traditional methods typically require explicit geometric representation of discontinuities,
leading to complex meshing and algorithmic challenges. In contrast, our method leverages distribution theory to embed the discontinuity conditions directly into the governing equations, allowing the problem to be
posed on a simpler, uncracked domain. Our distributional approach is conceptually aligned with the rigorous framework of {Special Functions of Bounded Variation} (SBV) which provides a natural setting for functions with jump discontinuities. We present the theoretical foundation and the corresponding finite element
formulation. The method demonstrates significant advantages in geometric simplification, mathematical
elegance, and numerical robustness compared to traditional techniques like the {eXtended Finite Element Method (XFEM)}, positioning itself as a novel numerical methodology that implements the SBV philosophy
using a standard finite element framework.
\end{abstract}

\begin{altabstract}
Cet article présente un nouveau cadre mathématique, dénommé mécanique des milieux pseudo-continus, destiné à traiter les discontinuités fortes -- telles que les fissures -- au sein de la mécanique des milieux continus. Les méthodes classiques nécessitent généralement une représentation géométrique explicite des discontinuités, ce qui entraîne des difficultés de maillage et des défis algorithmiques complexes. En revanche, notre méthode exploite la théorie des distributions pour intégrer directement les conditions de discontinuité dans les équations gouvernantes, ce qui permet de poser le problème sur un domaine plus simple, continu. Nous présentons une formulation variationnelle mixte qui sépare la solution en une partie régulière C\textsuperscript{1} continue et un champ de saut indépendant défini sur la surface de la discontinuité. L'approche distributionnelle est conceptuellement alignée sur le cadre rigoureux des {Fonctions Spéciales à Variation Bornée} (SBV), qui fournit un cadre naturel pour les fonctions présentant des discontinuités. Cet article se concentre sur les fondements théoriques de la méthode, incluant la formulation mixte, l'interprétation du cadre SBV et la discrétisation par éléments finis correspondante. Les exemples numériques et les études de validation ne font pas partie du présent travail ; ils seront présentés dans un prochain article dédié aux aspects d'implémentation et aux problèmes de référence. La méthode présente des avantages significatifs en termes de simplification géométrique par rapport aux techniques traditionnelles.
\end{altabstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\dateposted{2026-07-20}
\begin{document}
\maketitle


\section*{Introduction}

Within the standard framework of continuum mechanics, the treatment of strong discontinuities (e.g., cracks) presents a significant challenge. The conventional approach involves isolating the discontinuous surface $\gamma$. The domain is split into $\Omega = \bar{\Omega}\backslash \gamma$, where $\bar{\Omega}$ is the uncracked domain, and the solution is sought separately in the subdomains. These solutions are then coupled through additional conditions imposed on $\gamma$, such as traction-free boundary conditions for a crack.

This work presents a paradigm shift. Instead of adapting the geometry to the discontinuity, we adapt the mathematical formulation. By employing the theory of distributions, we can represent the discontinuous solution as a generalized function (a distribution) defined over the entire, simpler domain $\bar{\Omega}$. This goal is conceptually aligned with the framework of Special Functions of Bounded Variation (SBV) \cite{Ambrosio2000}, which was developed precisely to provide a rigorous functional-analytic setting for functions with jump discontinuities.

\goodbreak
While the SBV framework provides the theoretical justification, our contribution is a novel numerical implementation for a predetermined discontinuity. We develop a mixed finite element strategy that uses:
\begin{itemize}
\item A regular part approximated by $C^1$-continuous elements on a fixed, simple mesh of $\bar{\Omega}$.
\item An independent jump field defined on the crack surface $\gamma$.
\end{itemize}
This approach avoids the need for complex, conforming discretizations or specialized enrichment functions used in methods like XFEM.



The paper is structured as follows: Section~\ref{sect2} introduces a canonical model problem. Section~\ref{sect3} develops the core theoretical framework of pseudo-continuous representations,
including a detailed explanation of the SBV derivative and regularity considerations.

Section~\ref{sect4} details the mixed finite element discretization, including the derivation of the traction-free condition and the explanation of the variational notation.

Section~\ref{sect5} discusses the advantages, limitations, extensions to interfaces and contact, and future work. Section~\ref{sect6} concludes the paper.


\section{A model problem}\label{sect2}

To illustrate the core ideas without the algebraic burden of elasticity, we use the Poisson equation, a prototypical elliptic problem. Consider a domain $\Omega \subset \bbR^n$ $(n=2,3)$ containing an internal crack, modeled as a hypersurface $\gamma \subset \bbR^n$ of dimension $n-1$. The cracked domain is $\Omega = \bar{\Omega} \setminus \gamma$, and its external boundary is $\Gamma = \partial \bar{\Omega}$.

We consider the following mixed boundary value problem. Find $u(\bfx)$, for $\bfx \in \Omega$, such that:
\begin{equation}\label{eq:1}
\begin{cases}
-\Delta u = f & \text{in } \Omega, \\
u = 0 & \text{on } \Gamma, \\
\frac{\partial u}{\partial n} = 0 & \text{on } \gamma^+ \quad\text{and}\quad \gamma^-.
\end{cases} %\qquad (1)
\end{equation}

Here, the homogeneous Neumann condition on $\gamma$ models a traction-free crack surface in a simplified setting.

It is well known that this problem has a unique solution in the Sobolev space $\widetilde{H_0^1}(\Omega)$~\cite{2}, defined~as:
\[
\widetilde{H_0^1}(\Omega) = \left\{ v \in H^1(\Omega) \,\middle|\, v|_{\Gamma} = 0 \right\},
\]
where $\widetilde{H_0^1}(\Omega)$ is the space of functions which are not required to be continuous across $\gamma$ and not required to vanish on the entire boundary $\partial \Omega$ since they are not required to vanish on $\gamma$, defined~as:
\[
\widetilde{H_0^1}(\Omega) = \left\{ v \in H^1(\Omega) \,\middle|\, v|_{\Gamma} = 0 \right\}.
\]

The solution satisfies the standard variational formulation:
\begin{equation}\label{eq:2}
\int_{\Omega} \nabla u \cdot \nabla v \rmd\bfx = \int_{\Omega} f v \rmd\bfx \quad \forall v \in \widetilde{H_0^1}(\Omega). %\qquad (2)
\end{equation}

A key feature of the solution is its \emph{jump} across the crack:
\[
[u]_\gamma := u|_{\gamma^+} - u|_{\gamma^-} \neq 0.
\]
Our goal is to capture this behavior without explicitly cutting the domain.

\pagebreak
\section{Pseudo-continuous representation using distribution theory}\label{sect3}

\subsection{Extension by distributions}

The objective is to define an extension $\bar{u}$ of the solution $u$ from the cracked domain $\Omega$ to the entire, uncracked domain $\bar{\Omega}$. This extension must encapsulate the discontinuous behavior of $u$ across $\gamma$. Since $\bar{u}$ is discontinuous, its classical derivatives do not exist on $\gamma$; we must therefore define them in the sense of distributions.

\subsection{Deriving the distributional formulation}

Starting from the solution $u \in H_0^1(\Omega)$ of Problem~\eqref{eq:1}, we consider its action on a test function $v \in H^2(\bar{\Omega})$. For $u$ sufficiently regular away from $\gamma$, integration by parts yields:
\[
-\int_{\Omega} \Delta u \cdot v \rmd\bfx = \int_{\Omega} \nabla u \cdot \nabla v \rmd\bfx - \int_{\gamma^+} \frac{\partial u}{\partial n^+} v \rmd\Haus^{n-1} - \int_{\gamma^-} \frac{\partial u}{\partial n^-} v \rmd\Haus^{n-1}.
\]

Noting that the unit normal vectors on opposite crack faces satisfy $\bfn^+ = -\bfn^- = \bfn$, and using the boundary conditions $\partial u/\partial n = 0$ on $\gamma^{\pm}$, we obtain the fundamental relation:
\begin{equation}\label{eq:3}
-\int_{\Omega} \Delta u \cdot v \rmd\bfx = -\int_{\Omega} u \Delta v \rmd\bfx - \int_{\gamma} [u] \frac{\partial v }{\partial n} \rmd\Haus^{n-1} \quad \forall v \in H^2(\bar{\Omega}).%\qquad (3)
\end{equation}
where $\Haus^{n-1}$ is the $(n-1)$-dimensional Hausdorff measure (surface measure on $\gamma$).

\subsection{The SBV framework and distributional derivatives}

The SBV (Special functions of Bounded Variation) framework provides the rigorous mathematical foundation for our approach. For a function $u \in SBV(\Omega)$, the distributional derivative $Du$ is a measure that decomposes into two parts:
\begin{equation}\label{eq:4}
Du = \nabla u \, \clL^n + [u] \, v _u \, \Haus^{n-1} \llcorner J_u
%\tag{4}
\end{equation}
where:
\begin{itemize}
\item $\nabla u$ is the \emph{approximate gradient} defined $\clL^n$-almost everywhere (the classical gradient where the function is smooth). This represents the regular, continuous part of the deformation.

\item $\clL^n$ denotes the $n$-dimensional Lebesgue measure (volume measure in $\bbR^n$). The term $\nabla u \, \clL^n$ is the part of the derivative that is absolutely continuous with respect to the Lebesgue measure.

\item $[u] = u^+ - u^-$ is the \emph{jump} of $u$ across the discontinuity set $J_u$ (the crack surface $\gamma$).\linebreak It represents the magnitude of the discontinuity.

\item $v _u$ is the \emph{unit normal vector} to the jump set $J_u$, oriented from the minus side to the plus side.

\item The notation $\Haus^{n-1} \llcorner J_u$ indicates the restriction of this measure to the jump set.

\item The product $[u] \, v _u \, \Haus^{n-1} \llcorner J_u$ is a \emph{singular measure} concentrated on the discontinuity surface. It captures the effect of the jump.
\end{itemize}

This decomposition allows us to define the extended solution $\bar{u}$ as a distribution on $\bar{\Omega}$:

\begin{defi}[Distributional Gradient]
The gradient of $\bar{u}$, denoted $\nabla \bar{u}$, is defined as the functional:
\begin{equation}\label{eq:5}
\langle \nabla \bar{u}, v \rangle = -\int_{\bar{\Omega}} \bar{u} \, \nabla \cdot v \rmd\bfx + \int_{\gamma} [u] (v \cdot \bfn)
\,d\Haus^{n-1} \quad \forall v \in H^1(\bar{\Omega}).
\end{equation}
\end{defi}

\begin{defi}[Distributional Laplacian]
The Laplacian of $\bar{u}$, denoted $\Delta \bar{u}$, is defined as the functional:
\begin{equation}\label{eq:6}
\langle \Delta \bar{u}, v \rangle = \int_{\bar{\Omega}} \bar{u} \, \Delta v \rmd\bfx + \int_{\gamma} \left(\left[\frac{\partial u}{\partial n}\right] v - [u] \frac{\partial v}{\partial n} \right) d\Haus^{n-1} \quad \forall v \in H^2(\bar{\Omega}). 
%\tag{6}
\end{equation}
\end{defi}

\subsection{Justification of the second-Order formulation}

In classical finite elements for the Poisson equation, one uses the symmetric form $\int \nabla u \cdot \nabla v$. That requires $u \in H^1$ and leads to continuity across elements. Here, $u$ is discontinuous across $\gamma$, and its gradient is not square-integrable on $\gamma$. To distribute the discontinuity into a surface term, we integrate by parts twice, transferring \emph{two derivatives} onto the test function $v $. This yields $\Delta v $ and a boundary term involving $[u] \partial v /\partial n$. This is not an arbitrary choice but a direct consequence of the distributional definition of $\Delta \bar{u}$ in the dual of $H^2$. It allows the jump to appear naturally as a measure, without modifying the mesh.


\subsection{Regularity considerations}

It is essential to clarify the regularity of the different fields involved:
\begin{itemize}
\item The \emph{physical solution} $u$ on the cracked domain $\Omega$ belongs to $H^1(\Omega)$.
\item The \emph{extended field} $\bar{u}$ is a \emph{distribution} on $\bar{\Omega}$. It belongs to $H^{-2}(\bar{\Omega})$, the dual space of $H^2_0(\bar{\Omega})$.
\item The \emph{test functions} $v$ are required to be in $H^2_0(\bar{\Omega})$ because the weak form involves $\Delta v$.
\end{itemize}

This is standard in distribution theory: we test the equation $\Delta \bar{u} = f$ against smooth functions~$v$, yielding:
\[
\langle \Delta \bar{u}, v \rangle = \langle f, v \rangle,
\]
where $\langle\,\cdot\,,\,\cdot\,\rangle$ denotes the duality pairing between $H^{-2}$ and $H^2_0$.



\subsection{The pseudo-continuous weak form}

Incorporating the boundary condition $\partial u/\partial n|_{\gamma} = 0$ (which implies $[\partial u/\partial n] = 0$), and starting from the definition of the distributional Laplacian, we obtain the pseudo-continuous weak form: find $\bar{u}$ such that:
\begin{equation}\label{eq:7}
-\int_{\bar{\Omega}} \bar{u} \, \Delta v \rmd\bfx = \int_{\bar{\Omega}} f v \rmd\bfx + \int_{\gamma} [\bar{u}] \frac{\partial v}{\partial n} \rmd\Haus^{n-1} \quad \forall v \in H^2_0(\bar{\Omega}). 
%\tag{7}
\end{equation}

\section{Mixed finite element discretization}\label{sect4}

\subsection{The mixed variational formulation}

The key insight of our numerical approach is to recognize that the jump $[u]$ is not a property of the regular part $\bar{u}_h$ itself, but an independent kinematic variable. We therefore introduce a mixed two-field variational scheme:
\begin{itemize}
\item A regular part $\bar{u}_h \in V_h \subset H_0^2(\bar{\Omega})$ approximated with $C^1$-continuous elements on the uncracked domain.
\item An independent jump field $j_\gamma$ defined only on the crack surface $\gamma$, representing the displacement discontinuity $[u]$.
\end{itemize}

The physical displacement $u$ is then reconstructed from the regular part and the independent jump field as:
\[
u(\bfx) = \underbrace{\bar{u}_h(\bfx)}_{\text{regular part} \,+\, \underbrace{j_\gamma(\bfx_\gamma) \, H_\gamma(\bfx)}_{\text{jump part}},}
\]
where $H_\gamma$ is a Heaviside function aligned with $\gamma$ (1 on the $+$ side and 0 on the $-$ side).

\subsubsection*{Theoretical status}

In the SBV framework, any function $u \in SBV(\Omega)$ decomposes uniquely into a regular $H^1$ component and a jump concentrated on $J_u$. Our approximation $\bar{u}_h \in H_0^2(\bar{\Omega})$ (continuous derivatives) approximates the regular part, and $j_\gamma$ approximates the jump amplitude. The equality is understood \emph{in the limit of mesh refinement} within a Galerkin approximation.

\subsection{The Discrete Variational Problem}

Let $V_h \subset H_0^2(\bar{\Omega})$ be a finite-dimensional subspace spanned by $C^1$-continuous basis functions: $V_h = \Span\{\phi_1, \phi_2, \dots, \phi_N\}$. Let $J_h \subset L^2(\gamma, \Haus^{n-1})$ be a finite-dimensional space for the jump field, spanned by basis functions $\{\psi_1, \psi_2, \Span\dots, \psi_M\}$ defined on $\gamma$. The associated norm is:
\[
\left\|j_\gamma\right\|^2_{L^2(\gamma)} = \int_\gamma j_\gamma^2 \rmd\Haus^{n-1}(\gamma).
\]

The mixed discrete problem is: find $(\bar{u}_h, j_\gamma) \in V_h \times J_h$ such that:
\begin{equation}\label{eq:8}
-\int_{\bar{\Omega}} \bar{u}_h \Delta v _h \rmd\bfx = \int_{\bar{\Omega}} f v _h \rmd\bfx + \int_{\gamma} j_\gamma \frac{\partial v _h}{\partial n} \rmd\Haus^{n-1}(\gamma) \quad \forall v _h \in V_h,
%\qquad (9)
\end{equation}
supplemented by an additional equation that links $j_\gamma$ to $\bar{u}_h$ via the traction condition on $\gamma$. For the traction-free crack, this additional equation enforces that the variational derivative with respect to $j_\gamma$ vanishes. The equation that enforces the traction-free condition is written as:
\begin{equation}\label{eq:9}
\int_{\gamma} \frac{\partial \bar{u}_h}{\partial n} \, \delta j_\gamma \rmd\Haus^{n-1}(\gamma) = 0 \quad \forall \delta j_\gamma \in J_h,
% \qquad (10)
\end{equation}
where $\delta j_\gamma$ represents a variation (or test function) for the jump field $j_\gamma$. This is the weak form of the condition $\partial u/\partial n = 0$ on $\gamma^{\pm}$.

\subsection{Linear system assembly}

Expressing $\bar{u}_h = \sum_{j=1}^N U_j \phi_j$ and $j_\gamma = \sum_{k=1}^M I_k \psi_k$, and using test functions $v _h = \phi_i$ and $\delta j_\gamma = \psi_l$, we obtain the coupled linear system:
\begin{equation}\label{eq:10}
\begin{bmatrix}
\bfK_{uu} & \bfK_{uj} \\
\bfK_{ju} & \bfK_{jj}
\end{bmatrix}
\begin{bmatrix}
\bfU\\
\bfJ
\end{bmatrix} =
\begin{bmatrix}
\bfF \\
\bfZO
\end{bmatrix},
%\qquad (11)
\end{equation}
where:
\begin{equation}\label{eq:11}
\begin{aligned}
(K_{uu})_{ij} &= -\int_{\bar{\Omega}} \phi_j \Delta \phi_i \rmd\bfx, \\[4pt]
(K_{uj})_{ik} &= \int_{\gamma} \psi_k \frac{\partial \phi_i}{\partial n} \rmd\Haus^{n-1}(\gamma), \\[4pt] (K_{ju})_{lj} &= \int_{\gamma} \frac{\partial \phi_j}{\partial n} \psi_l \rmd\Haus^{n-1}(\gamma), \\[4pt] (K_{jj})_{lk} &= \text{(terms from the traction condition)}, \\[4pt]
F_i &= \int_{\bar{\Omega}} f \phi_i \rmd\bfx.
\end{aligned}% \qquad (12)
\end{equation}

\subsection{Implementation remarks}
\begin{description}
\item[$C^1$-continuity] The formulation necessitates $C^1$-continuous elements (e.g., Hermite, Argyris, or isogeometric analysis elements) because the weak form involves second-order derivatives of test functions.
\item[Localized effect] The coupling terms $\int_{\gamma} \psi_k \frac{\partial \phi_i}{\partial n} \rmd\Haus^{n-1}(\gamma)$ are non-zero only for basis functions $\phi_i$ whose supports intersect the crack $\gamma$. This leads to a highly localized modification of the stiffness matrix.
\item[Simple meshing] The mesh is generated for the uncracked domain $\bar{\Omega}$. The crack $\gamma$ is simply an internal surface where these coupling integrals are computed; the mesh does not need to conform to it.
\end{description}

\section{Discussion}\label{sect5}

\subsection{Comparison with XFEM and crack tip singularities}

XFEM's ability to capture the $1/\sqrt{r}$ stress singularity via asymptotic enrichment is a key advantage that we acknowledge~\cite{4}. Our current formulation does not model this singularity explicitly. It is better suited for situations where:
\begin{itemize}
\item The crack tip lies outside the region of interest (e.g., problems dominated by interface behavior).
\item The far-field solution is the primary concern.
\item Multiple or intersecting cracks make enrichment strategies cumbersome.
\end{itemize}

We note that our method captures singularities implicitly through the solution of the weak form with jump terms, but it does not guarantee the correct asymptotic behavior near the tip without sufficient mesh refinement. A promising direction is a hybrid approach: using our method for the jump away from the tip and XFEM-type enrichment near the tip, combining the strengths of both techniques.

\subsection{Extensions to interfaces and ontact}

The distributional framework is general and can be extended to other types of discontinuities:
\begin{description}
\item[{Material interfaces}] continuity of displacement $[u] = 0$ with possible jumps in flux. This requires modifying the surface integral terms to enforce displacement continuity while allowing flux jumps determined by material properties.
\item[{Contact problems}] inequality constraints on the jump field (no interpenetration) and the introduction of contact pressure as a Lagrange multiplier. The weak form would be extended to include variational inequalities and additional Lagrange multiplier fields.
\end{description}

\subsection{Towards Evolving Discontinuities}

The extension to propagating cracks is natural within the same distributional setting:
\begin{itemize}
\item Represent the moving crack $\gamma(t)$ via a level-set function $\phi(\bfx,t)$.
\item The jump field $j_{\gamma(t)}$ evolves according to a fracture criterion (e.g., energy release rate).
\end{itemize}

\subsubsection{Handling moving bases for $J_h(t)$}

When the crack $\gamma(t)$ evolves, the support of the jump field changes~\cite{3}. The space $J_h(t) \subset L^2(\gamma(t), \Haus^{n-1})$ is therefore time-dependent. We propose two practical approaches:
\begin{description}
\item [Explicit surface remeshing] At each time step, $\gamma(t)$ is discretized independently by a surface mesh, and a new basis $\{\psi_k(t)\}$ is constructed. This is analogous to standard moving-mesh methods but without re-meshing the bulk domain.
\item [Level-set + fixed basis] Define a fixed set of shape functions $\{\psi_k(\bfx)\}$ on the background mesh of~$\bar{\Omega}$. At each time step, the integrals over $\gamma(t)$ are evaluated via a level-set function $\phi(\bfx,t)$ (e.g., $\gamma(t)=\{\phi=0\}$). The restriction of $\psi_k$ to $\gamma(t)$ provides a natural basis for $J_h(t)$. This avoids rebuilding the basis entirely.
\end{description}

A detailed implementation and comparison will be presented in a follow-up paper focused on numerical aspects.

\section{Conclusion}\label{sect6}

The pseudo-continuous media framework represents a significant shift in computational mechanics. By leveraging distribution theory and SBV spaces, it transforms problems with complex geometric discontinuities into mathematically equivalent mixed formulations on simplified domains. The complexity is moved from the geometry (meshing) to the formulation (incorporating jumps via independent fields).

Our mixed finite element formulation, grounded in rigorous distribution theory and the SBV framework, provides a solid theoretical foundation for treating strong discontinuities without conforming meshes. The method:
\begin{itemize}
\item Uses a non-conforming mesh without geometric alignment to the discontinuity.
\item Maintains mathematical consistency by treating the discontinuity intrinsically through an independent jump field.
\item Avoids the need for additional enriched degrees of freedom required by XFEM, though with the trade-off of requiring $C^1$ continuity.
\item Offers a clear pathway for extensions to interfaces, contact, and propagating cracks.
\end{itemize}

The theoretical developments presented here pave the way for a new class of numerical methods for problems with strong discontinuities. The numerical implementation and validation of the method will be the subject of a forthcoming paper, which will demonstrate its practical capabilities and confirm the theoretical predictions
%~\cite{2,3,4}
.



\printCOI

\printbibliography
\end{document}