Convolution Hierarchical Deep-learning Neural Networks (C-HiDeNN): finite elements, isogeometric analysis, tensor decomposition, and beyond

Abstract

This paper presents a general Convolution Hierarchical Deep-learning Neural Networks (C-HiDeNN) computational framework for solving partial differential equations. This is the first paper of a series of papers devoted to C-HiDeNN. We focus on the theoretical foundation and formulation of the method. The C-HiDeNN framework provides a flexible way to construct high-order $C^n$ approximation with arbitrary convergence rates and automatic mesh adaptivity. By constraining the C-HiDeNN to build certain functions, it can be degenerated to a specification, the so-called convolution finite element method (C-FEM). The C-FEM will be presented in detail and used to study the numerical performance of the convolution approximation. The C-FEM combines the standard $C^0$ FE shape function and the meshfree-type radial basis interpolation. It has been demonstrated that the C-FEM can achieve arbitrary orders of smoothness and convergence rates by adjusting the different controlling parameters, such as the patch function dilation parameter and polynomial order, without increasing the degrees of freedom of the discretized systems, compared to FEM. We will also present the convolution tensor decomposition method under the reduced-order modeling setup. The proposed methods are expected to provide highly efficient solutions for extra-large scale problems while maintaining superior accuracy. The applications to transient heat transfer problems in additive manufacturing, topology optimization, GPU-based parallelization, and convolution isogeometric analysis have been discussed.

Appendices

Appendix A Brief review of HiDeNN-FEM

The idea of HiDeNN-FEM is to use Deep-learning Neural Networks (DeNN) (We added the “e” instead of just using DNN because we prefer the acronym HiDeNN) to reconstruct the FE shape function by constraining the weights and biases with mesh coordinates, which reads

$$\begin{aligned} \displaystyle {\mathcal {N}}_{{I}}(\varvec{x};\varvec{x}_I^*,\varvec{{\mathcal {A}}}):={\mathcal {F}}_{{I}}(\varvec{x};\varvec{w}(\varvec{x}_I^*),\varvec{b}(\varvec{x}_I^*),\varvec{{\mathcal {A}}}) \end{aligned}$$

(39)

where ${\mathcal {F}}_{{I}}$ stands for the fully connected DeNN structures with weights $\varvec{w}$, biases $\varvec{b}$, and the activation function $\varvec{{\mathcal {A}}}$. ${\mathcal {N}}_{{I}}$ denotes the FE shape function for the node at position $\varvec{x}_I^*$. Assuming a domain $\varOmega$ is discretized by n points, we can write the HiDeNN-FEM approximation as

$$\begin{aligned} \displaystyle u^h(\varvec{x})=\sum _{I=1}^{n}{\mathcal {N}}_{{I}}(\varvec{x};\varvec{x}_I^*,\varvec{{\mathcal {A}}})u_I \end{aligned}$$

(40)

where $u_I$ is the discretized nodal solution of the problem, $u^h$ is the approximated solution function. Considering the vector notation $\varvec{{\mathcal {N}}}=[{\mathcal {N}}_1,\dots ,{\mathcal {N}}_{n}]$ and $\varvec{{u}}=[{u}_1,\dots ,{u}_{n}]^T$, the equation (40) can be simplified as

$$\begin{aligned} \displaystyle u^h(\varvec{x})=\varvec{{\mathcal {N}}}\left( \varvec{x};\varvec{x}^*,\varvec{{\mathcal {A}}}\right) \varvec{{u}} \end{aligned}$$

(41)

The detailed construction of such shape functions using DeNN can be found in [38]. Figure 22 illustrates a detailed architecture of the partially connected DeNN, as an example of HiDeNN-FEM. It should be noticed that the FE shape function is only one of the choices, the HiDeNN structure allows one to easily switch from one to another by releasing the constraints on the weights $\varvec{w}$ and biases $\varvec{b}$.

Fig. 22

HiDeNN-FEM shape functions [38]. Weights and biases are predefined as constants. ReLU (rectified linear unit) activation function is used in the HiDeNN architecture

Since the HiDeNN shape function serves the same role as FE shape function, the derivatives and integration of the shape function can be implemented in exactly the same way as FEM. Finally, the solution of HiDeNN-FEM can be obtained through an optimization problem in which both the displacement field and the mesh coordinates are simultaneously optimized. The problem reads

$$\begin{aligned} \displaystyle u^{\text {HiDeNN-FEM}}=\underset{u_I,\varvec{x}_I^*\in \varOmega \setminus \partial \varOmega }{\arg \min }\ \varPi (u^h(\varvec{x})) \end{aligned}$$

(42)

where $\varOmega$ denotes the entire domain, $\partial \varOmega$ denotes the boundary, $\varPi$ denotes the potential (when it exists) or the residual of the problem. Hence HiDeNN-FEM provides a new way to adaptively modify the mesh and reduces to FEM when the mesh is fixed. Moreover, it is shown that the HiDeNN-FEM has a function approximation space and gives more accurate results than FEM due to the flexibility of the HiDeNN framework [39].

Below, we summarize the key features of the HiDeNN-FEM method [38]:

• Partially connected DeNN with constrained weights and biases

• A flexible framework for function approximation and solving PDEs with automatic r-h-p adaptivity

• More accurate than traditional FEM

• Reduces to FEM when freezing the mesh coordinates

The developed novel Convolution HiDeNN method can be seen as an enhancement for highly smooth solutions and improved convergence rates without increasing the degrees of freedom (DoFs).

Appendix B Derivation of radial basis interpolation function

Let us consider the following 1D case for illustration purposes.

$$\begin{aligned} \begin{aligned} u^{\text {C-FEM}}({\xi })=&\sum _{i\in A^e}{N}_{{i}}({\xi })\sum _{j\in A^i_s} {W}^{\varvec{\xi }_i}_{a,j} ({\xi }) u_j\\ \end{aligned} \end{aligned}$$

(43)

Then we can define the following interpolation as part of the approximation centering around the i-th node in the element domain.

$$\begin{aligned} u^{i}({\xi })=\sum _{j\in A^i_s} {W}^{{\xi }_i}_{a,j} ({\xi }) u_j, \end{aligned}$$

(44)

where the supporting node set of W is $A^i_s$ with a given patch size s. Figure 23 illustrates a 1D convolution element with patch size $s=1$ and the two-node shape functions for $N_i$. Now the question is how to compute the W based on the given supporting nodes.

Fig. 23

Supporting nodes for a 1D convolution element with patch size $s=1$

Assuming the nodal solution value for the 4 nodes in Fig. 23 is $[u_1,u_2,u_3,u_4]$, we illustrate the radial basis interpolation procedure for the part centering around $i=2$. In this case, the parametric coordinates for the support nodes are $\{-3,-1,1\}$. Then we can consider the radial basis interpolation $u^{i=2}({\xi })$ has the following form

$$\begin{aligned} u^{i}({\xi })=\varvec{\Psi }_a (\xi ){\textbf{k}}+{\textbf{p}}(\xi ){\textbf{l}}, \end{aligned}$$

(45)

where $\varvec{\Psi }_a$ is a defined kernel function, which can be the reproducing kernel or cubic spline kernel [3, 4] with the dilation parameter a, ${\textbf{p}}(\xi )$ is the polynomial basis vector of p-th order, ${\textbf{k}}=[k_1, k_2, k_3]^T$ and ${\textbf{l}}=[l_1, l_2, l_3]^T$ are the coefficient vector that helps to enforce the reproducing condition and Kronecker delta property. We give here a specific example for $\varvec{\Psi }_a$ and ${\textbf{p}}(\xi )$ using a cubic spline kernel and a second-order polynomial.

$$\begin{aligned}&\varvec{\Psi }_a (\xi )=[{\Psi }_a (\xi -\xi _1),\ {\Psi }_a (\xi -\xi _2) , \ {\Psi }_a (\xi -\xi _3)]\nonumber \\&\text {where}\quad {\Psi }_a (\xi -\xi _I): = {\Psi }_a (z) \quad \text {with} \quad z=\frac{|\xi -\xi _I|}{a}\nonumber \\&={\left\{ \begin{array}{ll} \frac{2}{3} -4z^2+4z^3 \quad \forall z \in [0, \frac{1}{2}]\\ \frac{4}{3} -4z+4z^2- \frac{4}{3}z^3 \quad \forall z \in [\frac{1}{2}, 1]\\ 0 \quad \forall z \in (1,+\infty ) \end{array}\right. }, \end{aligned}$$

(46)

and

$$\begin{aligned} \begin{aligned} {\textbf{p}}=[1,\ \xi ,\ \xi ^2] \end{aligned} \end{aligned}$$

(47)

Now we can compute ${\textbf{k}}$ and ${\textbf{l}}$ by enforcing the below conditions.

$$\begin{aligned} \begin{aligned} {\left\{ \begin{array}{ll} u^{i}({\xi _1})= u_1 \\ u^{i}({\xi _2})= u_2 \\ u^{i}({\xi _3})= u_3 \\ \sum {\textbf{k}}= 0\\ {[}\xi _1, \xi _2, \xi _3]\ {\textbf{k}}=0\\ {[}\xi _1^2, \xi _2^2, \xi _3^2]\ {\textbf{k}}=0\\ \end{array}\right. }, \end{aligned} \end{aligned}$$

(48)

Solving the above equations gives the solution to ${\textbf{k}}$ and ${\textbf{l}}$, which reads

$$\begin{aligned} \begin{aligned} {\left\{ \begin{array}{ll} {\textbf{k}}{=}{\textbf{K}}{\textbf{u}}\\ {\textbf{l}}{=}{\textbf{L}}{\textbf{u}}\\ \end{array}\right. }, \end{aligned} \end{aligned}$$

(49)

with

$$\begin{aligned} \begin{aligned} {\left\{ \begin{array}{ll} {\textbf{u}}{=}[u_1,\ u_2,\ u_3]^T\\ {\textbf{L}}{=}({\textbf{P}}^T{\textbf{R}}_0{\textbf{P}})^{-1}{\textbf{P}}^T{\textbf{R}}_0^{-1}\\ {\textbf{K}}{=}{\textbf{R}}_0^{-1}({\textbf{I}}-{\textbf{P}}{\textbf{L}}) \end{array}\right. }, \end{aligned} \end{aligned}$$

(50)

and

$$\begin{aligned} {\left\{ \begin{array}{ll} {\textbf{R}}_0{=}\begin{pmatrix} \varvec{\Psi }_a (\xi _1)\\ \varvec{\Psi }_a (\xi _2)\\ \varvec{\Psi }_a (\xi _3) \end{pmatrix}{=}\begin{pmatrix} {\Psi }_a (\xi _1{-}\xi _1) &{}\ {\Psi }_a (\xi _1{-}\xi _2) &{}\ {\Psi }_a (\xi _1{-}\xi _3)\\ {\Psi }_a (\xi _2{-}\xi _1) &{}\ {\Psi }_a (\xi _2{-}\xi _2) &{}\ {\Psi }_a (\xi _2{-}\xi _3)\\ {\Psi }_a (\xi _3{-}\xi _1) &{}\ {\Psi }_a (\xi _3{-}\xi _2) &{}\ {\Psi }_a (\xi _3-\xi _3) \end{pmatrix}\\ \\ {\textbf{P}}{=}\begin{pmatrix} {\textbf{p}} (\xi _1)\\ {\textbf{p}} (\xi _2)\\ {\textbf{p}} (\xi _3) \end{pmatrix}=\begin{pmatrix} 1 &{}\ \xi _1 &{}\ \xi _1^2\\ 1 &{}\ \xi _2 &{}\ \xi _2^2\\ 1 &{}\ \xi _3 &{}\ \xi _3^2\\ \end{pmatrix}\\ \end{array}\right. }\!\!\!, \end{aligned}$$

(51)

Finally, the radial basis interpolation with the computed coefficients reads

$$\begin{aligned} u^{i}({\xi }) {=} \varvec{\Psi }_a (\xi ){\textbf{k}}{+}{\textbf{p}}(\xi ){\textbf{l}}= & {} \varvec{\Psi }_a (\xi ){\textbf{K}}{\textbf{u}}+{\textbf{p}}(\xi ){\textbf{L}}{\textbf{u}}\nonumber \\= & {} (\varvec{\Psi }_a (\xi ){\textbf{K}}+{\textbf{p}}(\xi ){\textbf{L}}){\textbf{u}}\nonumber \\= & {} {W}^{{\xi }_i}_{a,1}(\xi ) u_1+{W}^{{\xi }_i}_{a,2} (\xi )u_2 \nonumber \\{} & {} \quad +{W}^{{\xi }_i}_{a,3}(\xi ) u_3\nonumber \\= & {} \sum _{j\in A^i_s} {W}^{{\xi }_i}_{a,j} ({\xi }) u_j , \end{aligned}$$

(52)

where ${W}^{{\xi }_i}_{a,j}$ is obtained by identifying the corresponding coefficient of $u_j$. By analogy, we can compute the other convolution patch functions W with the support $A^{i=3}_s$. Detailed mathematical derivation and analysis of the radial basis interpolation can be found in [77].

Appendix C Illustration of a 1D convolution element

For a better understanding of the convolution shape function $\tilde{{N}}_k$, we illustrate here a 1D convolution approximation with ${N}_{{i}}$ chosen to be linear and patch size $s=1$, see again Fig. 23 for the supporting nodes of one convolution element.

$$\begin{aligned} u^{\text {C-FEM}}({\xi })=\sum _{i\in A^e}{N}_{{i}}({\xi })\sum _{j\in A^i_s} {W}^{{\xi }_i}_{a,j} ({\xi }) u_j \end{aligned}$$

(53)

In this case, the FE shape function nodal support set $A^e=\{2,3\}$, and then $A^{i=2}_s=\{1,2,3\}$, $A^{i=3}_s=\{2,3,4\}$. The equation (53) becomes

$$\begin{aligned} u^{\text {C-FEM}}({\xi })= & {} \sum _{i\in A^e}{N}_{{i}}({\xi })\sum _{j\in A^i_s} {W}^{{\xi }_i}_{a,j} ({\xi }) u_j \nonumber \\= & {} {N}_{{2}}(\xi ){W}^{{\xi }_2}_{a,1}(\xi )u_1+({N}_{{2}}(\xi ){W}^{{\xi }_2}_{a,2}(\xi )\nonumber \\{} & {} \quad +{N}_{{3}}(\xi ){W}^{{\xi }_3}_{a,2}(\xi ))u_2\nonumber \\{} & {} \quad +({N}_{{2}}(\xi ){W}^{{\xi }_2}_{a,3}(\xi )\nonumber \\{} & {} \quad +{N}_{{3}}(\xi ){W}^{{\xi }_3}_{a,3}(\xi ))u_3+{N}_{{3}}(\xi ){W}^{{\xi }_3}_{a,4}(\xi )u_4 \nonumber \\= & {} \sum _{k\in A_s} \tilde{{N}}_k({\xi }) u_k, \end{aligned}$$

(54)

where $A_s=\bigcup _{i\in A^e} A^{i}_s =\{1,2,3,4\}$. Therefore, there are in total 4 convolution shape functions for $s=1$. If $s=2$, we can expect 6 shape functions, as shown in Fig. 5.

Appendix D Modification of the natural coordinates for irregular meshes

The modification of the natural (parametric) coordinates of the patch nodes can be done according to the distance ratio in the physical domain. As shown in Fig. 24, if the mesh is irregular in the physical space, the corresponding patch nodes in parametric space are adapted. The neighbouring nodes $x_1$ and $x_4$ are clearly far from the element in the physical space. In order to reflect this distancing information in parametric domain. The coordinates $\xi _1$ and $\xi _4$ are modified as below

$$\begin{aligned} \xi _{1}=\xi _{2}-|\xi _2 - \xi _3|\frac{|x_1-x_2|}{|x_2-x_3|}, \quad \xi _{4}=\xi _{3}+|\xi _2 - \xi _3|\frac{|x_3-x_4|}{|x_2-x_3|}, \end{aligned}$$

(55)

that is

$$\begin{aligned} \frac{|\xi _{1}-\xi _{2}|}{|\xi _2 - \xi _3|}=\frac{|x_1-x_2|}{|x_2-x_3|}, \quad \frac{|\xi _{3}-\xi _{4}|}{|\xi _2 - \xi _3|}=\frac{|x_3-x_4|}{|x_2-x_3|}, \end{aligned}$$

(56)

where the right-hand-side of the above equation is the physical distance ratio and left-hand-side is the one in parametric space. By enforcing this condition, we are able to keep the nonuniform distribution of the nodes in parametric space. Here, the original physical element domain is $[x_2, x_3]$ and the corresponding parametric element domain $[\xi _2, \xi _3]$ remains the same as before, i.e., $[-1,1]$.

Fig. 24

Modified parametric coordinates for a 1D convolution element with patch size $s=1$

From our experience this modification is necessary and can help achieve the optimal convergence rates when irregular meshes are used. Similar concept can be used in 2D/3D irregular meshes.

Appendix E Error estimate of C-HiDeNN interpolation

Following the procedure in the previous work on the meshfree reproducing property [8, 27], we can derive the error estimate of C-HiDeNN interpolation, based on the reproducing property.

Theorem 1

Let us assume an exact solution $u(x)\in C^{p+1}$$(\bar{\varOmega })\cap H^{p+1} (\varOmega )$, where $\varOmega$ is a bounded open set in ${\mathbb {R}}^n$. Define an interpolation operator ${\mathcal {I}}$

$$\begin{aligned} {\mathcal {I}} u(x)=\sum _I \tilde{N}_I(x) u(x_I), \end{aligned}$$

(57)

where the interpolation operator satisfies the reproducing property

$$\begin{aligned} {\mathcal {I}} x^m = \sum _I \tilde{N}_I(x) (x_I)^m = x^m, m=0,1,2,\dots ,p. \end{aligned}$$

(58)

In the refinement, the mesh is given and refined at the same scale. The shape function $\tilde{N}_I(x)$ can be regarded as a scaled reference shape function $\tilde{N}^0(\xi )$ in the parent domain with $\xi =(x-x_I)/h$. Suppose the boundary $\partial \varOmega$ is smooth enough, then the following interpolation estimate holds

$$\begin{aligned} \Vert u(x)-{\mathcal {I}} u \Vert _{H^1(\varOmega )} \le C_1 h^p \Vert u\Vert _{H^{p+1}(\varOmega )}. \end{aligned}$$

(59)

where h is the mesh size and $C_1$ is a constant independent of h.

Proof

We have

$$\begin{aligned} u(x)-{\mathcal {I}} u= & {} u(x)-\sum _I \tilde{N}_I(x) u(x_I) = u(x)\nonumber \\{} & {} \quad -\sum _{I\in \varLambda (x)} \tilde{N}_I(x) u(x_I), \end{aligned}$$

(60)

where $\varLambda (x)=\{ I| x\in \text {supp}\{\tilde{N}_I\}\cap \varOmega \}$ is an index set supporting all the shape functions covering x.

Take Taylor expansion of $u(x_I)$ at x:

$$\begin{aligned} u(x_I)= & {} u(x)+\sum _{m=1}^{p} \dfrac{1}{m!}(x_I-x)^m D^m u(x)\nonumber \\{} & {} \quad {+} \dfrac{1}{(p{+}1)!} (x_I{-}x)^{p{+}1} D^{p{+}1} u(x{+}\theta (x_I{-}x)),\nonumber \\ \end{aligned}$$

(61)

where $D^m u$ represents the m-th order derivative of u(x), and $0<\theta <1$ depends on x and $x_I$. Substituting (61) into (60) yields

$$\begin{aligned}{} & {} \!\!\!u(x)-{\mathcal {I}} u \end{aligned}$$

(62)

$$\begin{aligned}{} & {} \quad {=} u(x){-}\sum _{I\in \varLambda (x)} \tilde{N}_I(x) u(x_I) \nonumber \\{} & {} \quad {=} u(x){-}\sum _{I\in \varLambda (x)} \tilde{N}_I(x) \bigg (u(x){+}\sum _{m=1}^{p} \dfrac{1}{m!}(x_I{-}x)^m D^m u(x) \nonumber \\{} & {} \qquad {+} \dfrac{1}{(p{+}1)!} (x_I{-}x)^{p{+}1} D^{p{+}1} u(x+\theta (x_I-x))\bigg ) \nonumber \\{} & {} \quad = u(x)-\bigg ( u(x)\sum _{I\in \varLambda (x)} \tilde{N}_I(x) \nonumber \\{} & {} \qquad + \dfrac{1}{m!} D^m u(x) \sum _{m=1}^{p} \sum _{I\in \varLambda (x)} \tilde{N}_I(x) (x_I-x)^m \nonumber \\{} & {} \qquad {+}\!\!\!\sum _{I\in \varLambda (x)}\!\!\!\dfrac{1}{(p{+}1)!} (x_I{-}x)^{p{+}1} D^{p{+}1} u(x{+}\theta (x_I{-}x)\!) \tilde{N}_I(x) \bigg )\nonumber \\ \end{aligned}$$

(63)

According to the reproducing property, we have

$$\begin{aligned} {\mathcal {I}} (y{-}x)^m= & {} \sum _I \tilde{N}_I(x) (y{-}x_I)^m {=} \sum _{I\in \varLambda (x)} \tilde{N}_I(x) (y{-}x_I)^m \nonumber \\= & {} (y-x)^m, m=1,2,\dots ,p, \end{aligned}$$

(64)

particularly, for $y=x$

$$\begin{aligned} \sum _{I\in \varLambda (x)} \tilde{N}_I(x) (x-x_I)^m = 0, m=1,2,\dots ,p. \end{aligned}$$

(65)

So the deviation between u(x) and its approximation ${\mathcal {I}}u$ is

$$\begin{aligned}{} & {} \!\!\!u(x)-{\mathcal {I}} u \nonumber \\{} & {} \quad {=}{-}\!\!\sum _{I\in \varLambda (x)} \dfrac{1}{(p{+}1)!} (x_I{-}x)^{p{+}1} D^{p{+}1} u(x{+}\theta (x_I{-}x)) \tilde{N}_I(x)\nonumber \\ \end{aligned}$$

(66)

It deduces the following estimate

$$\begin{aligned}{} & {} | u(x)-{\mathcal {I}} u | \le \dfrac{1}{(p+1)!} \sum _{I\in \varLambda (x)} |x_I-x|^{p+1} |D^{p+1}\nonumber \\{} & {} \quad u(x+\theta (x_I-x))| |\tilde{N}_I(x)| \end{aligned}$$

(67)

Denote the support size of shape function $\tilde{N}_I(x)$ as rh, the maximum value of $|\tilde{N}_I(x)|$ is M. Note that the absolute value of C-HiDeNN shape function is generally less than 1 (in the element domain). The following estimate holds

$$\begin{aligned} | u(x){-}{\mathcal {I}} u | {\le } \dfrac{(rh)^{p{+}1}}{(p{+}1)!}M \sum _{I\in \varLambda (x)} |D^{p+1} u(x{+}\theta (x_I{-}x))|\nonumber \\ \end{aligned}$$

(68)

We can conclude that $\exists \ C_0, \text {with}\ 0< C_0 < \infty$, such that

$$\begin{aligned} \Vert u(x){-}{\mathcal {I}} u \Vert _{L^2} {\le } C_0 h^{p+1} \Vert u\Vert _{H^{p{+}1}(\varOmega )}. \end{aligned}$$

(69)

Next, we consider the derivative of $u(x)-{\mathcal {I}} u$. Taking the derivative of ${\mathcal {I}} u$ yields

$$\begin{aligned} D^1 {\mathcal {I}} u = \sum _I D^1 \tilde{N}_I(x) u(x_I). \end{aligned}$$

(70)

The derivative of $u(x)-{\mathcal {I}} u$ is

$$\begin{aligned}{} & {} D^1 u(x)- D^1 {\mathcal {I}} u \end{aligned}$$

(71)

$$\begin{aligned}{} & {} \quad = D^1 u(x)-\sum _{I\in \varLambda (x)} D^1 \tilde{N}_I(x) u(x_I) \nonumber \\{} & {} \quad = D^1 u(x)-\sum _{I\in \varLambda (x)} D^1 \tilde{N}_I(x) \bigg (u(x)+\sum _{m=1}^{p} \dfrac{1}{m!}(x_I-x)^m\nonumber \\{} & {} \qquad \times D^m u(x) {+} \dfrac{1}{(p{+}1)!} (x_I{-}x)^{p{+}1} D^{p{+}1} u(x{+}\theta (x_I{-}x))\bigg ) \nonumber \\{} & {} \quad = u(x)-\bigg ( u(x)\sum _{I\in \varLambda (x)} D^1 \tilde{N}_I(x) + \dfrac{1}{m!} D^m u(x) \nonumber \\{} & {} \qquad \times \sum _{m=1}^{p} \sum _{I\in \varLambda (x)} D^1 \tilde{N}_I(x) (x_I-x)^m \nonumber \\{} & {} \qquad + \sum _{I\in \varLambda (x)} \dfrac{1}{(p+1)!} (x_I-x)^{p+1} D^{p+1}\nonumber \\{} & {} \qquad \times u(x+\theta (x_I-x)) D^1 \tilde{N}_I(x) \bigg ) \end{aligned}$$

(72)

Taking the derivative of (64) with respect to x yields

$$\begin{aligned} D^1 {\mathcal {I}} (y-x)^m= & {} \sum _{I\in \varLambda (x)} D^1 \tilde{N}_I(x) (y-x_I)^m \nonumber \\= & {} - m (y-x)^{m-1}, m=1,2,...,p \end{aligned}$$

(73)

For $y=x$, we have

$$\begin{aligned} \sum _{I\in \varLambda (x)} D^1 \tilde{N}_I(x) (x-x_I)^m = - \delta _{1,m}, m=1,2,...,p \end{aligned}$$

(74)

This leads to

$$\begin{aligned}{} & {} D^1 u(x)- D^1 {\mathcal {I}} u\nonumber \\{} & {} \quad = - \sum _{I\in \varLambda (x)} \dfrac{1}{(p+1)!} (x_I-x)^{p+1} D^{p+1}\nonumber \\{} & {} \quad \times u(x+\theta (x_I-x)) D^1 \tilde{N}_I(x) \end{aligned}$$

(75)

The reference shape function $\tilde{N}^0(\xi )$ is defined in the parent domain with $\xi =(x-x_I)/h$, independent of the element size h. Thus

$$\begin{aligned} D^1 \tilde{N}_I(x) = \dfrac{1}{h} D^1_{\xi } \tilde{N}^0(\xi ). \end{aligned}$$

(76)

The estimate for the derivative of $u(x)-{\mathcal {I}} u$ becomes

$$\begin{aligned}{} & {} |D^1 u(x)- D^1 {\mathcal {I}} u|\nonumber \\{} & {} \quad = \left| \sum _{I\in \varLambda (x)} \dfrac{1}{(p+1)!h} (x_I-x)^{p+1} D^{p+1}\right. \nonumber \\{} & {} \quad \left. u(x+\theta (x_I-x)) D^1_{\xi } \tilde{N}^0(\xi ) \right| \nonumber \\{} & {} \quad \le \dfrac{r^{p+1}h^p}{(p+1)!}K \sum _{I\in \varLambda (x)} |D^{p+1} u(x+\theta (x_I-x))|, \end{aligned}$$

(77)

where K is the maximum value of $|D^1_\xi \tilde{N}^0(\xi )|$. We can conclude that $\exists \ C_1 \text {with}\ 0< C_1 < \infty$, such that

$$\begin{aligned} \Vert u(x)-{\mathcal {I}} u \Vert _{H^1} \le C_1 h^{p} \Vert u\Vert _{H^{p+1}(\varOmega )}. \end{aligned}$$

(78)

. $\square$

For a linear second-order elliptic boundary value problem, we can use Céa’s inequality [27], therefore

$$\begin{aligned} \Vert u^{\text {C{-}HiDeNN}}(\varvec{x}){-}u^{\text {Ext}}(\varvec{x}) \Vert _{H^1}\le & {} c \inf _{v^h\in V^h}\Vert v^h(\varvec{x}){-}u^{\text {Ext}}(\varvec{x}) \Vert _{H^1}\nonumber \\\le & {} c\Vert {\mathcal {I}}u^{\text {Ext}}{-}u^{\text {Ext}}(\varvec{x}) \Vert _{H^1}\nonumber \\\le & {} c \times C_1 h^p \Vert u^{\text {Ext}} \Vert _{H^{p+1}}, \end{aligned}$$

(79)

where $V^h$ represents the C-HiDeNN interpolation space.

Appendix F Illustration of C-PGD/TD solution procedure

The C-PGD solution can be solved through the so-called alternating fixed point algorithm as traditional PGD method [40], or through a minimization problem [44]. Let us consider a 3D Poisson problem as below

$$\begin{aligned} \nabla ^2 u(x,y,z)+b(x,y,z)=0 \ \text {in}\ \varOmega \quad \text {with}\ \nabla u.\varvec{n}=0 \ \text {on}\ \partial \varOmega \end{aligned}$$

(80)

where $\nabla$ denotes the gradient operator, b is the body source term with the assumption that $b(x,y,z)=b_x(x)b_y(y)b_z(z)$. The domain $\varOmega$ is a regular domain, $\varvec{n}$ is the normal vector on the surface. The weak form of the problem can be obtained by multiplying both sides of the above equation by the test function $\delta u$, which reads

$$\begin{aligned} \int _\varOmega (\nabla \delta u)^T \nabla u\ d\varOmega -\int _\varOmega \delta u\ b\ d\varOmega =0 \end{aligned}$$

(81)

Assuming the solution can be accurately approximated by M modes and $M-1$ modes are already computed, the following decomposition for u and $\delta u$ read

$$\begin{aligned} u(x,y,z)= & {} \sum _{m=1}^{M-1} u^{(m)}_x(x) u^{(m)}_y(y) u^{(m)}_z(z)\nonumber \\{} & {} \quad +u^{(M)}_x(x) u^{(M)}_y(y) u^{(M)}_z(z) \end{aligned}$$

(82)

and

$$\begin{aligned} \delta u(x,y,z)= & {} \delta u^{(M)}_x u^{(M)}_y u^{(M)}_z + u^{(M)}_x \delta u^{(M)}_y u^{(M)}_z\nonumber \\{} & {} \quad + u^{(M)}_x u^{(M)}_y \delta u^{(M)}_z \end{aligned}$$

(83)

For notation simplification, we can omit the superscript M in the above equations. Therefore,

$$\begin{aligned} u(x,y,z)=\sum _{m=1}^{M-1} u^{(m)}_x(x) u^{(m)}_y(y) u^{(m)}_z(z) +u_x(x) u_y(y) u_z(z) \end{aligned}$$

(84)

and

$$\begin{aligned} \delta u(x,y,z)= \delta u_x u_y u_z + u_x \delta u_y u_z+ u_x u_y \delta u_z \end{aligned}$$

(85)

The convolution approximation is then applied to each of the separated 1D functions $u_x(x), u_y(y), u_z(z)$, which reads

$$\begin{aligned} \begin{aligned} u_x(x)=\tilde{\varvec{N}}(x)\varvec{u}_x,\ u_y(y)=\tilde{\varvec{N}}(y)\varvec{u}_y,\ u_z(z)=\tilde{\varvec{N}}(z)\varvec{u}_z \end{aligned} \end{aligned}$$

(86)

where $\tilde{\varvec{N}}$ is the convolution shape function vector formed by the 1D convolution functions in a patch domain. $\varvec{u}_x, \varvec{u}_y, \varvec{u}_z$ are the associated nodal solution vectors in three directions. Similarly, the same convolution approximation can be applied to $\delta {u}_x, \delta {u}_y, \delta {u}_z$ and ${u}_x^{(m)}, {u}_y^{(m)}, {u}_z^{(m)}$.

The $u_x(x), u_y(y), u_z(z)$ can be computed by alternatively fixing two of the functions. For example, we can assume $u_y, u_z$ are given by assumed values, then $\delta u_y=0$ and $\delta u_z=0$, $\delta u(x,y,z)= \delta u_x u_y u_z$, the weak form (81) for solving $u_x$ under the decomposition reads

$$\begin{aligned} \begin{aligned}&\int _\varOmega (\nabla (\delta u_x u_y u_z))^T \nabla \left( \sum _{m=1}^{M-1} u^{(m)}_x u^{(m)}_y u^{(m)}_z {+}u_x u_y u_z\right) \ d\varOmega \\&\quad -\int _\varOmega \delta u_x u_y u_z \ b\ d\varOmega =0 \end{aligned} \end{aligned}$$

(87)

with

$$\begin{aligned} \nabla (\delta u_x u_y u_z)= & {} \left[ \frac{\partial \delta u_x}{\partial x} u_y u_z,\ \delta u_x \frac{\partial u_y}{\partial y}u_z, \delta u_x u_y \frac{\partial u_z}{\partial z}\right] ^T \nonumber \\= & {} \left[ \tilde{\varvec{B}}_x\delta \varvec{u}_x \tilde{\varvec{N}}_y\varvec{u}_y \tilde{\varvec{N}}_z\varvec{u}_z,\ \tilde{\varvec{N}}_x\delta \varvec{u}_x \tilde{\varvec{B}}_y\varvec{u}_y \tilde{\varvec{N}}_z\varvec{u}_z,\right. \nonumber \\{} & {} \left. \tilde{\varvec{N}}_x\delta \varvec{u}_x \tilde{\varvec{N}}_y\varvec{u}_y \tilde{\varvec{B}}_z\varvec{u}_z\right] ^T\nonumber \\ \end{aligned}$$

(88)

and

$$\begin{aligned} \nabla ( u^{(m)}_x u^{(m)}_y u^{(m)}_z)= & {} \left[ \frac{\partial u^{(m)}_x}{\partial x} u^{(m)}_y u^{(m)}_z,\ u^{(m)}_x \frac{\partial u^{(m)}_y}{\partial y}u^{(m)}_z,\right. \nonumber \\{} & {} \left. \quad u^{(m)}_x u^{(m)}_y \frac{\partial u^{(m)}_z}{\partial z}\right] ^T\nonumber \\= & {} \left[ \tilde{\varvec{B}}_x\varvec{u}_x^{(m)} \tilde{\varvec{N}}_y\varvec{u}_y^{(m)} \tilde{\varvec{N}}_z\varvec{u}_z^{(m)}\right. \nonumber \\ {}{} & {} \left. \quad \tilde{\varvec{N}}_x\varvec{u}_x^{(m)} \tilde{\varvec{B}}_y\varvec{u}_y^{(m)} \tilde{\varvec{N}}_z\varvec{u}_z^{(m)}, \right. \nonumber \\ {}{} & {} \left. \quad \ \tilde{\varvec{N}}_x\varvec{u}_x^{(m)} \tilde{\varvec{N}}_y\varvec{u}_y^{(m)} \tilde{\varvec{B}}_z\varvec{u}_z^{(m)}\right] ^T \end{aligned}$$

(89)

where $\tilde{\varvec{B}}_x=\partial \tilde{\varvec{N}}_x/\partial x$, $\tilde{\varvec{B}}_y=\partial \tilde{\varvec{N}}_y/\partial y$, $\tilde{\varvec{B}}_z=\partial \tilde{\varvec{N}}_z/\partial z$ Therefore, the Eq. (87) becomes

$$\begin{aligned}{} & {} \sum _{m=1}^{M-1}\int _\varOmega \delta \varvec{u}_x^T\tilde{\varvec{B}}_x^T\tilde{\varvec{B}}_x\varvec{u}_x^{(m)}\varvec{u}_y^T\tilde{\varvec{N}}_y^T\tilde{\varvec{N}}_y\varvec{u}_y^{(m)}\varvec{u}_z^T\tilde{\varvec{N}}_z^T\tilde{\varvec{N}}_z\varvec{u}_z^{(m)}\ d\varOmega \nonumber \\{} & {} \quad {+}\!\sum _{m{=}1}^{M{-}1}\int _\varOmega \delta \varvec{u}_x^T\tilde{\varvec{N}}_x^T\tilde{\varvec{N}}_x\varvec{u}_x^{(m)}\varvec{u}_y^T\tilde{\varvec{B}}_y^T\tilde{\varvec{B}}_y\varvec{u}_y^{(m)}\varvec{u}_z^T\tilde{\varvec{N}}_z^T\tilde{\varvec{N}}_z\varvec{u}_z^{(m)}\ d\varOmega \nonumber \\{} & {} \quad +\sum _{m=1}^{M-1}\int _\varOmega \delta \varvec{u}_x^T\tilde{\varvec{N}}_x^T\tilde{\varvec{N}}_x\varvec{u}_x^{(m)}\varvec{u}_y^T\tilde{\varvec{N}}_y^T\tilde{\varvec{N}}_y\varvec{u}_y^{(m)}\varvec{u}_z^T\tilde{\varvec{B}}_z^T\tilde{\varvec{B}}_z\varvec{u}_z^{(m)}\ d\varOmega \nonumber \\{} & {} \quad +\int _\varOmega \delta \varvec{u}_x^T\tilde{\varvec{B}}_x^T\tilde{\varvec{B}}_x\varvec{u}_x\varvec{u}_y^T\tilde{\varvec{N}}_y^T\tilde{\varvec{N}}_y\varvec{u}_y\varvec{u}_z^T\tilde{\varvec{N}}_z^T\tilde{\varvec{N}}_z\varvec{u}_z\ d\varOmega \nonumber \\{} & {} \quad +\int _\varOmega \delta \varvec{u}_x^T\tilde{\varvec{N}}_x^T\tilde{\varvec{N}}_x\varvec{u}_x\varvec{u}_y^T\tilde{\varvec{B}}_y^T\tilde{\varvec{B}}_y\varvec{u}_y\varvec{u}_z^T\tilde{\varvec{N}}_z^T\tilde{\varvec{N}}_z\varvec{u}_z\ d\varOmega \nonumber \\{} & {} \quad +\int _\varOmega \delta \varvec{u}_x^T\tilde{\varvec{N}}_x^T\tilde{\varvec{N}}_x\varvec{u}_x\varvec{u}_y^T\tilde{\varvec{N}}_y^T\tilde{\varvec{N}}_y\varvec{u}_y\varvec{u}_z^T\tilde{\varvec{B}}_z^T\tilde{\varvec{B}}_z\varvec{u}_z\ d\varOmega \nonumber \\{} & {} \quad -\int _\varOmega \delta \varvec{u}_x^T\tilde{\varvec{N}}_x^T b_x\varvec{u}_y^T\tilde{\varvec{N}}_y^T b_y\varvec{u}_z^T\tilde{\varvec{N}}_z^T b_z\ d\varOmega \nonumber \\{} & {} \quad =0 \end{aligned}$$

(90)

The final discretized form for solving $u_x$ is

$$\begin{aligned}{} & {} \sum _{m=1}^{M-1}\tilde{\varvec{K}}_{xx}\varvec{u}_x^{(m)}\varvec{u}_y^T\tilde{\varvec{M}}_{yy}\varvec{u}_y^{(m)}\varvec{u}_z^T\tilde{\varvec{M}}_{zz}\varvec{u}_z^{(m)}\nonumber \\{} & {} \quad +\sum _{m=1}^{M-1}\tilde{\varvec{M}}_{xx}\varvec{u}_x^{(m)}\varvec{u}_y^T\tilde{\varvec{K}}_{yy}\varvec{u}_y^{(m)}\varvec{u}_z^T\tilde{\varvec{M}}_{zz}\varvec{u}_z^{(m)}\nonumber \\{} & {} \quad +\sum _{m=1}^{M-1}\tilde{\varvec{M}}_{xx}\varvec{u}_x^{(m)}\varvec{u}_y^T\tilde{\varvec{M}}_{yy}\varvec{u}_y^{(m)}\varvec{u}_z^T\tilde{\varvec{K}}_{zz}\varvec{u}_z^{(m)}\nonumber \\{} & {} \quad +\tilde{\varvec{K}}_{xx}\varvec{u}_x\varvec{u}_y^T\tilde{\varvec{M}}_{yy}\varvec{u}_y\varvec{u}_z^T\tilde{\varvec{M}}_{zz}\varvec{u}_z\nonumber \\{} & {} \quad +\tilde{\varvec{M}}_{xx}\varvec{u}_x\varvec{u}_y^T\tilde{\varvec{K}}_{yy}\varvec{u}_y\varvec{u}_z^T\tilde{\varvec{M}}_{zz}\varvec{u}_z\nonumber \\{} & {} \quad +\tilde{\varvec{M}}_{xx}\varvec{u}_x\varvec{u}_y^T\tilde{\varvec{M}}_{yy}\varvec{u}_y\varvec{u}_z^T\tilde{\varvec{K}}_{zz}\varvec{u}_z\nonumber \\{} & {} \quad -\tilde{\varvec{Q}}_x\varvec{u}_y^T\tilde{\varvec{Q}}_y \varvec{u}_z^T\tilde{\varvec{Q}}_z\nonumber \\{} & {} \quad =0 \end{aligned}$$

(91)

where $\tilde{\varvec{K}}_{xx}=\int _{\varOmega _x} \tilde{\varvec{B}}_x^T\tilde{\varvec{B}}_x \ dx$, $\tilde{\varvec{K}}_{yy}=\int _{\varOmega _y} \tilde{\varvec{B}}_y^T\tilde{\varvec{B}}_y \ dy$, $\tilde{\varvec{K}}_{zz}=\int _{\varOmega _z} \tilde{\varvec{B}}_z^T\tilde{\varvec{B}}_z \ dz$, $\tilde{\varvec{M}}_{xx}=\int _{\varOmega _x} \tilde{\varvec{N}}_x^T\tilde{\varvec{N}}_x \ dx$, $\tilde{\varvec{M}}_{yy}=\int _{\varOmega _y} \tilde{\varvec{N}}_y^T\tilde{\varvec{N}}_y \ dy$, $\tilde{\varvec{M}}_{zz}=\int _{\varOmega _z} \tilde{\varvec{N}}_z^T\tilde{\varvec{N}}_z \ dz$, with $\varOmega =\varOmega _x\times \varOmega _y\times \varOmega _z$.

Rearranging the above equation leads to the following linear system of equations

$$\begin{aligned} \tilde{\varvec{K}}\varvec{u}_x=\tilde{\varvec{Q}} \end{aligned}$$

(92)

with

$$\begin{aligned} \tilde{\varvec{K}}= & {} \tilde{\varvec{K}}_{xx}\varvec{u}_y^T\tilde{\varvec{M}}_{yy}\varvec{u}_y\varvec{u}_z^T\tilde{\varvec{M}}_{zz}\varvec{u}_z+\tilde{\varvec{M}}_{xx}\varvec{u}_y^T\tilde{\varvec{K}}_{yy}\varvec{u}_y\varvec{u}_z^T\tilde{\varvec{M}}_{zz}\varvec{u}_z\nonumber \\{} & {} \quad +\tilde{\varvec{M}}_{xx}\varvec{u}_y^T\tilde{\varvec{M}}_{yy}\varvec{u}_y\varvec{u}_z^T\tilde{\varvec{K}}_{zz}\varvec{u}_z\nonumber \\ \tilde{\varvec{Q}}= & {} \tilde{\varvec{Q}}_x\varvec{u}_y^T\tilde{\varvec{Q}}_y \varvec{u}_z^T\tilde{\varvec{Q}}_z{-}\!\!\sum _{m=1}^{M-1}\tilde{\varvec{K}}_{xx}\varvec{u}_x^{(m)}\varvec{u}_y^T\tilde{\varvec{M}}_{yy}\varvec{u}_y^{(m)}\varvec{u}_z^T\tilde{\varvec{M}}_{zz}\varvec{u}_z^{(m)}\nonumber \\{} & {} \quad -\sum _{m=1}^{M-1}\tilde{\varvec{M}}_{xx}\varvec{u}_x^{(m)}\varvec{u}_y^T\tilde{\varvec{K}}_{yy}\varvec{u}_y^{(m)}\varvec{u}_z^T\tilde{\varvec{M}}_{zz}\varvec{u}_z^{(m)}\nonumber \\{} & {} \quad -\sum _{m=1}^{M-1}\tilde{\varvec{M}}_{xx}\varvec{u}_x^{(m)}\varvec{u}_y^T\tilde{\varvec{M}}_{yy}\varvec{u}_y^{(m)}\varvec{u}_z^T\tilde{\varvec{K}}_{zz}\varvec{u}_z^{(m)}\nonumber \\ \end{aligned}$$

(93)

By solving the above equation, an estimate of $\varvec{u}_x$ and therefore $u_x$ can be obtained. Similarly, we can solve for $u_y$ and $u_z$. This alternative fixed point procedure should be repeated until the convergence the product of $u_x u_y u_z$. This is the solution for one C-PGD mode and can be used for computing incrementally all the modes by varying M from 1 to a given number. The total number of modes can be determined by the convergence criterion: $\Vert \sum _{m=1}^{M-1} u^{(m)}_x(x) u^{(m)}_y(y) u^{(m)}_z(z) -u^{(M)}_x(x) u^{(M)}_y(y) u^{(M)}_z(z)\Vert$ is small enough.

By adopting the TD definition in [39], the C-TD can use the same solution procedure to give a rough estimate number of modes and then optimize all the current modes together. The advantages in doing so is to reduce the necessary number of modes. In our work, we are interested in the final convergent accuracy of C-PGD/TD with comparison to traditional FEM and PGD/TD. The number of modes is ignored. Hence, we do not distinguish the C-PGD and C-TD in this paper as they both can give a better accuracy than traditional methods with the proposed convolution approximation.

Appendix G Error bound analysis of different methods

Assuming a convex potential energy $\varPi$ exits for a problem, the numerical solution to the problem can be obtained by solving the following minimization problem

$$\begin{aligned} \displaystyle u^{h}=\underset{u^{h*}\in V^h }{\arg \min }\ \varPi (u^{h*}(\varvec{x})) \end{aligned}$$

(94)

where $V^h$ denotes a given approximation space. Now, consider two approximation spaces $V^h_1$, $V^h_2$ with $V^h_1 \subset V^h_2$. Their corresponding solutions by solving the above problem are respectively $u^h_1$ and $u^h_2$. It is easy to know that the optimized potential energy has following relationship

$$\begin{aligned} \varPi (u^{\text {Ext}})\le \varPi (u^h_2)\le \varPi (u^h_1) \end{aligned}$$

(95)

where $u^{\text {Ext}}$ is the exact solution. By the convexity of the problem, we have

$$\begin{aligned} \Vert u^h_2-u^{\text {Ext}})\Vert _E \le \Vert u^h_1-u^{\text {Ext}})\Vert _E \end{aligned}$$

(96)

By analogy, denoting the approximation spaces of TD with infinite modes, C-TD with infinite modes, FEM, C-FEM, C-HiDeNN, DeNN respectively by $V^{\text {TD}}$, $V^{\text {C-TD}}$, $V^{\text {FEM}}$, $V^{\text {C-FEM}}$, $V^{\text {C-FEM}}$, $V^{\text {C-HiDeNN}}$, $V^{\text {DeNN}}$, with the following relationship

$$\begin{aligned} V^{\text {TD}}\subset V^{\text {FEM}} \subset V^{\text {C-TD}} \subset V^{\text {C-FEM}} \subset V^{\text {C-HiDeNN}} \subset V^{\text {DeNN}} \end{aligned}$$

(97)

The optimized potential energy follows

$$\begin{aligned} \varPi (u^{\text {Ext}}){} & {} \le \varPi (u^{\text {DeNN}}) \nonumber \\{} & {} \le \varPi (u^{\text {C-HiDeNN}}) \nonumber \\{} & {} \le \varPi (u^{\text {C-FEM}}) \le \varPi (u^{\text {C-TD}})\nonumber \\{} & {} \le \varPi (u^{\text {FEM}}) \le \varPi (u^{\text {TD}}) \end{aligned}$$

(98)

Therefore, we have

$$\begin{aligned} \Vert u^{\text {DeNN}}-u^{\text {Ext}}\Vert _E{} & {} \le \Vert u^{\text {C-HiDeNN}}-u^{\text {Ext}}\Vert _E \nonumber \\{} & {} \le \Vert u^{\text {C-FEM}}-u^{\text {Ext}}\Vert _E \nonumber \\{} & {} \le \Vert u^{\text {C-PGD/TD}}-u^{\text {Ext}}\Vert _E \nonumber \\{} & {} \le \Vert u^{\text {FEM}}-u^{\text {Ext}}\Vert _E \nonumber \\{} & {} \le \Vert u^{\text {PGD/TD}}-u^{\text {Ext}}\Vert _E \end{aligned}$$

(99)

Here, since we are only interested in the approximation space of TD/PGD, we do not distinguish the two methods in this analysis.