How to use RT03d finite elements?

General Discussion
1

Hey y’all!
As prep for my master’s thesis, I’m trying to implement a numerical method for the heat equation using Raviart-Thomas finite elements. While the documentation has a nice example on how to use RT0 finite elements in 2D, I’m having a lot of trouble trying to get them to work in 3D. I naïvely tried to adapt the code in the LaplaceRT.edp example file as such:

Code 
load "msh3"

// Parameters
func gd = 1.;
func g1n = 1.;
func g2n = 1.;
func g3n = 1.;

// Mesh
mesh3 Th = cube(10, 10, 10);

// Fespace
fespace Vh(Th, RT03d);
Vh [u1, u2, u3];
Vh [v1, v2, v3];

fespace Ph(Th, P0);
Ph p, q;

// Problem
problem laplaceMixte ([u1, u2, u3, p], [v1, v2, v3, q], solver=GMRES, eps=1.0e-10, tgv=1e30, dimKrylov=150)
    = int3d(Th)(
        p*q*1e-15 //this term is here to be sure
        // that all sub matrix are inversible (LU requirement)
        + u1*v1
        + u2*v2
        + u3*v3
        + p*(dx(v1)+dy(v2)+dz(v3))
        + (dx(u1)+dy(u2)+dz(u3))*q
    )
    + int3d(Th) (
        q
    )
    - int2d(Th, 1, 2, 3, 4, 5)(
        gd*(v1*N.x + v2*N.y + v3*N.z)
    )
    + on(6, u1=g1n, u2=g2n, u3=g3n)
    ;

// Solve
laplaceMixte;

// Plot
plot([u1, u2, u3], coef=0.1, wait=true, value=true);
plot(p, fill=1, wait=true, value=true);

The GMRES solver is not converging, and even if I cut it off prematurely, the plots are empty. What’s going wrong here? And advice would be greatly appreciated. I’m using FreeFEM++ version 4.14.

2

Why are you using GMRES?

3

I only changed the parts of the example code that were absolutely necessary to go from 2D to 3D. Using the default solver options (i.e. passing no additional arguments to the problem constructor) still leads to empty plots.

4

Even with solver = sparsesolver?

5

Yup, still empty plots.

6

If you do cout << u1[].min << " " << u1[].max << endl; you’ll see that the solution is clearly nonzero. FreeFEM is not a very suitable tool for 3D visualization, so I would just output the solution using savevtk and do any post-processing in ParaView.

1 Like