# Recovery of smoothed continuous stresses using the Z^2 error estimator

I have trouble implementing the $$Z^2$$ error estimator (The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique by Zienkiewicz and Zhu). For this, I am working with Abaqus 6.14 and use a remeshing scheme to optimize the mesh density and solution quality.

I have to minimize the minimization condition $$F(\boldsymbol{a})=\sum_{i=1}^n\,\left(\sigma_h(x_i,y_i)-\sigma_p^*(x_i,y_i)\right)^2$$ where $$\boldsymbol{\sigma_h}$$ represents the stresses that I have at my integration points with the coordinates $$(x_i, y_i)$$. I can easily get them from my output file. (Please note that the stresses in the summation are scalar components of the respective stress tensors.)

I have a problem with the smoothed continuous stresses $$\boldsymbol{\sigma^*_p}$$ that are given by $$\sigma_p^*=\boldsymbol{P}(x_i,y_i)\,\boldsymbol{a}$$.

$$\boldsymbol{P}(x_i,y_i)$$ represents a set of polynomial terms which is predefined and $$\boldsymbol{a}$$ is a set of unknown parameters. I need $$\boldsymbol{a}$$ to calculate my smoothed continuous stresses. In my case, I use triangular elements with one integration point. Therefore, (I think that) I only need a linear expansion and so I have set $$\boldsymbol{P}=[1,x,y]$$ and $$\boldsymbol{a}=[a_1, a_2, a_3]^T$$, where $$(x_i,y_i)$$ are the coordinates of the current integration point.

Since I want to minimize $$F(\boldsymbol{a})$$, it is valid to write $$\sum_{i=1}^n\,\sigma_h(x_i,y_i)=\sum_{i=1}^n\,\sigma_p^*(x_i,y_i)$$.

After some fiddeling, the final equation to calculate $$\boldsymbol{a}$$ takes the form $$\boldsymbol{a}=\boldsymbol{A}^{-1}\,\boldsymbol{b}$$, where $$\boldsymbol{A}=\sum_{i=1}^n\,\boldsymbol{P^T}(x_i,y_i)\boldsymbol{P}(x_i,y_i)$$ and $$\boldsymbol{b}=\sum_{i=1}^n\,\boldsymbol{P^T}(x_i,y_i)\sigma_h(x_i,y_i)$$.

My problem lies within the fact that I cannot invert $$\boldsymbol{A}$$, since $$\det{\boldsymbol{(A)}}=0$$.

I have tried multiple ways to formulate $$\boldsymbol{P}$$, like $$\boldsymbol{P}=[1,x,y]$$, $$\boldsymbol{P^T}=[1,x,y]^T$$ and $$\boldsymbol{a}=[a_1, a_2, a_3]^T$$ and

$$\boldsymbol{P}=\begin{bmatrix}1&x&y\\1&x&y\\1&x&y\end{bmatrix}$$, $$\boldsymbol{P^T}=\begin{bmatrix}1&1&1\\x&x&x\\y&y&y\end{bmatrix}$$ and $$\boldsymbol{a}=\begin{bmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{bmatrix}$$.

I know that I am overlooking something, but I cannot say what it is. I really hope that I could precisely show where my problem lies and that somebody can help me with this.

I have posted a previous question that maybe offers a little bit more context to this question, if needed or interested.

Define the $$n$$ x $$3$$ matrix $$A$$ where each row $$A_i=<1,x_i,y_i>$$, where $$n$$ is the number of evaluation points in your patch. In my version, $$x_i$$ and $$y_i$$ are the absolute coordinates, not relative to the central node as in the book. Presumably, this is what the authors mean by $$P(x_i,y_i)$$. Then, define the $$n$$ x $$6$$ matrix of right-hand sides $$B$$ where each row $$B_i = σ_i$$, your "raw" $$6$$ x $$1$$ Voigt stress vectors calculated by FEM.
Now you have the over-determined system $$AX=B$$, where $$X$$ is a $$3$$ x $$6$$ matrix. Each of the columns of $$X$$ holds the basis function of a component of $$σ^∗$$.
Obviously, you can't invert $$A$$ as it is an $$n$$ x $$3$$ matrix (ignoring the case where $$n$$ = 3), but you can solve this system in a least-squares sense. If you aren't familiar with solving these systems directly with SVD or QR decomposition, the following works as well (this appears to be what the authors were trying to convey, albeit in an inefficient and convoluted manner): $$X = (A^TA)^{-1}A^TB$$ where $$σ^∗_k=X^TP(x_k,y_k)$$ As noted in the book, you should not follow this procedure for nodes located on boundaries and interfaces and use the basis functions generated for an adjacent patch and use some sort of averaging when the boundary node belongs to multiple patches.