Data-Driven Numerical Method Construction Via Physics-Informed Hard-Constrained Manifold Learning: Finite Elements as a Demonstration

1. Introduction

Data-driven learning and Artificial Intelligence (AI) have been increasingly used in scientific computing to solve partial differential equations (PDEs), with most existing approaches focusing on approximating solutions to prescribed PDE problems. This raises a different question: can data be used not only to approximate solutions, but also to construct a reusable numerical method itself? In the present work, we explore this possibility through finite elements and refer to the resulting paradigm as data-driven numerical method construction. Specifically, instead of constructing an element stiffness matrix entirely through analytical derivations based on prescribed interpolation functions, variational formulations, numerical integration rules, stabilization strategies, and remedies for numerical pathologies, we consider whether high-quality finite elements can be constructed through a direct data-driven paradigm while retaining the fundamental mechanical requirements of finite element theory. More importantly, the resulting data-driven finite element should retain the defining computational properties of a conventional finite element: it should depend only on local element information, remain independent of the global loading and boundary conditions, and be reusable and assemblable in different finite-element models.

Finite Element Method (FEM) provides a representative setting for investigating data-driven numerical method construction because its basic computational components, i.e., the finite elements, are inherently local, reusable, and assemblable. Over more than half a century of intensive development, FEM has become one of the most important numerical tools in engineering science, with broad applications in structural mechanics, fluid dynamics, electromagnetics, and multiphysics problems (Hughes, 2000; Zienkiewicz et al., 2013). The accuracy, convergence, and robustness of finite element analyses depend strongly on the quality of the underlying element technology. Traditional finite elements are usually constructed through analytical derivations based on interpolation functions, variational principles, and numerical quadrature (Bathe, 2014). Over the past decades, substantial progress has been made in the development of isoparametric elements, conforming and nonconforming elements, mixed elements, and various enhanced formulations. To alleviate numerical difficulties such as shear locking, volumetric locking, curvature locking, and hourglass modes, techniques including reduced integration, enhanced assumed strain methods, assumed natural strain methods, and stabilization schemes have been developed (Bathe and Dvorkin, 1986; Simo and Rifai, 1990; Simo and Armero, 1992). These developments have produced a large body of mature element formulations and corresponding high-quality element stiffness data. They therefore provide useful reference data for examining whether such element-level numerical characteristics can be learned within a mechanically admissible data-driven framework. At the same time, the construction of high-quality elements remains highly expertise-intensive. This challenge becomes particularly acute for three-dimensional solid-shell elements, which must simultaneously address thickness Poisson locking, transverse shear locking, membrane locking, and hourglass instabilities (Abed-Meraim and Combescure, 2002; Hauptmann and Schweizerhof, 1998; Reese, 2007).

Existing AI-assisted approaches in scientific computing can be distinguished from the objective pursued here according to what is learned. At the level of complete boundary-value problems or computational subdomains, existing studies usually focus on learning solution fields or local response operators, which can be taken as data-driven PDE solutions. Representative examples include PINNs (Haghighat et al., 2021; Raissi et al., 2019), HiDeNN-FEM (Liu et al., 2023), DeepFEM (Dong et al., 2023), the deep finite element method (Xiong et al., 2025), and the neural-operator element method (Ouyang et al., 2025). These methods are effective for solving or accelerating prescribed PDE problems, or for constructing reusable operators for particular classes of computational regions. Nevertheless, their learned quantities are generally associated with prescribed PDE problems, computational domains, boundary conditions, coefficient fields, or governing equations. They therefore do not directly provide a finite element determined only by local element geometry and material parameters that can be assembled into a finite element system in the same manner as a conventional finite element. At the element level, another class of studies introduces data-driven techniques into particular components of established finite-element formulations while retaining the classic construction framework under variational formulations and interpolation functions. Hence, they can be taken as AI-assisted classical numerical methods. For example, Jung et al. (2020) used neural networks to generate the strain-displacement matrix, whereas the subsequent stiffness construction still relies on isoparametric mapping, Gaussian quadrature, and element-specific correction procedures. Related studies have used nodal coordinates and material parameters to assist the design of four-node quadrilateral or eight-node hexahedral elements, or to approximate the element stiffness matrix for plane-stress elements (Jia et al., 2020; Nath et al., 2023; Subhash et al., 2024). Other studies have introduced data-driven models into specific components of conventional numerical formulations. Tandale et al. (2022) used neural networks to approximate element internal forces and tangent stiffness matrices, while Kim et al. (2026) and Yu et al. (2023) developed learned quadrature schemes for beam and enriched solid elements, respectively. Enabe and Provasi (2025) combined the virtual element method with deep learning to predict displacement fields. Jung et al. (2022) used deep learning to identify optimal bending directions for alleviating shear locking, and Park et al. (2026) subsequently improved the iterative procedure through a load-decomposition strategy. Therefore, the two lines, which are data-driven PDE solutions and AI-assisted classical FEM construction, are different from the goal of this paper. Here, the learned object is the finite element itself: a local stiffness operator determined only by element geometry and material parameters and intended to be assembled and reused in different structural problems.

Direct data-driven construction of a numerical method requires more than accurate data fitting: the learned object must satisfy the mathematical and physical properties required for a valid numerical method. For finite elements, the availability of high-quality element stiffness data does not by itself make direct data-driven finite-element construction feasible. An element stiffness matrix is not an arbitrary symmetric matrix whose entries can be determined solely by minimizing the discrepancy from reference data. Finite element theory imposes physics-informed hard constraints on admissible element stiffness matrices, including matrix symmetry, the rigid-body zero-energy constraint, the constant strain or constant generalized-strain consistency constraint, and the energy stability constraint. These requirements are directly related to the fundamental consistency, stability, and convergence properties of finite element calculations. Therefore, a small stiffness-matrix prediction error alone does not imply a mechanically admissible finite element: even a small violation of the null-space or consistency conditions can destroy basic element properties. The direct data-driven FEM construction problem should instead be formulated as a regression problem subject to the prescribed physics-informed hard constraints. Existing physics-informed learning approaches commonly introduce physical requirements into the training objective through residual or penalty terms (Karniadakis et al., 2021; Lu et al., 2021; Raissi et al., 2019). Such soft-constraint approaches can reduce constraint violations but cannot guarantee their exact satisfaction. Constraint-projection-based prediction-correction (Valente et al., 2025) provides another strategy for satisfying hard constraints, but the correction of the predicted quantities may compromise the original approximation capability of the learning model. Therefore, the central difficulty is not merely how to learn the element stiffness matrix from data, but how to restrict data-driven learning to the mechanically admissible stiffness-matrix space.

To realize this data-driven FEM construction, the present work develops a physics-informed hard-constrained manifold learning (HCML) framework, referred to as HCML-FE. Instead of constructing finite elements through explicit analytical formulations, element construction is formulated as a mapping from local element geometry and material parameters to a mechanically admissible element stiffness matrix. The key idea is to represent the admissible element stiffness matrices through a low-dimensional parameterization of the constrained stiffness-matrix manifold rather than learning the complete element stiffness matrix directly in its original high-dimensional space. Specifically, the rigid-body zero-energy and constant strain or constant generalized-strain consistency constraints are analytically transformed into linear relations among the stiffness-matrix entries. Symbolic Gaussian elimination is then employed to identify the independent stiffness parameters and derive an element stiffness matrix reconstruction operator. Within the resulting HCML framework, a neural network is used only to learn the nonlinear mapping from the element input parameters to the independent stiffness parameter vector, while the complete element stiffness matrix is subsequently reconstructed through the deterministic reconstruction operator. Consequently, the prescribed linear physics-informed hard constraints are satisfied by construction, rather than approximately imposed through penalty terms in the training objective.

The proposed HCML framework is not tied to a specific element type or source of reference data. In the present study, the reference stiffness matrices of four-node plane-stress elements are generated using an independently developed finite-element program, whereas those of the eight-node solid-shell element are extracted from commercial finite-element software. An important structural result is that the prescribed linear constraints reduce the stiffness description to only 3 independent parameters for the four-node plane-stress element and 78 for the eight-node solid-shell element. For each element type, the number and locations of these independent stiffness entries are fixed and do not depend on the admissible element geometry or material parameters. The proposed framework is validated through rigid-body zero-energy tests, constant-strain and constant generalized-strain patch tests, mesh-refinement and mesh-distortion studies, and a series of two- and three-dimensional structural benchmarks. The numerical results show that the resulting HCML-FEs exhibit convergence behavior and characteristic element responses consistent with the corresponding reference elements and can learn element technologies encoded in high-quality stiffness data. The learned solid-shell element can also be assembled with conventional solid elements in mixed-element analysis. These results demonstrate the data-driven numerical-method construction: physics defines mechanical admissibility, and data determine the element-specific numerical characteristics, i.e., element technologies.

The contribution of this work is threefold. First, finite-element construction is formulated as a data-driven numerical-method construction problem rather than a prescribed PDE-solution problem, whose differences are shown in Figure 1. Second, the linear mechanical constraints on an admissible element stiffness matrix are reduced analytically to a low-dimensional independent-parameter representation, enabling exact satisfaction of these constraints during learning. Third, the framework is demonstrated on element formulations with different numerical technologies encoded in stiffness operators, and the resulting learned elements retain the locality, reusability, and assemblability required of conventional finite elements.

The remainder of this paper is organized as follows. Section 2 introduces the general formulation of the proposed direct data-driven FEM construction framework. Section 3 specializes the framework to four-node plane-stress quadrilateral elements and eight-node solid-shell elements. Section 4 describes the data generation procedure, neural-network architecture, training strategy, and prediction-accuracy assessment. Section 5 presents the implementation details and the whole flowchart. Section 6 verifies the performance of the proposed HCML-based finite elements. Finally, Section 7 concludes the paper and discusses future research directions.

2. General Framework for Data-Driven Finite Element Method Construction

2.1. Problem Statement and Physics-Informed Hard Constraints of Element Stiffness Matrix

To establish the general framework to data-driven FEM construction, we first introduce the element input space, output space, and admissible element stiffness matrix set.

Definition 1 (Element Input Parameter Space): For a given element type and degree-of-freedom definition, the element input parameters include the set of geometric and material parameters. The geometric parameters are represented by the vector $\mathbf{g}\in \mathcal{G}\subseteq {\mathbb{R}}^{{\mathrm{n}}_{\mathrm{g}}}$, where $\mathcal{G}$ denotes the admissible geometric configuration space and ${\mathrm{n}}_{\mathrm{g}}$ denotes the number of geometric parameters, including nodal coordinates, thickness parameters, or other geometry-related quantities. The material properties are represented by the vector $\mathbf{m}\in \mathcal{M}\subseteq {\mathbb{R}}^{{\mathrm{n}}_{\mathrm{m}}}$, where $\mathcal{M}$ denotes the admissible material parameter space and ${\mathrm{n}}_{\mathrm{m}}$ denotes the number of material parameters, including Young's modulus, Poisson's ratio, or others. The complete element input vector is defined as $\mathbf{x}=\left(\mathbf{g},\mathbf{m}\right)\in \mathcal{X}$, where the element input space is given by $\mathcal{X}=\mathcal{G}\times \mathcal{M}\subseteq {\mathbb{R}}^{{\mathrm{n}}_{\mathrm{g}}+{\mathrm{n}}_{\mathrm{m}}}$. Each admissible element input parameter $\mathbf{x}\in \mathcal{X}$ determines a corresponding element stiffness matrix.

Definition 2 (Symmetric Element Stiffness Matrix Space): For an element with degrees of freedom derived from the standard Galerkin formulation of linear elasticity, the element stiffness matrix is symmetric. The ambient element stiffness matrix space is therefore defined as

$$\mathcal{Y}≔{\mathbb{S}}^{n_d}=\left\{\mathbf{K}\in {\mathbb{R}}^{n_d\times n_d}\mid {\mathbf{K}}^{\mathrm{T}}=\mathbf{K}\right\}.$$

(1)

Finite element theory imposes physics-informed hard constraints that restrict feasible element stiffness matrix to a proper subset of ${\mathbb{S}}^{n_d}$ to guarantee the convergence and numerical stability of finite element calculations. These constraints include rigid-body zero-energy constraint, constant strain or constant generalized-strain consistency constraint, and energy stability constraint

Definition 3 (Mechanically Admissible Element Stiffness Matrix): The mechanically admissible element stiffness matrix set ${\mathcal{K}}_{\mathrm{a}\mathrm{d}}$ is defined by

$${\mathcal{K}}_{\mathrm{a}\mathrm{d}}≔\left\{\mathbf{K}\in \mathcal{Y}\mid \mathbf{K}\hspace{0.25em}\mathrm{s}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{s}\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{s}\mathbf{C}\mathbf{o}\mathbf{n}.1-\mathbf{C}\mathbf{o}\mathbf{n}.3\right\}.$$

(2)

Unless otherwise specified, $\mathbf{K}$ hereafter denotes an admissible element stiffness matrix, i.e., $\mathbf{K}\in {\mathcal{K}}_{\mathrm{a}\mathrm{d}}.$

The physics-informed hard constraints “$\mathbf{C}\mathbf{o}\mathbf{n}.1-\mathbf{C}\mathbf{o}\mathbf{n}.3$” defined in Eq. (2) are specified as follows.

Constraint 1 (Rigid-body zero-energy constraint): Let

$$\mathbf{R}\left(\mathit{x}\right)=\left[{\mathit{r}}_1,{\mathit{r}}_2,\dots ,{\mathit{r}}_{n_r}\right]\in {\mathbb{R}}^{n_d\times n_r},$$

(3)

where the columns of $\mathbf{R}\left(\mathit{x}\right)$ form a basis of the rigid-body displacement space for $\mathit{x}\in \mathcal{X}$, $n_r$ is the number of rigid-body modes. For planar elasticity problems, $n_r=3$, whereas for three-dimensional elasticity problems, $n_r=6$. Since rigid-body motions should not generate strain energy, the element stiffness matrix must satisfy the following Constraint 1:

$$\mathbf{K}\mathbf{R}\left(\mathit{x}\right)=0.$$

(4)

Equivalently, the above matrix equation can be written in component form as:

$${\displaystyle \sum _{j=1}^{n_d}}{K}_{ij}{R}_{jk}\left(\mathit{x}\right)=0, i=1,\dots ,n_d, k=1,\dots ,n_r.$$

(5)

Eq. (3) and (4) provide linear homogeneous equality constraints on the entries of the element stiffness matrix.

Constraint 2 (Constant strain or constant generalized-strain consistency constraint): Let

$$\mathbf{S}\left(\mathit{x}\right)=\left[{\mathit{s}}_1,{\mathit{s}}_2,\dots ,{\mathit{s}}_{n_s}\right]\in {\mathbb{R}}^{n_d\times n_s},$$

(6)

where the columns of $\mathbf{S}\left(\mathit{x}\right)$ denote the nodal displacement vector for $\mathit{x}\in \mathcal{X}$ corresponding to constant strain or constant generalized strain states, where $n_s$ is the number of such vectors. For each column in $\mathbf{S}\left(\mathit{x}\right)$, the corresponding equivalent nodal force vector can be obtained from the theoretical stress field. The linear relation between nodal displacement and nodal force under constant strain or constant generalized-strain can be represented as

$${\mathbf{F}}_S\left(\mathit{x}\right)=\mathcal{P}\left(\mathit{x}\right)\mathbf{S}\left(\mathit{x}\right),$$

(7)

where $\mathcal{P}\left(\mathit{x}\right)$ is a geometry and material dependent linear operator, ${\mathbf{F}}_S\left(\mathit{x}\right)\in {\mathbb{R}}^{n_d\times n_s}$ denotes the equivalent nodal force matrix. Hence, Constraint 2 can be expressed as

$$\mathbf{K}\mathbf{S}\left(\mathit{x}\right)={\mathbf{F}}_S\left(\mathit{x}\right).$$

(8)

Equivalently, the above matrix equation can be written in component form as:

$${\displaystyle \sum _{j=1}^{n_d}}{K}_{ij}{S}_{jk}\left(\mathit{x}\right)=F_{S,ik}\left(\mathit{x}\right), i=1,\dots ,n_d, k=1,\dots ,n_s.$$

(9)

Eq. (6) and (7) provide linear non-homogeneous equality constraints on the entries of the element stiffness matrix.

Combining Constraints 1 and 2 constitutes the linear constraint system with respect to the element stiffness matrix $\mathbf{K}$:

$$\mathbf{K}\mathbf{C}\left(\mathit{x}\right)=\mathbf{D}\left(\mathit{x}\right), {\mathbf{K}}^{\mathrm{T}}=\mathbf{K},$$

(10)

Where

$$\mathbf{C}\left(\mathit{x}\right)=\left[\mathbf{R}\left(\mathit{x}\right)\hspace{0.25em} \mathbf{S}\left(\mathit{x}\right)\right], \mathbf{D}\left(\mathbf{x}\right)=\left[0\hspace{0.25em}{\mathbf{F}}_S\left(\mathit{x}\right)\right].$$

In the following, we refer to Eq. (10) as the physics-informed hard constraints

Constraint 3 (Energy stability constraint): The admissible element stiffness matrix must be positive semi-definite:

$${\mathbf{u}}^{\mathrm{T}}\mathbf{K}\mathbf{u}\ge 0, \forall \mathbf{u}\in {\mathbb{R}}^{n_d}.$$

(11)

Moreover, the strain energy should be strictly positive for any non-rigid-body displacement mode, i.e., $\forall \mathbf{u}\notin \mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left(\mathbf{R}\right),{\mathbf{u}}^{\mathrm{T}}\mathbf{K}\mathbf{u}>0.$ Equivalently, the null space of $\mathbf{K}$ should coincide with the rigid-body displacement space,

$$\mathrm{N}\mathrm{u}\mathrm{l}\mathrm{l}\left(\mathbf{K}\right)=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left(\mathbf{R}\right),$$

(12)

which leads to the rank condition

$$\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\mathbf{K}\right)=n_d-n_r.$$

(13)

This condition excludes spurious zero-energy modes. Furthermore, as proven by Stummel’s (1979) theory, satisfaction of the rigid-body zero-energy constraint and constant-strain patch test and the exclusion of spurious zero-energy modes can theoretically guarantee the convergence of finite element calculation under mesh refinement.

In the present HCML parameterization, Constraints 1 and 2 constitute the linear equality constraints that are enforced exactly by the reconstruction operator. Constraint 3 concerns stability and the exclusion of spurious zero-energy modes; it is examined separately through the spectrum of the reconstructed stiffness matrix.

Definition 4 (Direct Data-Driven FEM Construction Problem): Suppose that the high-quality finite elements implicitly establish an unknown target mapping ${\mathcal{F}}^{\mathrm{*}}:\mathcal{X}\to \mathcal{Y}$. The direct data-driven FEM construction problem is to construct a parameterized learned mapping ${\mathcal{F}}_{\mathsf{\theta }}:\mathcal{X}\to \mathcal{Y}$ to approximate the unknown target mapping ${\mathcal{F}}^{\mathrm{*}}$,and strictly satisfy:

$${\mathcal{F}}_{\mathsf{\theta }}\left(\mathit{x}\right)\in {\mathcal{K}}_{\mathrm{a}\mathrm{d}}\left(\mathit{x}\right), \forall \mathit{x}\in \mathcal{X}.$$

(14)

Here, the symbol $\mathsf{\theta }$ denotes the learnable parameters, such as the weights and biases of a neural- network.

Given the training data set with $N_s$ samples, denoted as $\mathcal{D}={\left\{\left({\mathit{x}}^{\left(s\right)},{\mathbf{K}}_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\left(s\right)}\right)\right\}}_{s=1}^{N_s},$ where ${\mathit{x}}^{\left(s\right)}\in \mathcal{X}$ is the $s$-th element input parameter and ${\mathbf{K}}_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\left(s\right)}\in {\mathcal{K}}_{\mathrm{a}\mathrm{d}}\left({\mathit{x}}^{\left(s\right)}\right)$ is the corresponding labeled element stiffness matrix extracted from high-quality reference data. The direct data-driven FEM construction problem can be modeled as a hard-constrained regression problem:

$$\begin{gathered} {\mathsf{\theta }}^{\mathrm{*}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{\theta }}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{0.33em}{\mathbb{E}}_{\left(\mathit{x},{\mathbf{K}}_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)\sim \mathcal{D}}\left[\mathcal{L}\left({\mathcal{F}}_{\mathit{\theta }}\left(\mathit{x}\right),{\mathbf{K}}_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)\right], \\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t} \mathrm{t}\mathrm{o}:{\mathcal{F}}_{\mathsf{\theta }}\left(\mathit{x}\right)\in {\mathcal{K}}_{\mathrm{a}\mathrm{d}}\left(\mathit{x}\right),\forall \mathit{x}\in \mathcal{X} \end{gathered}$$

(15)

Here, $\mathcal{L}\left(\cdot ,\cdot \right)$ is a loss function that measures the discrepancy between the predicted element stiffness matrix and the labeled element stiffness matrix.

As shown in the following section, ${\mathcal{K}}_{\mathrm{a}\mathrm{d}}\left(\mathit{x}\right)$ is a low-dimensional manifold embedded in the high-dimensional symmetric element stiffness matrix space ${\mathbb{S}}^{n_d}$, and is therefore referred to as the mechanically admissible element stiffness matrix manifold. From this perspective, the direct data-driven FEM construction problem can be formulated as a manifold learning problem.

2.2. Algebraic Representation of Physics-Informed Hard Constraints

To transform the mechanically admissible element stiffness matrix manifold into a computable and embeddable mathematical structure, the linear physics-informed hard constraints Eq. (10) are reformulated as explicit algebraic equations with respect to the element stiffness matrix entries. To be specific, the symmetric element stiffness matrix $\mathbf{K}\in {\mathbb{S}}^{n_d}$ is vectorized by stacking its lower triangular part under the prescribed row-wise order:

$$\mathit{k}\left(\mathit{x}\right)=\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{h}\left(\mathbf{K}\right)\in {\mathbb{R}}^{n_k}, n_k={\displaystyle \frac{n_d\left(n_d+1\right)}{2}},$$

(16)

where $\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{h}\left(\cdot \right)$ denotes the half-vectorization operator that collects the lower triangular entries, including the diagonal entries, of a symmetric matrix. Under this vectorization, the symmetry condition ${\mathbf{K}}^{\mathrm{T}}=\mathbf{K}$ is satisfied by construction.

Eq. (10) can then be equivalently rewritten as a linear algebraic system with respect to the vector $\mathit{k}\left(\mathit{x}\right)\in {\mathbb{R}}^{n_k}$:

$$\mathbf{A}\left(\mathit{x}\right)\mathit{k}\left(\mathit{x}\right)=\mathit{d}\left(\mathit{x}\right),$$

(17)

here, $\mathbf{A}\left(\mathit{x}\right)$ is the coefficient matrix determined by the displacement mode matrix $\mathbf{C}$, the prescribed vectorization rule $\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{h}\left(\cdot \right)$, and the element input parameter $\mathit{x}$. The vector $\mathit{d}\left(\mathit{x}\right)$ is the corresponding right-hand-side vector obtained from $\mathbf{D}\left(\mathit{x}\right)$. For a given element type, displacement mode matrix, and admissible input $\mathit{x}$, both $\mathbf{A}\left(\mathit{x}\right)$ and $\mathit{d}\left(\mathit{x}\right)$ can be constructed analytically.

Taking all entries in $\mathit{k}\left(\mathit{x}\right)$ as symbols, symbolic Gaussian elimination is then performed for Eq. (17), which analytically and precisely determines the independent stiffness entries and the dependent stiffness entries. The dependent stiffness entries are expressed as linear combinations of the independent stiffness parameters, leading to an analytical reconstruction relation for the complete element stiffness matrix. Since this procedure is carried out at the symbolic analytical level, for a given element type, the number and locations of the independent stiffness entries, as well as the analytical relations between the independent and dependent entries, are valid for any admissible element in $\mathcal{X}$. Hence, we can give the following definition.

Definition 5 (Independent Stiffness Parameter Vector). Based on the symbolic Gaussian elimination result of Eq. (17), independent stiffness entries from $\mathit{k}\left(\mathit{x}\right)$ can be identified and assembled into a vector $\mathit{p}\in {\mathbb{R}}^{n_p},$ where $n_p$ is the dimension of the null space of matrix $\mathbf{A}\left(\mathit{x}\right)$. Both independent parameter vector $\mathit{p}\in {\mathbb{R}}^{n_p}$ and number $n_p$ are fixed for any $\mathit{x}\in \mathcal{X}$. The general solution of the constraint system Eq. (10) can be written as

$$\mathit{k}\left(\mathit{x}\right)=\mathbf{B}\left(\mathit{x}\right)\mathit{p}+\mathit{b}\left(\mathit{x}\right),$$

(18)

where $\mathbf{B}\left(\mathit{x}\right)\in {\mathbb{R}}^{n_k\times n_p}$ is a basis matrix of the solution space, and $\mathit{b}\left(\mathit{x}\right)\in {\mathbb{R}}^{n_k}$ is a particular solution. Both $\mathbf{B}\left(\mathit{x}\right)$ and $\mathit{b}\left(\mathit{x}\right)$ are derived by solving the linear constraint system symbolically.

Definition 6 (Element Stiffness Matrix Reconstruction Operator). Based on the symbolic representation of the solution space, a deterministic element stiffness matrix reconstruction operator is defined as

$$\mathcal{Q}\left(\mathit{p};\mathit{x}\right)={\mathrm{m}\mathrm{a}\mathrm{t}}_{\mathrm{s}\mathrm{y}\mathrm{m}}\left(\mathbf{B}\left(\mathit{x}\right)\mathit{p}+\mathit{b}\left(\mathit{x}\right)\right),$$

(19)

where ${\mathrm{m}\mathrm{a}\mathrm{t}}_{\mathrm{s}\mathrm{y}\mathrm{m}}\left(\cdot \right)$ denotes the inverse operation of $\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{h}\left(\cdot \right)$, which reconstructs a symmetric matrix from its independent lower triangular entries.

For the plane-stress elements and three-dimensional solid-shell elements considered in this work, this matrix reconstruction operator can be written

$${\mathbf{K}}^{2\mathrm{D}}={\mathcal{Q}}_{2\mathrm{D}}\left(\mathit{p};\mathit{x}\right)$$

(20)

or

$${\mathbf{K}}^{3\mathrm{D}}={\mathcal{Q}}_{3\mathrm{D}}\left(\mathit{p};\mathit{x}\right)$$

(21)

By the above derivation and construction, the reconstructed element stiffness matrix satisfies linear physics-informed hard constraints Eq. (10). As a result, we successfully transform the abstract mechanically admissible stiffness matrix manifold into an explicit and computable algebraic structure.

2.3. Hard-Constrained Manifold Learning (HCML) Framework

Based on Section 2.2, the data-driven FEM construction problem can be solved by the manifold learning framework. That is to say, instead of directly predicting the full element stiffness matrix, the neural-network predicts only the independent parameter vector $\mathit{p}$, which is subsequently mapped to the complete element stiffness matrix through the deterministic matrix reconstruction operator $\mathcal{Q}$. In this learning framework, the operator $\mathcal{Q}$ is embedded after the neural-network output layer and can be referred to as the hard-constrained reconstruction operator.

Let ${\mathit{p}}_{\mathsf{\theta }}\left(\mathit{x}\right)={\mathcal{N}}_{\mathsf{\theta }}\left(\mathit{x}\right)$ denote the output of a neural-network, where ${\mathcal{N}}_{\mathsf{\theta }}$ is a parameterized neural mapping with learnable parameters $\mathsf{\theta }$, and the output ${\mathit{p}}_{\mathsf{\theta }}\left(\mathit{x}\right)\in {\mathbb{R}}^{n_p}$ is the predicted independent stiffness parameter vector.

For plane-stress elements, the complete parameterized learned mapping is reconstructed as

$${\mathcal{F}}_{\mathsf{\theta }}^{ml}\left(\mathit{x}\right)={\mathcal{Q}}_{2\mathrm{D}}\left({\mathit{p}}_{\mathsf{\theta }}\left(\mathit{x}\right);\mathit{x}\right),$$

(22)

For three-dimensional solid or solid-shell elements, the complete parameterized learned mapping is

$${\mathcal{F}}_{\mathsf{\theta }}^{ml}\left(\mathit{x}\right)={\mathcal{Q}}_{3\mathrm{D}}\left({\mathit{p}}_{\mathsf{\theta }}\left(\mathit{x}\right);\mathit{x}\right).$$

(23)

The supervised learning objective loss has been defined in Eq. (15). It can be evaluated either by comparing the reconstructed full element stiffness matrix with the corresponding reference stiffness matrix or by comparing the predicted independent stiffness parameter vector with its reference counterpart.

The proposed manifold learning framework reduces the learning task from a high-dimensional element stiffness matrix regression problem with constraints to a lower-dimensional independent-parameter unconstrained regression problem. This reduction simplifies the training and optimization process without penalty terms.

2.4. Soft-Constraint Baseline for Full Element Stiffness Matrix Prediction

To clarify the role of the proposed framework and further demonstrate the necessity of enforcing the prescribed linear physics-informed hard constraints, a representative soft-constraint-based learning strategy is introduced for comparison. In this soft-constraint baseline, the physical hard constraints are not embedded into the output structure of the neural-network. Instead, they are imposed approximately by adding penalty terms to the loss function. Different from the proposed physics-informed hard-constrained learning framework, the soft-constraint baseline directly predicts all entries of the complete element stiffness matrix ${\widehat{\mathbf{K}}}_{\mathsf{\theta }}$, rather than the independent stiffness parameter vector $\mathit{p}$.

Based on the above consideration, the total loss function of the soft-constraint baseline is defined as

$${\mathcal{L}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}=\mathcal{L}+{\lambda }_{\mathrm{T}}\mathcal{L}\left({\widehat{\mathbf{K}}}_{\mathsf{\theta }},{\left({\widehat{\mathbf{K}}}_{\mathsf{\theta }}\right)}^{\mathrm{T}}\right)+{\lambda }_{\mathrm{R}}{\mathcal{L}}_{\mathrm{R}}\left({\widehat{\mathbf{K}}}_{\mathsf{\theta }}\mathbf{R}\right)+{\lambda }_{\mathrm{S}}{\mathcal{L}}_{\mathrm{S}}\left({\widehat{\mathbf{K}}}_{\mathsf{\theta }}\mathbf{S}-{\mathbf{F}}_S\left(\mathit{x}\right)\right),$$

(24)

where the first term measures the discrepancy between the predicted element stiffness matrix and the reference element stiffness matrix, the second term penalizes the violation of matrix symmetry, the third term penalizes the residual forces induced by rigid-body zero-energy modes, and the fourth term penalizes the residual associated with the constant strain or constant generalized-strain consistency requirement. The weights ${\lambda }_{\mathrm{T}}$, ${\lambda }_{\mathrm{R}}$, and ${\lambda }_{\mathrm{S}}$ control the relative contributions of data fitting and the penalty-based enforcement of the prescribed physical constraints. The neural-network architecture, model parameters, training samples, and optimizer are kept identical to those used in the proposed hard-constrained manifold learning framework. This baseline is introduced not as an alternative finite-element formulation, but to examine whether small data-fitting errors are sufficient to preserve mechanically essential element properties.

3. Specialized Manifold Parameterization for Two Representative Element Types

This section specializes the general methodology discussed in Section 2 to two representative element types: four-node plane-stress quadrilateral element and 3D eight-node solid-shell element. For each element type, the canonical configuration, independent stiffness parameters, and element stiffness matrix reconstruction operator are derived explicitly.

3.1. Four-Node Plane-Stress Quadrilateral Element

3.1.1. Canonical Configuration and Element Input Parameters for Plane-Stress Quadrilateral Element

As shown in Figure 2, we establish the canonical configuration by fixing the first edge of the quadrilateral along the global x-axis with unit length. The nodal coordinates are given by

$${\mathbf{X}}_1^{2\mathrm{D}}=\left(0,0\right), {\mathbf{X}}_2^{2\mathrm{D}}=\left(1,0\right), {\mathbf{X}}_3^{2\mathrm{D}}=\left(x_3,y_3\right), {\mathbf{X}}_4^{2\mathrm{D}}=\left(x_4,y_4\right)$$

where ${\mathbf{X}}_i^{2\mathrm{D}}$ denotes the coordinate vectors of node $i$ in the two-dimensional space.

Each node possesses two displacement degrees of freedom, corresponding to the displacement components in the x- and y-directions. Accordingly, the element displacement vector is arranged as

$${\mathit{u}}^{2\mathrm{D}}={\left[u_1,u_2,u_3,u_4,v_1,v_2,v_3,v_4\right]}^{\mathrm{T}}.$$

(25)

where $u_i$​ and $v_i$ denote the x- and y-direction displacement components of node $i$, respectively. This ordering is adopted throughout the following element stiffness matrix parameterization.

The element is assumed to be composed of a linear elastic material characterized by Young's modulus $E$ and Poisson's ratio $\mu$. Since $E$ acts only as a scaling parameter for the linear elastic element stiffness matrix, it is normalized to $E=1$, while $\mu$ is retained as the material input. Hence, the element input parameters of the element can be defined as

$${\mathit{x}}^{2\mathrm{D}}=\left(\mu ,x_3,y_3,x_4,y_4\right)\in {\mathcal{X}}^{2\mathrm{D}}.$$

(26)

The subsequent construction of physics-informed hard constraints and stiffness reconstruction relations are based on this canonical configuration. In practical computations, an admissible four-node plane-stress quadrilateral element with arbitrary position, orientation, and size can be mapped to this canonical configuration through translation, rotation, and scaling transformations, as shown in Figure 3.

3.1.2. Hard Constraint System for the Quadrilateral Element Stiffness Matrix

Rigid-Body Zero-Energy Constraints

The rigid-body displacement modes consist of translations along the $x$- and $y$-directions and an in-plane rigid-body rotation, which are assembled into the matrix

$${\mathbf{R}}^{2\mathrm{D}}\left({\mathit{x}}^{2\mathrm{D}}\right)=\left[\begin{array}{c}{\mathbf{R}}_x^{2\mathrm{D}}\\ {\mathbf{R}}_y^{2\mathrm{D}}\end{array}\right]\in {\mathbb{R}}^{8\times 3},$$

(27)

where

$${\mathbf{R}}_x^{2\mathrm{D}}=\left[\begin{array}{ccc}1& 0& 0\\ 1& 0& 0\\ 1& 0& -y_3\\ 1& 0& -y_4\end{array}\right],{\mathbf{R}}_y^{2\mathrm{D}}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 1& 1\\ 0& 1& x_3\\ 0& 1& x_4\end{array}\right]$$

The first two columns of ${\mathbf{R}}^{2\mathrm{D}}\left({\mathit{x}}^{2\mathrm{D}}\right)$ represent rigid-body translations in the $x$-and $y$-directions, respectively, while the third column represents an in-plane rigid-body rotation.

The rigid-body zero-energy constraints can then be expressed by Eq. (4), this equation yields 24 homogeneous linear constraint equations with respect to the entries of the element stiffness matrix ${\mathbf{K}}^{2\mathrm{D}}$, although these equations are not necessarily linearly independent.

Constant Strain Consistency Constraint

The constant-strain displacement field of the plane-stress quadrilateral element is expressed in terms of the three independent strain components ${\epsilon }_{xx}$, ${\epsilon }_{yy}$, and ${\gamma }_{xy}$ as

$$u\left(x,y\right)={\epsilon }_{xx}x+{\displaystyle \frac{1}{2}}{\gamma }_{xy}y, v\left(x,y\right)={\epsilon }_{yy}y+{\displaystyle \frac{1}{2}}{\gamma }_{xy}x,$$

(28)

This representation corresponds to the constant strain vector

$$\mathsf{\epsilon }={\left[{\epsilon }_{xx},{\epsilon }_{yy},{\gamma }_{xy}\right]}^{\mathrm{T}}.$$

(29)

The nodal displacement modes associated with the three-unit constant-strain components are assembled into

$${\mathbf{S}}^{2\mathrm{D}}({\mathit{x}}^{2\mathrm{D}})=\left[\begin{array}{c}{\mathbf{S}}_x^{2\mathrm{D}}\\ {\mathbf{S}}_y^{2\mathrm{D}}\end{array}\right]\in {\mathbb{R}}^{8\times 3},$$

(30)

where

$${\mathbf{S}}_x^{2\mathrm{D}}=\left[\begin{array}{ccc}0& 0& 0\\ 1& 0& 0\\ x_3& 0& y_3/2\\ x_4& 0& y_4/2\end{array}\right], {\mathbf{S}}_y^{2\mathrm{D}}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 1/2\\ 0& y_3& x_3/2\\ 0& y_4& x_4/2\end{array}\right]$$

For each unit constant-strain mode, the corresponding equivalent nodal force vector is obtained from the theoretical constant stress field through boundary traction integration. The equivalent nodal force matrix associated with the three constant-strain modes is expressed as:

$${\mathbf{F}}_S^{2\mathrm{D}}({\mathit{x}}^{2\mathrm{D}})={\displaystyle \frac{1}{2\left(1-{\mu }^2\right)}}{\left({\mathbf{A}}_S^{2\mathrm{D}}\left(\mu \right)\hspace{0.17em}{\mathbf{N}}^{2\mathrm{D}}\left({\mathit{x}}^{2\mathrm{D}}\right)\otimes \left[\begin{array}{cc}1& -1\end{array}\right]\right)}^{\mathrm{T}}$$

(31)

where

$${\mathbf{A}}_S^{2\mathrm{D}}\left(\mu \right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& \left(1-\mu \right)/2\end{array}\right],$$

and

$${\mathbf{N}}^{2\mathrm{D}}\left({\mathit{x}}^{2\mathrm{D}}\right)=\left[\begin{array}{cccc}-y_4& -\mu \left(1-x_4\right)& -y_3& -\mu x_3\\ -\mu y_4& -\left(1-x_4\right)& -\mu y_3& -x_3\\ -\left(1-x_4\right)& -y_4& -x_3& y_3\end{array}\right],$$

Here, $\otimes$ denotes the Kronecker product.

The constant strain consistency constraint can then be expressed by Eq. (8). This relation provides 24 non-homogeneous linear equality constraints with respect to the entries of the element stiffness matrix ${\mathbf{K}}^{2\mathrm{D}}$, although these equations are not necessarily linearly independent.

By jointly enforcing the rigid-body zero-energy constraints and constant strain consistency constraint, the complete linear constraint system for the four-node plane-stress quadrilateral element can be written as

$${\mathbf{K}}^{2\mathrm{D}}{\mathbf{C}}^{2\mathrm{D}}\left({\mathit{x}}^{2\mathrm{D}}\right)={\mathbf{D}}^{2\mathrm{D}}\left({\mathit{x}}^{2\mathrm{D}}\right),$$

(32)

where

$${\mathbf{C}}^{2\mathrm{D}}\left({\mathit{x}}^{2\mathrm{D}}\right)=\left[{\mathbf{R}}^{2\mathrm{D}}\hspace{0.25em}{\mathbf{S}}^{2\mathrm{D}}\right]\in {\mathbb{R}}^{8\times 6}, {\mathbf{D}}^{2\mathrm{D}}\left({\mathit{x}}^{2\mathrm{D}}\right)=\left[0\hspace{0.25em}{\mathbf{F}}_S^{2\mathrm{D}}\right]\in {\mathbb{R}}^{8\times 6}.$$

This system provides the algebraic basis for deriving the independent stiffness parameters and stiffness matrix reconstruction operator for the four-node plane-stress quadrilateral element.

3.1.3. Identification of Independent Stiffness Parameters and Construction of Reconstruction Operator for the Quadrilateral Element

According to the vectorized linear constraint formulation introduced in Section 2.2, the linear physics-informed hard constraints Eq. (10) for the four-node plane-stress quadrilateral element can be rewritten as an algebraic system with respect to the independent entries of the symmetric element stiffness matrix. For the canonical configuration defined in Section 3.1.1, the solution space of this constraint system has dimension $n_p^{2\mathrm{D}}=3$.

The symbolic computation library SymPy in Python is employed to solve the linear constraint system Eq. (10) automatically. This symbolic solution process itself constitutes a rigorous mathematical derivation for determining the independent stiffness entries. For the present two-dimensional plane-stress quadrilateral element, the independent parameter vector is chosen as

$${\mathit{p}}^{2\mathrm{D}}={\left[{\widehat{K}}_{7,7},{\widehat{K}}_{7,8},{\widehat{K}}_{8,8}\right]}^{\mathrm{T}}.$$

(33)

The locations of the independent stiffness entries are shown in Figure 4. The number and locations of the independent entries are fixed for quadrilateral elements of arbitrary shapes, and the locations of the independent stiffness entries are shown in the following with any element type.

Based on the symbolic solution of the constraint system, the full element stiffness matrix can be reconstructed from the independent parameter vector by the element-specific hard-constrained reconstruction operator:

$${\mathbf{K}}^{2\mathrm{D}}={\mathcal{Q}}_{2\mathrm{D}}\left({\mathit{p}}^{2\mathrm{D}};{\mathit{x}}^{2\mathrm{D}}\right).$$

(34)

Equivalently, each stiffness entry can be expressed in the affine form

$$K_{i,j}^{2\mathrm{D}}=K_{i,j}^{\left({\mathbf{x}}^{2\mathrm{D}}\right)}+K_{i,j}^{\left(\mathrm{7,7}\right)}{\widehat{K}}_{\mathrm{7,7}}+K_{i,j}^{\left(\mathrm{7,8}\right)}{\widehat{K}}_{\mathrm{7,8}}+K_{i,j}^{\left(\mathrm{8,8}\right)}{\widehat{K}}_{\mathrm{8,8}},$$

(35)

here, ${\widehat{K}}_{7,7}$,${\widehat{K}}_{7,8}$ and ${\widehat{K}}_{8,8}$ are independent parameters, $K_{i,j}^{\left({\mathit{x}}^{2\mathrm{D}}\right)}$denotes the term from the non-homogeneous part in Eq. (10) and is irrelevant with independent parameters. The terms $K_{i,j}^{\left(\mathrm{7,7}\right)}$,$K_{i,j}^{\left(\mathrm{7,8}\right)}$and $K_{i,j}^{\left(\mathrm{8,8}\right)}$are coefficient functions associated with the three independent stiffness parameters, respectively. To illustrate the resulting reconstruction relation, representative stiffness entries are listed as follows：

$$\begin{gathered} K_{6,6}^{2\mathrm{D}}={\displaystyle \frac{x_3{y}_4+y_3{x}_4-y_3}{2y_3^2\left({\mu }^2-1\right)}}-{\displaystyle \frac{y_4^2}{y_3^2}}{\widehat{K}}_{8,8}, \\ K_{6,7}^{2\mathrm{D}}={\displaystyle \frac{\mu }{2\left({\mu }^2-1\right)}}-{\displaystyle \frac{y_4}{y_3}}{\widehat{K}}_{7,8}, \\ \mathrm{br}-\mathrm{to}-\mathrm{break} K_{6,8}^{2\mathrm{D}}=-{\displaystyle \frac{x_3}{2y_3\left({\mu }^2-1\right)}}-{\displaystyle \frac{y_4}{y_3}}{\widehat{K}}_{8,8}. \end{gathered}$$

(36)

The complete symbolic reconstruction relations are provided in Appendix B. This reduction from 36 independent entries of a general symmetric 8×8 matrix to only three free stiffness parameters quantifies how strongly the prescribed mechanical constraints restrict the admissible stiffness description.

3.2. Eight-Node Solid-Shell Element

3.2.1. Canonical Configuration and Element Input Parameters for the Solid-Shell Element

As illustrated in Figure 5, a three-dimensional eight-node solid-shell element constructed by midsurface extrusion is considered. The element is constructed by extending the two-dimensional quadrilateral canonical configuration along the thickness direction. The midsurface adopts the same in-plane canonical configuration as the four-node plane-stress quadrilateral element, while the thickness direction is aligned with the global z-axis. The total thickness of the element is denoted by $t$. The bottom surface consists of nodes 1–4 and is located at $z=-t/2$, whereas the top surface consists of nodes 5–8 and is located at $z=t/2$. Accordingly, the nodal coordinates of the bottom surface are given by

$${\mathbf{X}}_1^{3\mathrm{D}}=\left(0,0,-t/2\right), {\mathbf{X}}_2^{3\mathrm{D}}=\left(1,0,-t/2\right), {\mathbf{X}}_3^{3\mathrm{D}}=\left(x_3,y_3,-t/2\right), {\mathbf{X}}_4^{3\mathrm{D}}=\left(x_4,y_4,-t/2\right)$$

and those of the top surface are given by

$${\mathbf{X}}_5^{3\mathrm{D}}=\left(0,0,t/2\right), {\mathbf{X}}_6^{3\mathrm{D}}=\left(1,0,t/2\right), {\mathbf{X}}_7^{3\mathrm{D}}=\left(x_3,y_3,t/2\right), {\mathbf{X}}_8^{3\mathrm{D}}=\left(x_4,y_4,t/2\right)$$

where ${\mathbf{X}}_i^{3\mathrm{D}}$ denote the coordinate vectors of node $i$ in the three-dimensional space. Each node possesses three displacement degrees of freedom, corresponding to the displacement components in the x-, y-, and z-directions. Accordingly, the element displacement vector is arranged as

$${\mathit{u}}^{3\mathrm{D}}={\left[u_1,\dots ,u_8,{ v}_1,\dots ,v_8,{ w}_1,\dots ,w_8\right]}^{\mathrm{T}}$$

(37)

The Young's modulus is normalized to $\mathrm{E}=1$, while $\mathsf{\mu }$ is retained as the material input. Hence, the element input parameters of the solid-shell element can be defined as

$${\mathit{x}}^{3\mathrm{D}}=\left(\mu ,t,x_3,y_3,x_4,y_4\right)\in {\mathcal{X}}^{3\mathrm{D}}$$

(38)

3.2.2. Hard Constraint System for the Solid-Shell Element Stiffness Matrix

Rigid-Body Zero-Energy Constraints

The six rigid-body displacement modes consist of three translations along the x-, y-, and z-directions and three rotations about the x-, y-, and z-axes. They can be assembled into the matrix:

$${\mathbf{R}}^{3\mathrm{D}}({\mathit{x}}^{3\mathrm{D}})=\left[\begin{array}{c}\begin{array}{c}{\mathbf{R}}_x^{3\mathrm{D}}\\ {\mathbf{R}}_y^{3\mathrm{D}}\end{array}\\ {\mathbf{R}}_z^{3\mathrm{D}}\end{array}\right]\in {\mathbb{R}}^{24\times 6}$$

(39)

where

$${\mathbf{R}}_x^{3\mathrm{D}}=\left[1_8{0}_8{0}_8{0}_8{\mathit{z}}_{\mathrm{n}}-{\mathit{y}}_{\mathrm{n}}\right]\in {\mathbb{R}}^{24\times 6}$$

$${\mathbf{R}}_y^{3\mathrm{D}}=\left[0_8{1}_8{0}_8-{\mathit{z}}_{\mathrm{n}}{0}_8-{\mathit{x}}_{\mathrm{n}}\right]\in {\mathbb{R}}^{24\times 6}$$

$${\mathbf{R}}_z^{3\mathrm{D}}=\left[0_8{0}_8{1}_8{\mathit{y}}_{\mathrm{n}}-{\mathit{x}}_{\mathrm{n}}{0}_8\right]\in {\mathbb{R}}^{24\times 6}$$

$${\mathit{x}}_{\mathrm{n}}={[\mathrm{0,1},x_3,x_4, 0, 1,x_3,x_4]}^{\mathrm{T}}$$

$${\mathit{y}}_{\mathrm{n}}={[\mathrm{0,0},y_3,y_4, 0, 0,y_3,y_4]}^{\mathrm{T}}$$

$${\mathit{z}}_{\mathrm{n}}={[-{\displaystyle \frac{t}{2}},-{\displaystyle \frac{t}{2}},-{\displaystyle \frac{t}{2 }},-{\displaystyle \frac{t}{2}}, {\displaystyle \frac{t}{2}}, {\displaystyle \frac{t}{2}},{\displaystyle \frac{t}{2}},{\displaystyle \frac{t}{2}}]}^{\mathrm{T}}$$

The first three columns of ${\mathbf{R}}^{3\mathrm{D}}$ represent rigid-body translations along the x-, y-, and z-directions, respectively, while the last three columns represent rigid-body rotations about the x-, y-, and z-axes, respectively.

The rigid-body zero-energy constraints can then be expressed as by Eq. (4), this equation yields 144 homogeneous linear constraint equations with respect to the entries of the element stiffness matrix ${\mathbf{K}}^{3\mathrm{D}}$, although these equations are not necessarily linearly independent.

Constant Generalized-Strain Consistence Constraints

After excluding rigid-body translations and rotations, in terms of the three-independent generalized-strain components ${\varphi }_1$, ${\varphi }_2$, and ${\varphi }_3$ as

$$\begin{array}{c}\begin{array}{c}u\left(x,y,z\right)=-z\left({\varphi }_{1}x+{\varphi }_{3}y\right)\\ v\left(x,y,z\right)=-z\left({\varphi }_{3}x+{\varphi }_{2}y\right)\end{array}\\ w\left(x,y,z\right)=-{\displaystyle \frac{\left({\varphi }_1{x}^2+{\varphi }_2{y}^2+{\varphi }_{3}xy\right)}{2}}\end{array},$$

(40)

This representation corresponds to the constant generalized-strain vector

$$\mathit{\varphi }={\left[\begin{array}{ccc}{\varphi }_1& {\varphi }_2& {\varphi }_3\end{array}\right]}^{\mathbf{T}},$$

(41)

where ${\varphi }_1$, ${\varphi }_2$, and ${\varphi }_3$ are constant generalized-strain modes associated with bending in the x-direction, bending in the y-direction, and twisting in the xy-direction, respectively.

The nodal displacement modes associated with ${\varphi }_1=1$, ${\varphi }_2=1$ and ${\varphi }_3=1$ are assembled into:

$${\mathbf{S}}^{3\mathrm{D}\mathrm{g}}({\mathit{x}}^{3\mathrm{D}})=\left[\begin{array}{c}\begin{array}{c}{\mathbf{S}}_x^{3\mathrm{D}\mathrm{g}}\\ {\mathbf{S}}_y^{3\mathrm{D}\mathrm{g}}\end{array}\\ {\mathbf{S}}_z^{3\mathrm{D}\mathrm{g}}\end{array}\right]\in {\mathbb{R}}^{24\times 3}$$

(42)

where

$${\mathbf{S}}_{\mathit{x}}^{3\mathrm{D}\mathrm{g}}=t\mathbf{H}{\mathbf{G}}_x , {\mathbf{S}}_{\mathit{y}}^{3\mathrm{D}\mathrm{g}}=t{\mathbf{H}\mathbf{G}}_y , {\mathbf{S}}_{\mathit{z}}^{3\mathrm{D}\mathrm{g}}={\displaystyle \frac{1}{2}}t\mathbf{H}{\mathbf{G}}_z ,$$

here

$$\mathbf{H}=\left[\begin{array}{c}-1\\ 1\end{array}\right] \otimes {\mathbf{I}}_4$$

and

$${\mathbf{G}}_x=\left[\begin{array}{ccc}0& 0& 0\\ 1& 0& 0\\ y_3& 0& x_3\\ y_4& 0& x_4\end{array}\right] , {\mathbf{G}}_y=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& y_3& x_3\\ 0& y_4& x_4\end{array}\right] , {\mathbf{G}}_z=\left[\begin{array}{ccc}0& 0& 0\\ 1& 0& 0\\ x_3^2& y_3^2& x_3{y}_3\\ x_4^2& y_4^2& x_4{y}_4\end{array}\right]$$

For each unit constant generalized-strain mode, the corresponding equivalent nodal force vector is obtained from the theoretical generalized stress resultants. Following the analytical boundary-integration form, the equivalent nodal force matrix associated with the three generalized-strain modes is written as

$${\mathbf{F}}_S^{3\mathrm{D}\mathrm{g}}\left({\mathit{x}}^{3\mathrm{D}}\right)={\displaystyle \frac{D_0}{4}}\left({\mathbf{I}}_3\otimes \mathbf{L}\right)\overline{\mathbf{G}}\mathit{\varphi }$$

(43)

where

$$D_0={\displaystyle \frac{E\left(2t^3\right)}{12\left(1-{\mu }^2\right)}}, \mathbf{L}=\left[\begin{array}{c}1\\ -1\\ -1\\ 1\end{array}\right]\otimes {\mathbf{I}}_2$$

$$\overline{\mathbf{G}}({\mathit{x}}^{3\mathrm{D}})=\left[\begin{array}{ccc}{y}_4& y_4& (1-\mu )(1-x_4)/2\\ y_3& y_3& (1-\mu )x_3/2\\ (1-x_4)\mu & 1-x_4& (1-\mu )y_4/2\\ x_3\mu & x_3& -(1-\mu )y_3/2\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$$

The constant generalized-strain consistency constraints can then be expressed by Eq. (8), which provides 72 non-homogeneous linear equality constraints with respect to the entries of the element stiffness matrix ${\mathbf{K}}^{3\mathrm{D}}$, although these equations are not necessarily linearly independent.

Membrane-Mode Constant-Strain Constraints

After excluding the rigid-body translations and the in-plane rigid-body rotation, the non-rigid displacement field can be expressed in terms of three independent membrane strain components as:

$$\begin{array}{c}\begin{array}{c}u\left(x,y,z\right)={\epsilon }_{xx}^{m}x+{\displaystyle \frac{1}{2}}{\gamma }_{xy}^{m}y\\ v\left(x,y,z\right)={\epsilon }_{yy}^{m}y+{\displaystyle \frac{1}{2}}{\gamma }_{xy}^{m}x\end{array}\\ w\left(x,y,z\right)=-{\displaystyle \frac{{\mu }^{\mathrm{*}}}{1-{\mu }^{\mathrm{*}}}}\left({\epsilon }_{xx}^m+{\epsilon }_{yy}^m\right)z+C_0\end{array}$$

(44)

This representation (46) corresponds to the membrane constant-strain vector

$${\mathsf{\epsilon }}^m={\left[{\epsilon }_{xx}^m,{\epsilon }_{yy}^m,{\gamma }_{xy}^m\right]}^{\mathrm{T}}$$

(45)

where the incompressibility condition is imposed in the z-direction, i.e., (${\mu }^{*}=0.5$), and ($C_0$) is a constant.

The nodal displacement modes associated with ${\epsilon }_{xx}^m=1$,${\epsilon }_{yy}^m=1$ and ${\gamma }_{xy}^m=1$ can be assembled into:

$${\mathbf{S}}^{3\mathrm{D}\mathrm{m}}({\mathit{x}}^{3\mathrm{D}})=\left[\begin{array}{c}\begin{array}{c}{\mathbf{S}}_x^{3\mathrm{D}\mathrm{m}}\\ {\mathbf{S}}_y^{3\mathrm{D}\mathrm{m}}\end{array}\\ {\mathbf{S}}_z^{3\mathrm{D}\mathrm{m}}\end{array}\right]\in {\mathbb{R}}^{24\times 3}$$

(46)

where

$${\mathbf{S}}_x^{3\mathrm{D}\mathrm{m}}=\left[\begin{array}{cc}{1}_2\otimes \mathbf{N}& 0\end{array}\right],{\mathbf{S}}_y^{3\mathrm{D}\mathrm{m}}=\left[\begin{array}{cc}0& 1_2\otimes \mathbf{N}\end{array}\right],{\mathbf{S}}_z^{3\mathrm{D}\mathrm{m}}=-{\displaystyle \frac{1}{2}}t\mathbf{H}\left[\begin{array}{cccccc}1& 0& 0& 0& 1& 0\end{array}\right]$$

and

$${\mathit{\varphi }}^{3\mathrm{D}\mathrm{m}}={\left[\begin{array}{cccccc}{\alpha }_1& {\alpha }_2& {\alpha }_3& {\beta }_1& {\beta }_2& {\beta }_3\end{array}\right]}^{\mathrm{T}}$$

$$\mathbf{H}=\left[\begin{array}{c}-1\\ 1\end{array}\right] \otimes {\mathbf{I}}_4$$

$$\mathbf{N}=\left[\begin{array}{ccc}0& 0& 1\\ 1& 0& 1\\ x_3& y_3& 1\\ x_4& y_4& 1\end{array}\right]$$

For each membrane-mode constant-strain case, the corresponding equivalent nodal force vector is obtained from the theoretical membrane stress resultants through boundary traction integration. Following the analytical boundary-integration form, the associated equivalent nodal force matrix is written as:

$${\mathbf{F}}_S^{3\mathrm{D}\mathrm{m}}\left({\mathit{x}}^{3\mathrm{D}}\right)=D_0\left({\mathbf{I}}_4\otimes \mathbf{L}\right)\mathbf{G}\mathbf{T}$$

(47)

where

$$\mathit{G}({\mathit{x}}^{3\mathrm{D}})=\left[\begin{array}{ccc}{y}_{4}t& 0& {\displaystyle \frac{\left(x_4-1\right)t}{2}}\\ -y_{3}t& 0& {\displaystyle \frac{x_{3}t}{2}}\\ 0& x_4-1& {\displaystyle \frac{y_{4}t}{2}}\\ 0& x_3& -{\displaystyle \frac{y_{3}t}{2}}\end{array}\right]$$

$$\mathbf{T}=\left[\begin{array}{cccc}1& 0& 0& \mu \\ \mu & 0& 0& 1\\ 0& 1-\mu & 1-\mu & 0\end{array}\right]$$

The membrane-mode constant-strain constraints for the solid-shell element stiffness matrix are therefore expressed by Eq. (8), which provides 72 non-homogeneous linear equality constraints with respect to the entries of the element stiffness matrix ${\mathbf{K}}^{3\mathrm{D}}$, although these equations are not necessarily linearly independent.

By jointly enforcing the rigid-body zero-energy constraints, the constant generalized-strain constraints, and the membrane-mode constant-strain constraints, the complete linear constraint system for the three-dimensional solid-shell element can be written as

$${\mathbf{K}}^{3\mathrm{D}}{\mathbf{C}}^{3\mathrm{D}}\left({\mathit{x}}^{3\mathrm{D}}\right)={\mathbf{D}}^{3\mathrm{D}}\left({\mathit{x}}^{3\mathrm{D}}\right)$$

(48)

$${\mathbf{C}}^{3\mathrm{D}}\left({\mathit{x}}^{3\mathrm{D}}\right)=\left[\begin{array}{ccc}{\mathbf{R}}^{3\mathrm{D}}& {\mathbf{S}}^{3\mathrm{D}\mathrm{g}}& {\mathbf{S}}^{3\mathrm{D}\mathrm{m}}\end{array}\right]\in {\mathbb{R}}^{24\times 12}$$

$${\mathbf{D}}^{3\mathrm{D}}\left({\mathit{x}}^{3\mathrm{D}}\right)=\left[\begin{array}{ccc}0& {\mathbf{F}}_S^{3\mathrm{D}\mathrm{g}}& {\mathbf{F}}_S^{3\mathrm{D}\mathrm{m}}\end{array}\right]\in {\mathbb{R}}^{24\times 12}$$

3.2.3. Identification of Independent Stiffness Parameters and Construction of the Reconstruction Operator for the Solid-Shell Element

For the eight-node solid-shell element, the same vectorized linear constraint formulation and symbolic elimination procedure described in Section 3.1.3 are adopted. Under the canonical configuration defined in Section 3.2.1, the physics-informed hard constraints lead to a linear algebraic system for the lower-triangular entries of the symmetric element stiffness matrix. The dimension of the resulting admissible solution space is $n_p^{3\mathrm{D}}=78$. Therefore, 78 independent stiffness parameters are required to parameterize the element stiffness matrix satisfying the prescribed linear hard constraints for the present solid-shell element.

Following the fixed symbolic elimination strategy used for the two-dimensional quadrilateral element, a stable and consistent set of free stiffness entries is selected. The independent stiffness parameter vector of the solid-shell element is written as

$${\mathit{p}}^{3\mathrm{D}}=\left[\begin{array}{cccccccccc}{\widehat{K}}_{10,10}& {\widehat{K}}_{10,11}& \dots & {\widehat{K}}_{11,11}& {\widehat{K}}_{11,15}& \dots & {\widehat{K}}_{11,24}& {\widehat{K}}_{15,15}& \dots & {\widehat{K}}_{24,24}\end{array}\right]$$

(49)

The locations of the independent stiffness entries are shown in Figure 6

Based on the symbolic solution of the constraint system, the full element stiffness matrix can be reconstructed from the independent parameter vector through the element-specific hard-constrained reconstruction operator:

$${\mathbf{K}}^{3\mathrm{D}}={\mathcal{Q}}_{3\mathrm{D}}\left({\mathit{p}}^{3\mathrm{D}};{\mathit{x}}^{3\mathrm{D}}\right)$$

(50)

Equivalently, the element stiffness matrix can be expressed as

$${\mathbf{K}}^{3\mathrm{D}}=K^{\left({\mathit{x}}^{3\mathrm{D}}\right)}+K^{\left(\mathrm{10,10}\right)}{\widehat{K}}_{\mathrm{10,10}}+K^{\left(\mathrm{10,10}\right)}{\widehat{K}}_{\mathrm{10,11}}+\dots +K^{\left(\mathrm{15,15}\right)}{\widehat{K}}_{\mathrm{15,15}}+\dots$$

(51)

here,${\widehat{K}}_{10,10},...$ are independent parameters, $K^{\left({\mathit{x}}^{3\mathrm{D}}\right)}$denotes the term from the non-homogeneous part in Eq. (51) and is irrelevant with independent parameters. The terms $K_{i,j}^{\left(\mathrm{10,10}\right)},...$ are coefficient functions associated with the 78 independent stiffness parameters, respectively. To illustrate the resulting reconstruction relation, representative stiffness entries are listed as follows:

$$K_{12,23}^{3\mathrm{D}}={\displaystyle \frac{-2t{x}_3+2t{x}_4}{x_3{y}_4-y_3{x}_4}}{\widehat{K}}_{15,24}+{\displaystyle \frac{-x_3{y}_4+y_3{x}_4-y_3+y_4}{x_3{y}_4-y_3{x}_4}}{\widehat{K}}_{17,24}+\dots$$

This reduction from 300 independent entries of a general symmetric 24×24 matrix to only 78 free stiffness parameters quantifies how strongly the prescribed mechanical constraints restrict the admissible stiffness description.

4. Data Generation and Neural-Network Training

4.1. Data Generation Based on Canonical Configuration

The training datasets are generated according to the canonical configuration established in Section 3.

For the four-node plane-stress quadrilateral element, geometric samples are generated under the two-dimensional canonical configuration defined in Section 3.1.1. The first edge is fixed as a unit edge aligned with the global x-axis, while the remaining two nodal coordinates $\left(x_3,y_3\right)$ and $\left(x_4,y_4\right)$ are randomly sampled within a prescribed admissible region.

Specifically, to ensure stable finite element computations and reliable sample quality, each candidate quadrilateral geometry is required to satisfy the following geometric quality constraints. First, the quadrilateral must be convex. Second, all interior angles must satisfy ${\theta }_i>{45}^{°}, i=1,\dots ,4$. Finally, to avoid excessive edge-length distortion, the edge-length ratio of each candidate quadrilateral is bounded by ${\displaystyle \frac{\mathrm{m}\mathrm{a}\mathrm{x}\left(l_j\right)}{\mathrm{m}\mathrm{i}\mathrm{n}\left(l_j\right)}}\le 4$, where $l_j$ denotes the length of the $j$-th edge. The locations of node 3 and node 4 are depicted in Figure 7. Rectangular elements are regarded as a special case.

For the eight-node solid-shell element, the midsurface geometry is generated using the same strategy as the four-node plane-stress quadrilateral element. Once a valid midsurface configuration is obtained, it is extruded along the midsurface normal to construct the solid-shell element. The thickness parameter $t$ is uniformly sampled as $t~U\left(\mathrm{0.1,1}\right)$. This unified sampling strategy ensures parameterization consistency between the two-dimensional quadrilateral-element dataset and the solid-shell element dataset.

All samples are assumed to follow an isotropic linear-elastic constitutive model. The Young's modulus is normalized to 1, while Poisson's ratio is retained as the material input parameter. Specifically, Poisson's ratio is uniformly sampled as $\mu ~U\left(\mathrm{0.2,0.45}\right)$.

Once the samples associated with the geometric and material parameters are defined according to the above sampling strategy, the corresponding element stiffness matrices are generated using independently developed finite-element programs for the plane-stress elements and commercial finite-element software for the solid-shell elements. Finally, the training dataset $\mathcal{D}={\left\{\left({\mathit{x}}^{\left(s\right)},{\mathbf{K}}_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\left(s\right)}\right)\right\}}_{s=1}^N$ comprising $N$ samples is established.

The related element types for numerical examples in this paper are summarized in Table 1. Here, Q4, Q4R, Q4I, and SS8 denote the conventional finite elements used as reference elements, whereas HCML-Q4, HCML-Q4R, HCML-Q4I, and HCML-SS8 denote the corresponding learned data-driven elements by Hard-Constrained Manifold Learning (HCML). The purpose of learning Q4, Q4R, and Q4I separately is not merely to reproduce different stiffness matrices, but to examine whether numerical characteristics associated with different element formulations can be inherited through the same HCML construction framework.

4.2. Neural-Network Architecture and Training

For the four-node plane-stress quadrilateral element and the three-dimensional eight-node solid-shell element, two fully connected feedforward neural-networks are constructed to accomplish the Hard-Constrained Manifold Learning (HCML) framework. Instead of predicting the complete element stiffness matrix, each network predicts only the independent stiffness parameter vector, and the complete element stiffness matrix is then obtained by the element stiffness matrix reconstruction operator in Section 3.

Because the input dimensions and the number of independent stiffness parameters differ between the two element types, separate network architectures are used. The architectures adopted in this work are summarized in Table 2.

For a training set containing $N$ samples, the mean squared error loss is defined in the independent-parameter space as

$${\mathcal{L}}_{\mathrm{M}\mathrm{S}\mathrm{E}}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}{‖{\widehat{\mathit{p}}}_i-{\mathit{p}}_{i,\mathrm{r}\mathrm{e}\mathrm{f}}‖}_2^2$$

(52)

where ${\widehat{\mathit{p}}}_i$ is the predicted independent stiffness parameter vector and ${\mathit{p}}_{i,\mathrm{r}\mathrm{e}\mathrm{f}}$ is the corresponding reference parameter vector extracted from the reference element stiffness matrix. Since the full element stiffness matrix is reconstructed through the element stiffness matrix reconstruction operator Eq. (19), the prescribed linear physics-informed hard constraints are satisfied by construction and are not imposed through additional penalty terms in the loss function. Further details of the training settings, together with the training and validation loss curves, are provided in Appendix C.

To further examine the predictive capability of the proposed model under continuous variations in element geometry, a regular unit-square element is taken as the baseline configuration. Only the vertical coordinate of node 4, $y_4$, is varied continuously within the range of 0.50-1.50. For the three types of plane-stress elements, the selected independent stiffness entries $K_{7,7},K_{7,8},K_{8,8}$​ are computed using both the trained neural-networks and reference elements. Their variations with respect to $y_4$ are then compared. The results show that the neural-network predictions almost overlap with the results of reference elements over the entire examined range, indicating that the proposed model can accurately capture the nonlinear variation of stiffness components induced by geometric distortion for different element types.

Figure 8. Comparison of the neural-network predictions and results of reference elements for the stiffness components of Q4, Q4R, and Q4I elements as $y_4$ aries.

Figure 8. Comparison of the neural-network predictions and results of reference elements for the stiffness components of Q4, Q4R, and Q4I elements as $y_4$ aries.

5. Implementation of the Direct Data-Driven Elements and Complete Flowchart

Figure 9 presents the complete computational workflow of the HCML-based direct data-driven elements, covering element input, neural-network prediction, element stiffness matrix reconstruction, and UEL-based finite element analysis. The element geometry and material parameters are first assembled into the input vector $\mathit{x}=\left(\mathit{g},\mathit{m}\right)$, from which the pretrained neural-network predicts the independent stiffness parameter vector $\mathit{p}$. The full element stiffness matrix is then reconstructed from $\mathit{p}$ using the reconstruction operator derived from the linear system of physics-informed hard constraints. The learned elements are implemented through the Abaqus user element (UEL) interface, which enables user-defined element formulations to provide the element stiffness matrix and residual vector to the finite element solver through a Fortran subroutine. Before the finite element analysis, the parameters of the pretrained neural network are fixed and embedded into the UEL subroutine. During each UEL call, the subroutine successively performs element-input evaluation, neural-network forward propagation, hard-constrained stiffness-matrix reconstruction, stiffness scaling, and local-to-global coordinate transformation. The resulting element stiffness matrix and residual vector are returned to Abaqus for global assembly and equilibrium solution under the prescribed loads and boundary conditions. Finally, auxiliary skin elements, which share nodes with the UEL mesh and serve only as post-processing carriers without contributing to the structural stiffness, are employed when required to visualize the resulting displacement and stress fields.

6. Numerical Examples

This section presents numerical examples to test the performance of the HCML-FEs. All the learned element types have been summarized in Table 1. Section 6.1 presents fundamental numerical tests, including the rigid-body zero-energy mode test and the first-order patch tests for both plane-stress and solid-shell elements. Section 6.2 evaluates the performance of the HCML-Q4, HCML-Q4R, and HCML-Q4I elements with respect to mesh-distortion and mesh-dependence tests, as well as representative two-dimensional structural benchmarks, including Cook's skew beam and the strap plate problem. Section 6.3 assesses the HCML-SS8 element using three-dimensional solid-shell benchmarks, namely the Barrel vault roof problem, the Raasch Challenge problem, the hyperboloid-shell, and the thin-walled pipe elbow. Unless otherwise specified, all numerical examples in this section are performed using a consistent unit system.

6.1. Fundamental Numerical Tests

6.1.1. Rigid-Body Zero-Energy Mode Test

The rigid-body zero-energy mode test is a standard element-level verification of the null-space property of an element stiffness matrix.

For the two-dimensional elements, HCML-Q4, HCML-Q4R, and HCML-Q4I are tested using three quadrilateral geometries, namely a rectangular quadrilateral, an arbitrary quadrilateral, and a distorted quadrilateral, as shown in Figure 10(a)-11(c). For the three-dimensional case, HCML-SS8 is tested using a representative solid-shell geometry as shown in Figure 10(d).

For the HCML-FEs, Young's modulus and Poisson's ratio are set to 2E+04 and 0.3, respectively. For each test case, the complete element stiffness matrix is generated using the proposed physics-informed hard-constrained data-driven finite element formulation. The eigenvalues of the stiffness matrix are then computed and sorted in ascending order, and the number of near-zero eigenvalues is compared with the theoretical number of rigid-body modes, as shown in Table 3. In addition, the near-zero eigenvalue with maximum absolute value and the first nonzero eigenvalue are also included in Table 3.

From the results, the number of detected near-zero eigenvalues is consistent with the theoretical number of rigid-body modes for all tested elements and geometric configurations. The HCML-Q4, HCML-Q4R, and HCML-Q4I elements exhibit three near-zero eigenvalues under the rectangular, arbitrary quadrilateral, and distorted quadrilateral configurations, corresponding to two rigid-body translational modes and one out-of-plane rigid-body rotational mode. In contrast, the HCML-SS8 element exhibits six near-zero eigenvalues, corresponding to the six rigid-body modes in three-dimensional space. For all test cases, the near-zero eigenvalues are several orders of magnitude smaller than the first nonzero eigenvalue, indicating a clear separation between rigid-body modes and deformational modes. These near-zero eigenvalues can be attributed to numerical round-off errors, and no additional spurious zero-energy modes are observed. These results demonstrate that the proposed HCML-FEs successfully pass the rigid-body zero-energy mode test.

As shown in the last row of the table, the soft constraints method may yield a relatively small element stiffness matrix prediction error. However, the predicted element stiffness matrix does not retain the zero eigenvalues associated with rigid-body motions and therefore fails the rigid-body zero-energy mode test. The soft-constraint baseline illustrates a central distinction between stiffness-matrix prediction accuracy and finite-element admissibility. Although the predicted stiffness entries may remain close to the reference values, the rigid-body null space is not preserved. Thus, small regression errors alone are insufficient to define a mechanically valid element, whereas the HCML reconstruction preserves the prescribed null-space constraint by construction. Similarly, a finite element constructed from a soft-constraint baseline cannot precisely pass the patch test discussed in the next subsection.

6.1.2. Patch Test Under Constant Strain and Constant Generalized-Strain

The plane-stress patch test is first performed using arbitrary quadrilateral elements to evaluate the HCML-Q4, HCML-Q4R, and HCML-Q4I. The MacNeal-Harder patch geometry (MacNeal and Harder, 1985) is adopted, as shown in Figure 11. The boundary nodes are prescribed according to a linear displacement field, while the interior nodal displacements are obtained from the finite element solution. Young's modulus and Poisson's ratio are set to 2E+04 and 0.3. For the arbitrary quadrilateral patch, the displacement components in the x- and y-directions are prescribed at the boundary nodes as:

$$u={10}^{-3}\left(x+{\displaystyle \frac{y}{2}}\right),v={10}^{-3}\left(y+{\displaystyle \frac{x}{2}}\right)$$

(53)

The analytical displacements prescribed by the linear displacement field and displacements computed with the HCML-FEs at the interior nodes are listed in Table 4. The displacement components are scaled by E-04.

The results show that the computed interior nodal displacements agree well with the analytical values, with relative differences on the order of E-08. This confirms that the arbitrary quadrilateral HCML-FE patch can reproduce the prescribed linear displacement field. Therefore, the arbitrary quadrilateral HCML-FEs pass the first-order plane-stress patch test.

The Constant generalized-strain patch test is then conducted for the three-dimensional HCML-SS8 element. As shown in Figure 12, a rectangular plate with uniform thickness is used as the test patch that is discretized by nine solid-shell elements with arbitrary quadrilateral midsurfaces. Young's modulus and Poisson's ratio are set to 2E+04 and 0.3. The external nodes are prescribed by a displacement field that produces a constant generalized-strain state, while the interior nodal displacements are solved from the finite element equilibrium equations with the proposed HCML-SS8. The imposed displacement field is given by:

$$u={10}^{-3}\left(x+{\displaystyle \frac{y}{2}}\right)z,v={10}^{-3}\left(y+{\displaystyle \frac{x}{2}}\right)z,w=-{10}^{-3}{\displaystyle \frac{x^2+xy+y^2}{2}}.$$

(54)

The analytical displacements prescribed by Eq. (54) and the calculated displacements by the proposed HCML-SS8 at the interior nodes are summarized in Table 5. The displacement components are scaled by E-06.

The calculated displacements of the interior nodes are in excellent agreement with the analytical values. This result shows that the HCML-SS8 element can successfully pass the constant generalized-strain patch test and reproduce the constant bending/torsion response required for solid-shell deformation.

6.2. Structural Numerical Examples Using the HCML-Q4, HCML-Q4R and HCML-Q4I Elements

6.2.1. Mesh-Distortion and Mesh-Dependence Tests

MacNeal's cantilever beam benchmark (MacNeal and Harder, 1985) is used to evaluate mesh-distortion sensitivity and the performance of HCML-FEs with arbitrary quadrilateral elements. As shown in Figure 13, the length, height, and thickness of the beam are $L=6$, $H=0.4,$ and $t=1$, respectively. Young's modulus and Poisson's ratio are set to 2E+04 and 0.3. The left end is clamped, and two loading cases are applied at the free end: (I) a tip shear force of $P=1$ and (II) a tip bending moment of $M=-0.2$. Three mesh patterns are examined: a regular mesh, a parallelogram mesh, and a trapezoidal mesh, denoted as Mesh I, Mesh II, and Mesh III, respectively. In this benchmark, we mainly compare the vertical displacement $v_{\mathrm{A}}$ at point A for two load cases and three mesh patterns. Here, the reference vertical displacements are $v_{\mathrm{A}}=-0.02$ for load $P$ and $v_{\mathrm{A}}=-0.001$ for load$M$. All the HCML-Q4, HCML-Q4R, and HCML-Q4I elements are examined and compared with their corresponding reference elements, respectively, as shown in Table 6.

Table 6 summarizes the vertical displacement $v_A$ obtained by reference elements and the HCML-FEs under different mesh levels. The reported displacement values are scaled by E-02. The benchmark is bending-dominated, and only one element is used through the thickness direction; therefore, it is sensitive to both the element formulation and mesh distortion. Except for Q4R and HCML-Q4R, the computational results of the HCML-FEs for the other two element types are close to those of the corresponding reference elements.

The results show that the present benchmark employs arbitrary quadrilateral meshes, and mesh distortion has a pronounced influence on the numerical response of low-order quadrilateral elements. Except for Q4R and HCML-Q4R, the HCML-FEs show good agreement with the corresponding reference elements and are able to reproduce the mesh-dependence characteristics of different reference formulations. Both Q4 and HCML-Q4 exhibit an overly stiff response, and mesh distortion further reduces the displacement magnitude. Q4I and HCML-Q4I provide relatively high accuracy on the regular mesh, but over-stiffening is still observed under distorted meshes, indicating that the incompatible-mode formulation improves bending performance but its effectiveness is sensitive to mesh distortion, which is like the nonconforming elements. In contrast, Q4R and HCML-Q4R exhibit evident hourglass-dominated responses in coarse distorted-mesh cases. The Q4R results also clarify the role of the training data in the proposed framework. The hard constraints define the mechanically admissible structure imposed on the learned stiffness matrix, but they do not remove numerical characteristics from the reference data that are not excluded by these constraints. Consequently, undesirable behavior contained in the reference element, such as the pronounced coarse-mesh sensitivity of the reduced-integration formulation, may also be inherited. In this sense, the physical constraints determine admissibility, whereas the reference data determine the remaining numerical character of the learned element.

6.2.2. Cook's Skew Beam

Cook et al.’s (2001) skew beam benchmark is used to evaluate the performance under combined geometric distortion and mesh refinement of the proposed HCML-FEs, as shown in Figure 14. Young's modulus and Poisson's ratio are set to 2E+04 and 0.3. The left edge is fully clamped, and a vertical upward concentrated load $P=1$ is applied at point A on the right end. The geometric inclination angle is denoted by $\alpha$. When $\alpha ={42.51}^{°}$, the vertical displacement at point A is considered, and its theoretical value is $v_{\mathrm{A}}^{ref}$ = 1.20E-04. For the Cook’s skew beam benchmark, the same number of mesh seeds, denoted by $N$, is assigned along each of the four edges. The domain is then discretized into arbitrary quadrilateral elements. It is evident that different values of the inclination angle $\alpha$ lead to different levels of element and mesh distortion.

When $\alpha ={42.51}^{°}$, Table 7 summarizes the maximum vertical displacement $v_A$​ obtained by reference elements and the HCML-Fes under different mesh levels. The reported displacement values are scaled by E-04. Figure 15 shows the displacement and von Mises stress contours of Q4I and HCML-Q4I for $N=4$. The results show that, in the Cook’s skew beam problem, the HCML-FEs can reproduce the typical convergence behavior of the corresponding reference elements. The Q4 and HCML-Q4 elements underestimate the displacement on coarse meshes, reflecting the relatively stiff response of bilinear quadrilateral elements in bending-dominated skew geometries. The Q4R and HCML-Q4R elements exhibit stronger coarse-mesh sensitivity, but their responses gradually stabilize with mesh refinement. In contrast, owing to their improved bending representation, the Q4I and HCML-Q4I elements remain closer to the benchmark reference solution over a wider range of mesh densities. Meanwhile, the displacement and von Mises stress contours of Q4I and HCML-Q4I shown in Figure 15 for $N=4$ exhibit generally consistent field distributions.

To further demonstrate the generalization capability of the HCML framework, Figure 16 compares the vertical displacement at point A and its relative error between Q4I and HCML-Q4I under different geometric inclination angles. Figure 16 shows that, as the geometric inclination angle $\alpha$ increases, some element interior angles gradually fall outside the coverage of the training dataset, with the minimum interior angle becoming smaller than ${45}^{°}$. For moderately out-of-domain geometries, the proposed HCML-Q4I still maintains acceptable prediction accuracy. However, the displacement error increases markedly as $\alpha$ further increases, indicating that the learned element possesses a certain degree of geometric generalization and extrapolation capability, although this capability is limited to a finite range of geometric distortion. By expanding the training domain, specifically by reducing the lower bound of the minimum interior angle to ${35}^{°}$, the relative errors between HCML-Q4I and Q4I are kept within approximately 1.1%, even when $\alpha$ deviates from ${45}^{°}$. Overall, the HCML-FEs not only accurately reproduce the responses of the corresponding reference elements and exhibit a reasonable convergence trend toward the benchmark reference solution, but also demonstrate a certain degree of geometric generalization and extrapolation capability. The results distinguish two forms of generalization. Within the sampled geometric domain, the learned element reproduces the reference response over continuously varying element shapes. Outside this domain, the mechanical constraints remain exactly satisfied, whereas the accuracy of the unconstrained stiffness information gradually deteriorates as the geometry moves farther from the training range.

6.2.3. Strap Plate Problem

The strap plate benchmark is used to assess the applicability of the proposed HCML-FEs to perforated structures with practical geometric complexity and pronounced local stress concentrations. As shown in Figure 17, the plate contains three circular holes. All nodes on the boundaries of the two left holes are fully constrained, and a concentrated load $q=1$ is applied at point A, located at the top of the right hole, in the negative y-direction. Three mesh densities labeled as Mesh I, Mesh II, and Mesh III in Figure 17 are adopted to examine the mesh-dependence behavior of the numerical response. The considerable variations in element shape and size within this mesh provide a stringent test of the geometric generalization capability of the HCML-FEs. Young's modulus and Poisson's ratio are set to 2E+04 and 0.3. The reference value of the vertical nodal displacement at the loaded point A is $v_{\mathrm{A}}=$ -5.374 E-04, which is obtained using a highly refined mesh.

Table 8 compares the maximum vertical nodal displacement at the loaded point A obtained by reference elements and the HCML-FEs under different mesh densities. The reported displacement values are scaled by E-04. Figure 18 presents the displacement and von Mises stress contours of Q4I and HCML-Q4I for Mesh I. The results show that, under different mesh densities, the learned elements remain in good agreement with their corresponding reference elements and exhibit consistent convergence trends. Figure 18 further indicates that HCML-Q4I and Q4I show highly consistent full-field displacement distributions, von Mises stress distributions, and stress-concentration regions near the loaded area. Overall, the strap plate benchmark demonstrates that the proposed HCML-FEs can be applied to engineering structures with complex geometric features and can accurately reproduce the mechanical responses of the corresponding reference elements.

6.3. Complex Numerical Examples for the HCML-SS8 Element

6.3.1. Barrel Vault Roof

The Barrel vault roof benchmark is used to evaluate the performance of the proposed learned solid-shell element in a cylindrical shell problem involving coupled membrane and bending deformation. As shown in Figure 19, the roof has a length $L=50$, radius $R=25$, and thickness $t=0.25$. Young's modulus and Poisson's ratio are set to 2E+11 and 0.3. The four boundary edges are fully clamped, and uniformly distributed loading is represented by equivalent nodal forces applied in the global z-direction. The model is discretized using 25×16×1 and 50×32×2 meshes, where the three values denote the numbers of elements in the longitudinal, circumferential, and thickness directions, respectively. For the initial 25×16×1 mesh, a unit vertical nodal load is applied to each node in the z-direction. After mesh refinement, the nodal force is adjusted according to the number of loaded nodes so that the total applied load remains unchanged.

Table 9 presents the relative errors in the maximum z-direction displacement between HCML-SS8 and SS8 for the two meshes. The displacement values are scaled by E-10. Figure 20 and Figure 21 compare the z-direction displacement and von Mises stress contours obtained using SS8 and HCML-SS8 on the 25×16×1 mesh.

The results demonstrate that the HCML-SS8 solution approaches the SS8 reference solution as the mesh is refined. The two elements exhibit consistent overall deformation modes, displacement gradients, and locations of the maximum-displacement region. The von Mises stress distributions on the top surface, midsurface, and bottom surface also show good agreement in their spatial patterns and principal stress-concentration regions. Together with the reduction in the displacement difference after mesh refinement, these results indicate that HCML-SS8 reproduces the dominant membrane–bending coupled response of the reference SS8 element within the investigated meshes.

6.3.2. Raasch Challenge

The Raasch Challenge is a classical shell benchmark commonly used to assess the capability of shell elements in representing coupled shear, bending, and twisting deformation in curved thin-strip structures. As shown in Figure 22, one end of the structure is rigidly clamped, while a nodal load with a magnitude of 1 is applied to each node at the free end in the positive z-direction. The material is assumed to be linearly elastic. Young's modulus and Poisson's ratio are set to 2E+11 and 0.3, respectively. The model is discretized using a structured SS8 mesh of 20×144×2, where 20, 144, and 2 denote the numbers of elements along the width, curved longitudinal, and thickness directions, respectively.

Figure 23 compares the z-direction displacement contours and the bottom-surface von Mises stress contours obtained using the SS8 element and the HCML-SS8 element. The reported displacement values are scaled by E-06. The two elements exhibit good agreement in terms of the overall deformation pattern, displacement distribution, and stress-field distribution. The relative error in the maximum nodal displacement in the z-direction is 2.83%, indicating that HCML-SS8 can accurately reproduce the displacement response of SS8.

Because the Raasch Challenge involves strongly coupled in-plane bending, transverse shear, and twisting deformations, this agreement demonstrates that HCML-SS8 effectively reproduces the capability of SS8 to represent the complex coupled response under the discretization conditions considered. Furthermore, the displacement and stress fields vary smoothly across the junction of the two circular segments, without evident local numerical oscillations, nonphysical discontinuities, or spurious flexibility. Overall, this benchmark confirms the applicability and reliability of HCML-SS8 for curved shell structures involving coupled deformation modes.

6.3.3. Hyperboloid-Shell

The model is a ${180}^{°}$ sector of a shell of revolution generated by rotating an offset semi-elliptical meridian about the y-axis. As shown in Figure 24, the shell has a waist diameter $D_{\mathrm{w}}=90$, an end diameter $D_{\mathrm{e}}=140$, and an axial length $L_{\mathrm{a}}=90$, and a uniform thickness $t=1$. The two longitudinal cut edges lying in the $z=0$ plane are fully clamped. A concentrated nodal force of 1E+04 is applied in the negative global z-direction at the node A located at the center of the inner surface. Young's modulus and Poisson's ratio are set to 2E+11 and 0.3, respectively. The model is discretized using coarse and refined structured meshes of 70×54×1 and 140×108×2, respectively, where the three numbers denote the numbers of elements in the circumferential, axial, and through-thickness directions. Thus, the refined mesh doubles the number of elements in each of the three directions relative to the coarse mesh.

Table 10 summarizes the maximum z-direction displacement at the loaded point A obtained by SS8 and HCML-SS8 with the coarse and refined meshes, together with the corresponding relative differences. Figure 25 compares the z-direction displacement contours and the bottom-surface von Mises stress contours obtained using the coarse 70×54×1 mesh. The reported displacement and stress values are scaled by E-06 and E+03, respectively.

To quantify the geometric deviation of the actual three-dimensional elements from the training samples, the angle between each thickness edge and the corresponding fitted midsurface normal is adopted as an indicator of thickness-direction distortion. For the ideal solid-shell training samples generated by normal extrusion of the midsurface, this deviation angle is $0^{°}$. A larger deviation angle indicates a more pronounced departure from the canonical normal-extrusion configuration. The thickness-edge deviation angles for the three solid-shell benchmark cases are summarized in Table 11.

Compared with the benchmarks in Section 6.3.1 and Section 6.3.2, the present problem is more challenging because the shell involves opposite curvatures along two orthogonal directions. As a result, the generated solid-shell elements are not uniformly aligned through the thickness and deviate more significantly from the canonical midsurface-extrusion configurations used in the training stage, as shown in Table 11. Nevertheless, the numerical results show that HCML-SS8 remains stable for this more complex geometry. Moreover, the relative difference between the displacements predicted by SS8 and HCML-SS8 decreases further as the mesh is refined from 70×54×1 to 140×108×2. Both the displacement and stress fields preserve the main deformation characteristics of the SS8 reference solution. These results further indicate that the HCML-SS8 retains a certain degree of robustness and limited out-of-distribution generalization capability under geometric conditions that are more complex than the training configurations.

6.3.4. Thin-Walled Pipe Elbow

To assess the applicability of the proposed element in mixed-element analysis, a thin-walled pipe elbow with two straight end sections is considered, as shown in Figure 26. This mixed-element example verifies a defining requirement of numerical-method construction emphasized in Section 1: the learned element can participate in the same global assembly process as conventional finite elements and is not tied to an isolated learned computational domain. The elbow has a centerline radius of $R=50$, inner and outer pipe radii of 20 and 21, respectively, and a corresponding wall thickness of $t=1$. Each straight end section has a length of 50. One end of the pipe is fully clamped, while a concentrated nodal force of $F_z=-1000$ is applied to each node on the annular cross-section at the other end. The two straight pipe sections are discretized using solid-shell elements, with 25×56×2 elements in each section, whereas the elbow is discretized using 33×56×2 Q8R elements. Young's modulus and Poisson's ratio are set to 2E+11 and 0.3.

Figure 27 compares the maximum nodal displacements in the z-direction and the top-surface von Mises stress contours obtained from the two mixed-element models. The reported displacement and stress values are scaled by E-05 and E+04, respectively. In both models, the elbow is discretized using Q8R elements, whereas the two straight pipe sections are modeled using either SS8 or HCML-SS8 elements. The maximum z-direction nodal displacement obtained from the HCML-SS8 and Q8R model differs from that predicted by the SS8 and Q8R model by only 0.24%, indicating close agreement in the predicted global structural flexibility. The two models also produce generally consistent principal stress-concentration regions, spatial distribution patterns, and stress-variation trends, with no evident local stress oscillations or nonphysical stress concentrations near the HCML-SS8 and Q8R interfaces.

7. Discussion and Conclusions

7.1. Discussion

The present results highlight a distinction between stiffness-matrix prediction and finite-element construction. A small regression error does not necessarily lead to a mechanically admissible element, as illustrated by the soft-constraint baseline. For data-driven numerical-method construction, the learned operator must therefore satisfy the structural requirements of the numerical method itself, rather than only approximate reference outputs. In the HCML framework, this is achieved by embedding the physical-informed hard constraints into the stiffness reconstruction through Symbolic derivation.

The results also clarify the different roles of physical(mechanics) and data. The hard constraints determine the admissible structure of the learned stiffness matrix, whereas the reference data determine the remaining element-specific numerical characteristics. This explains why elements learned from Q4, Q4R, Q4I and SS8 retain different numerical behaviors even though they are constructed through the same HCML procedure. The framework should therefore not be interpreted as an automatic improvement of the reference formulation. In this sense, mechanics defines what the learned element must satisfy, while data determine which element is constructed within those requirements.

The use of existing finite-element formulations as reference-data sources should be understood as a controlled setting for examining this construction problem. Their numerical characteristics are known in advance, making it possible to assess whether the HCML framework can reproduce element-level behavior while preserving the prescribed mechanical constraints and the locality, reusability, and assemblability required for finite-element analysis. Importantly, HCML does not rely on the variational formulation, interpolation functions, quadrature rules, stabilization procedures, or other internal details of the reference data; it requires only suitable element-level stiffness information. The origin of that information is therefore conceptually separate from the HCML construction itself.

7.2. Conclusion and Outlook

This work investigated whether data-driven learning can be used to construct a reusable numerical method rather than only approximate the solution of a prescribed boundary-value problem. Using finite elements as a representative setting, we developed a physics-informed Hard-Constrained Manifold Learning (HCML) framework. Under the prescribed linear constraints, symbolic elimination reduces the stiffness description to three independent parameters for the four-node plane-stress element and 78 for the eight-node solid-shell element. A neural network predicts only these independent quantities, while the complete stiffness matrix is recovered through a deterministic reconstruction operator. Numerical tests further show that the resulting HCML-FEs retain the characteristic responses of different reference formulations under bending, locking, mesh distortion, and mesh refinement, while remaining reusable and assemblable in different two- and three-dimensional structural problems.

The present work focuses on the construction problem once suitable element-level reference stiffness information is available. Further development may investigate how such information can be obtained independently of existing finite-element formulations. In particular, deriving discrete element information directly from continuum-mechanics analytical solutions, experimental observations, or other physically consistent sources requires an additional continuum-to-discrete construction procedure and constitutes a separate problem from the HCML formulation considered here. Extensions to nonlinear material behavior, finite deformation, and additional element families also remain to be investigated.