[jakubstary]

Hello,

I compute eikonal equation, which gives me a matrix of potential on the domain and I need to compute gradient of this potential. I use following approximation:

grad_i = 1/|C_i| * sum((phi_i + phi_j)/2 * normal_ij),

where C_i is finite volume, sum is over all neighbors of cell i, phi_i is potential in cell i and normal_ij is outward unit normal to surface between cell i and cell j. My implementation works fine on rectangular grid, but on triangular mesh it gives me bad results. Don't you know, where could be mistake on triangular grid?

My implementation in Matlab for triangular grid, where variable v are finite volumes, p are points and phi is potential:

for i=1:num_v
    neighbors = pdeent(v,i); %neighbors 
    neighbors(neighbours == i) = []; %remove cell i from neighbors 
    num_nbs = size(neighbours,2); %neighbors count
    if(num_nbs ~= 3) %on boundary zero
        continue;
    end
    
    for j=1:num_nbs
        nb_v = neighbours(j); %adjacent cell id
        edge_points = intersect(v(1:3,i),v(1:3,nb_v)); %edge points
        
        %edge size
        x = p(1,edge_points(1)) - p(1,edge_points(2));
        y = p(2,edge_points(1)) - p(2,edge_points(2));
        edges_h(1,j) = sqrt(x^2 + y^2);
        
        %computation of outward unit normal
        p3_i = v(1:3,i); p3_i(p3_i==edge_points(1)) = []; p3_i(p3_i==edge_points(2)) = [];
        normal = -flipud(p(:,edge_points(2)) - p(:,edge_points(1))).*[-1;1]; %normal
        e2_i =  p(:,p3_i) - p(:,edge_points(1)); %vector for control of outward
        if(dot(normal,e2_i) > 0) %if not outward
            normal = -normal; %switch
        end
        unit_normal(1:2,j) = normal/norm(normal); %unit normal

        average(j,1) = (phi(i) + phi(nb_v))/2;
    end
    
    gradient(1:2,i) = (unit_normal(1:2,:)*average)/max(edges_h);
end

Thank you.
Jakub

[mprinkey]

Your formulation is wrong for a Gauss gradient calculation. It's likely that the regular mesh that you used before hid the error.

The Gauss gradient formula turns the gradient averaged over the cell (a volume integral) into a surface integral of the function value over the boundary. Dividing by the cell volume is correct--that part comes from the volume integral. But the summation in the formula is approximating the surface integral. That means that the normal_ij in your formula is NOT a unit vector. It is the face area vector and should have a magnitude that is equal to the face's physical area.

I haven't debugged your code, but if you use the face area vector in place of the unit normal vectors, the formulation should be correct. Good luck.

PS. Have a look here: http://www.cfd-online.com/Wiki/Gradient_computation

Also, you are using a pretty sloppy interpolation for the face value in just taking the average of the two. Generally, it is better to at least use linear interpolation. In fact, even linear interpolation without correction may be inadequate if your mesh is skewed.