Issues with hanging nodes

[Matteo Menessini]

Hello there,

I am working on a code to solve the homogeneous Poisson problem in 2D to numerically estimate the potential and approximate the electric field between two electrodes with non-standard geometry (at the moment, one electrode is circular and the other is a flat plate, but the aim is to extend it to more complex geometries).

Homogeneous Dirichlet conditions are imposed at both electrodes with homogeneous Neumann conditions imposed at the rest of the boundary. The potential is then determined as the sum of the solution of the problem with the aforementioned conditions (named U_h), plus a test function (named R_g) which is equal to a potential value V_max at the edge of the circular electrode and goes to zero at a given distance from the electrode decreasing with an inverse square law. The potential U thus obtained (U = U_h + R_g) is smooth, respects the boundary conditions and seems physical. The electric field, obtained in post-processing using the DataPostProcessorVector class, is also physical: discontinuous at cell edges but decreasing monotonically as the radial distance from the circular electrode increases.

However, upon adaptive grid refinement, while the potential remains smooth, the electric field has presents non-physical spikes in its magnitude. The spike perfectly corresponds to the first layer of cells that have been refined (at any refinement cycle, a new spike appears where there is an increse in cell level), indicating that it is an issue related to how the gradient estimator in DataPostProcessorVector interacts with hanging nodes. See attached images for more details (the pictures show the field near the circular electrode and the graph of the potential and electric field magnitude in the radial direction from the electrode center up to 3R, where R is the circle radius).

In the code, hanging nodes:

- are created during setup and after each mesh refinement

- are condensed in the system matrix and rhs used to obtain U_h with the function distribute_local_to_global during assembly

- are distributed to U_h after solving the system

- after the potential is obtained and U = U_h + R_g is evaluated, constraints are distributed again to U

The homogeneous Dirichlet constraints are applied to the system matrix and rhs after the assembly using MatrixTools:apply_boundary_values().

The function R_g is evaluated in the numerical domain using VectorTools::interpolate(). It is noted that also distributing constaints to R_g does not appear change the results, as does using VectorTools::project() rather than VectorTools::interpolate().

The gradient of the function U_h obtained through DataPostProcessorVector does not present the spikes even after grid refinement, nor does the gradient of R_g. However, their sum U = U_h + R_g develops this issue.

I have also tried "manually" approximating the gradient looping on all active cells using fe_values.get_function_gradients() and fe_face_values.get_function_gradients() on the boundary but the issue remains the same.

Any ideas on what I might be doing wrong? Feel free to request any additional information that you might require, I am new to the forum and quite new to coding and C++, hence while I tried to be as thorough as possible with my explanation I surely missed some important details.

Thanks in advance.

[Wolfgang Bangerth]

Matteo,
you lose one order of convergence between looking at the value of the solution
and looking at the gradient. So I'm not surprised that it isn't very accurate.
In your case, the solution is smooth -- so I would suggest using second or
third (polynomial) order elements. I bet you will get substantially more
accurate answers this way.

I do not believe that what you see is a "bug". You constrain the solution at
hanging nodes, which leaves the space at adjacent cells too small to provide
much accuracy. *Asymptotically* this doesn't matter, but if your mesh is still
too coarse, it does, in particular when looking at derivatives of the solution.

As a separate question, since your domain appears to be circularly symmetric,
are you using an appropriate manifold description in the triangulation, and an
appropriate mapping in DataOut?

Best
W.

On 11/17/23 03:13, Matteo Menessini wrote:
> *** Caution: EXTERNAL Sender ***

> --
> The deal.II project is located at http://www.dealii.org/
> <http://www.dealii.org/>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <https://groups.google.com/d/forum/dealii?hl=en>
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com
> <mailto:dealii+un...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/b077f049-a592-4f16-bd6e-b6f35649dfb7n%40googlegroups.com <https://groups.google.com/d/msgid/dealii/b077f049-a592-4f16-bd6e-b6f35649dfb7n%40googlegroups.com?utm_medium=email&utm_source=footer>.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/