Pasquale Africa
Dear deal.II community,
@sdigre and I encountered an issue when computing a finite element solution to a Darcy equation using Raviart-Thomas(0) elements on 3D distorted (non-Cartesian) meshes.
The minimal steps to reproduce the issue (summed up in the code attached) are inspired by the step-20 tutorial, extended trivially to 3D geometries.
We first instantiate the 3D template version of the solver MixedLaplaceProblem<3>: so far everything works, providing the same 2D solution extruded across the third dimension z.
Now if we introduce a slight random mesh distortion, keeping the boundaries fixed, i.e. adding GridTools::distort_random(0.1, triangulation, true) after the mesh creation, then problems arise: the L^2 error on the velocity is about 6 times larger than the one computed on the undistorted mesh, the z component of the velocity is non-zero, whereas the x and y components show an unphysical pattern.
The same behavior happens importing an externally generated 3D cube mesh using any of the the GridIn reader formats.
The 2D version of the code runs smoothly for any kind of mesh.
Trying to debug the code, we noticed that the RT0 basis function associated with a face dof of a distorted cell is not normal to the face itself as it is expected to be, which makes us suspect a bug in the mapping of the Raviart-Thomas functions from the reference cell.
Thanks in advance.
Best regards,
Pasquale and Simone
Pasquale Africa
We have been investigating the error and we suspect that it is related to the issue
#7970
that has been assigned to the milestone "Release 9.2".
In particular I think that it is due to the missing implementation of get_face_sign_change_rt for 3D triangulations
(see https://github.com/dealii/dealii/blob/master/source/fe/fe_poly_tensor.cc#L111).
We will be thinking of a solution to that function, do you have any ideas or suggestions to start from?
Thanks.