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

Question

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.

Answer 1

I just implemented this method as detailed in the book by Zienkiewicz and Zhu (The Finite Element Method - Its Basis and Fundamentals) for linear tetrahedral elements. This is an excellent book, but was totally confused by the notation used in the superconvergent patch recovery section. I understood the procedure as very similar to defining shape functions, except over a patch of elements with a polynomial one degree higher than you used for your original calculations. The following should work for your problem (linear triangles with linear displacement formulation; recovering linear stresses).

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.

What I ended up doing was actually recovering quadratic displacements and derived the linear stresses by differentiating the resulting displacements. This so-called "polynomial preserving" scheme was much more appealing to me, and you end up improving the accuracy (or, at least, appearance) of all possible derived quantities of displacement in a single pass.