5
$\begingroup$

I am in the process of implementing adaptive mesh refinement for a finite element code that solves the Poisson equation. I have had some trouble finding good references on deciding which elements to refine. I have come across the Kelly error estimator from the deal.ii library:

\begin{equation} \eta^{2} = \Sigma_{F\epsilon{}\partial{K}}C_{F}\int_{\partial{K}_{F}}(\nabla{u_{h}}\dot{}n)^{2}do \end{equation}

I have two questions about using this error estimator:

1) My first question is about how to solve this integral. Let me run through what I am currently doing. Consider the triangular element:

enter image description here

where $u_{1}$, $u_{2}$, and $u_{3}$ is the approximate solution at each node, i.e. $u_{h} = [u_{1},u_{2},u_{3}]$. Now to solve this integral over the first face $F_{1}$, i.e.:

enter image description here

I calculate $(\nabla{u_{h}}\dot{}n)^{2}$ as:

\begin{equation} (\nabla{u_{h}}\dot{}n)^{2} = [(u_{1}\frac{\partial{N_{1}}}{\partial{x}} + u_{2}\frac{\partial{N_{2}}}{\partial{x}} + u_{3}\frac{\partial{N_{3}}}{\partial{x}})n_{x} + (u_{1}\frac{\partial{N_{1}}}{\partial{y}} + u_{2}\frac{\partial{N_{2}}}{\partial{y}} + u_{3}\frac{\partial{N_{3}}}{\partial{y}})n_{y}]^{2} \end{equation}

where for the linear triangle $\frac{\partial{N_{1}}}{\partial{x}}$, $\frac{\partial{N_{2}}}{\partial{x}}$, $\frac{\partial{N_{3}}}{\partial{x}}$, $\frac{\partial{N_{1}}}{\partial{y}}$... are all constant.

Thus to calculate $\eta$ I first parameterize face 1 using:

\begin{equation} x = x_{1}(1-s) + x_{2}s\\ y = y_{1}(1-s) + y_{2}s \end{equation}

I can then calculate the line integral as (for face $F_{1}$):

\begin{equation} \int_{0}^{1}(\nabla{u_{h}}\dot{}n)^{2}\sqrt{(dx/ds)^{2}+(dy/ds)^{2}}ds =(\nabla{u_{h}}\dot{}n)^{2}\sqrt{(x_{2}-x_{1})^{2}+(y_{2}-y_{1})^{2}} \end{equation}

where again for linear triangles $(\nabla{u_{h}}\dot{}n)^{2}$ is constant. My first question is whether this how we correctly calculate this line integral on each face?

2) Assuming that I am indeed calculating $\eta^{2}$ correctly, how do we use this to determine which elements should be refined? Do you find the maximum $\eta^{2}$ and then just refine only elements that are a certain percentage of the maximum?

$\endgroup$
  • 1
    $\begingroup$ The formula for the error estimator (and others below) is wrong in that it is not the square of the normal derivative you should be integrating, but the square of the jump of the normal derivative across cell faces. $\endgroup$ – Wolfgang Bangerth Dec 7 '15 at 0:56
2
$\begingroup$

As mentioned in a comment below the question, the formula you state is incorrect: it should read

\begin{equation} \eta^{2} = \Sigma_{F\epsilon{}\partial{K}}C_{F}\int_{\partial{K}_{F}}([\nabla{u_{h}]}\dot{}n)^{2}do \end{equation} where $[\cdot]$ is the jump of the quantity across the face.

Now for your question:

  • Because for linear elements, the gradient is constant. Consequently, if you integrate over a single face $f \subset \partial K_F$, you have that \begin{equation} \int_{f}([\nabla{u_{h}]}\dot{}n)^{2}do = ([\nabla{u_{h}]}\dot{}n)^{2} \int_{f}do = ([\nabla{u_{h}]}\dot{}n)^{2} |f|. \end{equation}

  • You use the error indicators computed for all cells to identify, for example, the 30% of cells with the largest indicators and then refine these. Alternatively, you can identify those cells with the largest indicators that together amount for 90% of the total error. Take a look, for example, at the documentation of the functions here: https://www.dealii.org/developer/doxygen/deal.II/namespaceGridRefinement.html#ad3b68e645838ebeb4f9c55352b56a0b3

$\endgroup$
  • $\begingroup$ Thanks Prof Bangerth. I cant believe I read that deal.ii page so many times and did not notice that it should be the jump across the face, which makes sense. One last question: In the case of triangles what should the $h_{F}$ in $C_{F}=h_{F}/2p_{F}$ be? Would just setting $h_{F}$ equal to the face length suffice? $\endgroup$ – James Dec 15 '15 at 5:47
  • $\begingroup$ Yes, for triangles you can just use the edge length. On quasi-uniform meshes, all measures of edge length and cell diameter are equivalent up to a constant, so you can also take the cell diameter or square root of the cell area if you want. $\endgroup$ – Wolfgang Bangerth Dec 18 '15 at 6:17
  • $\begingroup$ Since you had trouble with that page, could you suggest ways in which it could be improved? $\endgroup$ – Wolfgang Bangerth Dec 18 '15 at 6:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.