true

xxxxxxxxxx
 
begin 
    ENV["LC_NUMERIC"]="C"
    using Pkg
    Pkg.activate(mktempdir())
    Pkg.add(["PyPlot","PlutoUI","Triangulate"])
    using PlutoUI,PyPlot, Triangulate,SparseArrays, Printf
    PyPlot.svg(true);
end

5.3 s

Contents

Implementation of the finite volume method

Geometrical data for finite volumes

Needed data

Calculation steps for the interface contributions

Steps to the implenentation

Triangle form factors

Boundary form factors

Matrix assembly

Graphical representation

Calculation example

Grid generation

Plotting the grid

Desired number of triangles

Solving the problem

Problem data

7.4 ms

Implementation of the finite volume method

Here, we specifically introduce the Voronoi finite volume method on triangular grids.

We discuss the implementation of the method for the problem

$\begin{aligned} - \nabla δ \cdot \nabla u & = f \\ δ \partial_{n} u + α u & = g \end{aligned}$

4.8 μs

Geometrical data for finite volumes

As seen in the previous lecture, we need to be able to calculate the contributions to the Voronoi cell data for each triangle.

3.2 μs

PA=[3, 3]

95.3 ms

995 ms

Needed data

Edge lengths $h_{k l}$ :

$a = | P_{B} P_{C} |, b = | P_{A} P_{C} |, c = | P_{A} P_{B} |$

Contributions to lengths of the interfaces between Voronoi cells $| σ_{k l} \cap T |$ – $s_{a}, s_{b}, s_{c}$ : length fo lines joining the corresponding edge centers $P_{B C}, P_{A C}, P_{A B}$ with the triangle circumcenter $P_{C C}$ .
Practically, we need the values of the ratios $\frac{σ_{k l}}{h_{k l}}$ :

$e_{a} = \frac{s_{a}}{a}, e_{b} = \frac{s_{b}}{b}, e_{c} = \frac{s_{c}}{c}$

Triangle contributions to the Voronoi cell areas around the respective triangle nodes $ω_{A} = | P_{A} P_{A B} P_{C C} P_{A C} |, ω_{B} = | P_{B} P_{B C} P_{C C} P_{A B} |, ω_{C} = | P_{C} P_{A C} P_{C C} P_{B C} |$

8.7 μs

Calculation steps for the interface contributions

We show the calculation steps for $e_{a}, ω_{a}$ , the others can be obtained via corresponding permutations.

Semiperimeter:

$s = \frac{a}{2} + \frac{b}{2} + \frac{c}{2}$

Square area (from Heron's formula):

$16 A^{2} = 16 s (s - a) (s - b) (s - c) = (- a + b + c) (a - b + c) (a + b - c) (a + b + c)$

Square circumradius:

$R^{2} = \frac{a^{2} b^{2} c^{2}}{(- a + b + c) (a - b + c) (a + b - c) (a + b + c)} = \frac{a^{2} b^{2} c^{2}}{16 A^{2}}$

Square of the Voronoi interface contribution via Pythagoras:

$s_{a}^{2} = R^{2} - {(\frac{1}{2} a)}^{2} = - \frac{a^{2} {(a^{2} - b^{2} - c^{2})}^{2}}{4 (a - b - c) (a - b + c) (a + b - c) (a + b + c)}$

Square of interface contribution over edge length:

$e_{a}^{2} = \frac{s_{a}^{2}}{a^{2}} = - \frac{{(a^{2} - b^{2} - c^{2})}^{2}}{4 (a - b - c) (a - b + c) (a + b - c) (a + b + c)} = \frac{(b^{2} + c^{2} - a^{2})^{2}}{64 A^{2}}$

Interface contribution over edge length:

$e_{a} = \frac{s_{a}}{a} = \frac{b^{2} + c^{2} - a^{2}}{8 A}$

Calculation of the area contributions

$ω_{a} = \frac{1}{4} c s_{c} + \frac{1}{4} b s_{b} = \frac{1}{4} (c^{2} e_{c} + b^{2} e_{b})$

10.2 μs

The sign chosen implies a positive value if the angle $α_{A} < \frac{π}{2}$ , and a negative value if it is obtuse. In the latter case, this corresponds to the negative length of the line between edge midpoint and circumcenter, which is exactly the value which needs to be added to the corresponding amount from the opposite triangle in order to obtain the measure of the Voronoi face.
If an edge between two triangles is not locally Delaunay, the summary contribution from the two triangles with respect to this edge will become negative.

4.4 μs

Steps to the implenentation

We describe a triangular discretization mesh by three arrays:

pointlist: $2 \times n_{p o i n t s}$ floating point array of node coordinates of the triangulations. pointlist[:,i] then contains the coordinates of point i.
trianglelist: $3 \times n_{n t r i}$ integer array describing which three nodes belong to a given triangle. trianglelist[:,i] then contains the numbers of nodes belonging to triangle i.
segmentlist: $2 \times n_{s e g s}$ integer array describing which two nodes belong to a given boundary segment. segmentlist[:,i] contains the numbers of nodes for boundary segment i.

Triangle form factors

For triangle itri, we want to calculate the corresponding form factors $e$ and $ω$ :

7.6 μs

trifactors! (generic function with 1 method)

xxxxxxxxxx
 
function trifactors!(ω, e, itri, pointlist, trianglelist)
    # Obtain the node numbers for triangle itri
    i1=trianglelist[1,itri]
    i2=trianglelist[2,itri]
    i3=trianglelist[3,itri]
    
    # Calculate triangle area: 
    #   Matrix of edge vectors
    V11= pointlist[1,i2]- pointlist[1,i1]
    V21= pointlist[2,i2]- pointlist[2,i1]
    
    V12= pointlist[1,i3]- pointlist[1,i1]
    V22= pointlist[2,i3]- pointlist[2,i1]
    
    V13= pointlist[1,i3]- pointlist[1,i2]
    V23= pointlist[2,i3]- pointlist[2,i2]
    
    #   Compute determinant 
    det=V11*V22 - V12*V21
   
    #   Area
    area=0.5*det
    
    # Squares of edge lengths
    dd1=V13*V13+V23*V23 # l32
    dd2=V12*V12+V22*V22 # l31
    dd3=V11*V11+V21*V21 # l21
        
    # Contributions to e_kl=σ_kl/h_kl
    e[1]= (dd2+dd3-dd1)*0.125/area
    e[2]= (dd3+dd1-dd2)*0.125/area
    e[3]= (dd1+dd2-dd3)*0.125/area
    
    # Contributions to ω_k
    ω[1]= (e[3]*dd3+e[2]*dd2)*0.25
    ω[2]= (e[1]*dd1+e[3]*dd3)*0.25
    ω[3]= (e[2]*dd2+e[1]*dd1)*0.25
end                              

66.6 μs

Boundary form factors

Here we need for an interface segment of two points $P_{A}, P_{B}$ the contributions to the intersection of the Voronoi cell boundary with the outer boundary which is just the half length: $γ_{A} = \frac{1}{2} | P_{A} P_{B} |, γ_{B} = \frac{1}{2} | P_{A} P_{B} |$

3.3 μs

bfacefactors! (generic function with 1 method)

xxxxxxxxxx
 
function bfacefactors!(γ,ibface, pointlist, segmentlist)
    i1=segmentlist[1,ibface]
    i2=segmentlist[2,ibface]
    dx=pointlist[1,i1]-pointlist[1,i2]
    dy=pointlist[2,i1]-pointlist[2,i2]
    d=0.5*sqrt(dx*dx+dy*dy)
    γ[1]=d
    γ[2]=d
end

22.2 μs

Matrix assembly

The matrix assembly consists of two loops, one over all triangles, and another one over the boundary segments.

The implementation hints at the possibility to work in different space dimensions

3.5 μs

assemble! (generic function with 1 method)

xxxxxxxxxx
 
function assemble!(matrix, # System matrix
                   rhs,    # Right hand side vector
                   δ,      # heat conduction coefficient 
                   f::Function, # Source/sink function
                   α,      # boundary transfer coefficient
                   g::Function,  # boundary condition function
                   pointlist,    
                   trianglelist,
                   segmentlist)
   
    num_nodes_per_cell=3;
    num_edges_per_cell=3;
    num_nodes_per_bface=2
    ntri=size(trianglelist,2)
    nbface=size(segmentlist,2)
    
    # Local edge-node connectivity
    local_edgenodes=[ 2 3; 3 1; 1 2]'
   
    # Storage for form factors
    e=zeros(num_nodes_per_cell)
    ω=zeros(num_edges_per_cell)
    γ=zeros(num_nodes_per_bface)
​
    # Initialize right hand side to zero
    rhs.=0.0
​
    # Loop over all triangles
    for itri=1:ntri
        trifactors!(ω,e,itri,pointlist,trianglelist)
​
        # Assemble nodal contributions to right hand side
        for k_local=1:num_nodes_per_cell
            k_global=trianglelist[k_local,itri]
            x=pointlist[1,k_global]
            y=pointlist[2,k_global]
            rhs[k_global]+=f(x,y)*ω[k_local]
        end
        
        # Assemble edge contributions to matrix
        for iedge=1:num_edges_per_cell
            k_global=trianglelist[local_edgenodes[1,iedge],itri]
            l_global=trianglelist[local_edgenodes[2,iedge],itri]
            matrix[k_global,k_global]+=δ*e[iedge]
            matrix[l_global,k_global]-=δ*e[iedge]
            matrix[k_global,l_global]-=δ*e[iedge]
            matrix[l_global,l_global]+=δ*e[iedge]
        end
    end
​
    # Assemble boundary conditions
    
    for ibface=1:nbface
        bfacefactors!(γ,ibface, pointlist, segmentlist)
        for k_local=1:num_nodes_per_bface
            k_global=segmentlist[k_local,ibface]
            matrix[k_global,k_global]+=α*γ[k_local]
            x=pointlist[1,k_global]
            y=pointlist[2,k_global]
            rhs[k_global]+=g(x,y)*γ[k_local]
        end
    end
end
​

66.7 μs

Graphical representation

It would be nice to have a graphical representation of the solution data. We can interpret the solution as a piecewise linear function on the triangulation: each triangle has three nodes each carrying one solution value.

On the other hand, a linear function of two variables is defined by values in three points. This allows to define a piecewise linear, continuous solution function. This approach is well known for the finite element method which we will introduce later.

3.8 μs

plot (generic function with 1 method)

xxxxxxxxxx
 
function plot(u, pointlist, trianglelist)
    cmap="coolwarm" # color map for color coding function values
    num_isolines=10 # number of isolines for plot
    ax=gca(); ax.set_aspect(1) # don't distort the plot
​
    # bring data into format understood by PyPlot
    x=view(pointlist,1,:)
    y=view(pointlist,2,:)
    t=transpose(triout.trianglelist.-1)
​
    # Many (50) filled contour lines give the impression of a smooth color scale
    tricontourf(x,y,t,u,levels=50,cmap=cmap)
    colorbar(shrink=0.5) # Put a color bar next to the plot
    
    # Overlay the plot with isolines
    tricontour(x,y,t,u,levels=num_isolines,colors="k")
end
​

33.3 μs

An alternative way of showing the result is the 3D plot of the function graph:

2.1 μs

plot3d (generic function with 1 method)

xxxxxxxxxx
 
function plot3d(u, pointlist, trianglelist)
    cmap="coolwarm" 
    fig=figure(2)
    x=view(pointlist,1,:)
    y=view(pointlist,2,:)
    t=transpose(triout.trianglelist.-1)
    plot_trisurf(x,y,t,u,cmap=cmap)
end

26.4 μs

Calculation example

Now we are able to solve our intended problem.

Grid generation

3.6 μs

make_grid (generic function with 1 method)

xxxxxxxxxx
 
function make_grid(;maxarea=0.01)
    triin=TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([-1.0 -1.0; 1.0 -1.0 ; 1.0 1.0 ; -1.0 1.0 ]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4])
    a=@sprintf("%f",maxarea)
    (triout, vorout)=triangulate("pqAa$(a)qQD", triin)
    triin, triout
end

69.7 μs

TriangulateIO(
pointlist=[-1.0 1.0 1.0 -1.0; -1.0 -1.0 1.0 1.0],
segmentlist=Int32[1 2 3 4; 2 3 4 1],
segmentmarkerlist=Int32[1, 2, 3, 4],
)

TriangulateIO(
pointlist=[-1.0 1.0 … 0.5 -0.5; -1.0 -1.0 … 1.0 -1.0],
pointmarkerlist=Int32[1, 1, 2, 3, 0, 4, 1, 2, 3, 0  …  4, 1, 2, 0, 0, 3, 4, 2, 3, 1],
trianglelist=Int32[18 13 … 12 15; 12 24 … 3 24; 5 7 … 23 14],
segmentlist=Int32[2 3 … 23 24; 16 22 … 3 1],
segmentmarkerlist=Int32[1, 2, 3, 4, 4, 1, 2, 3, 4, 1, 2, 3, 4, 2, 3, 1],
)

xxxxxxxxxx
 
triin,triout=make_grid(maxarea=4.0/desired_number_of_triangles)

352 ms

Plotting the grid

In the triout data structure, we indeed see a pointlist, a trianglelist and a segmentlist.

We use the plot_in_out function from Triangulate.jl to plot the grid.

Plot grid:

23.1 ms

696 ms

Number of points: 24, number of triangles: 30.

5.7 ms

Desired number of triangles

From the desired number of triangles, we can claculate a value fo the maximum area constraint passed to the mesh generator: Desired number of triangles: 20

6.3 ms

Solving the problem

Problem data

3.5 μs

f (generic function with 1 method)

xxxxxxxxxx
 
f(x,y)=sinpi(x)*cospi(y)

14.9 μs

g (generic function with 1 method)

xxxxxxxxxx
 
 g(x,y)=0

11.7 μs

δ: 1.0 α: 1.0

9.8 ms

solve_example (generic function with 1 method)

xxxxxxxxxx
 
function solve_example(triout)
    n=size(triout.pointlist,2)
    matrix=spzeros(n,n)
    rhs=zeros(n)
    assemble!(matrix,rhs,δ,f,α,g,triout.pointlist,triout.trianglelist, triout.segmentlist)
    sol=matrix\rhs
end 

22.1 μs

solution

Float641

0.0207156

-0.0121475

-0.010301

0.0245238

0.0162066

-0.0152359

0.00557976

0.0377689

0.0104351

0.0121834

0.0184025

0.0169708

0.00189377

-0.0563245

-0.00840472

-0.0464402

0.00999763

0.0832824

0.0155037

0.0643528

0.00921857

0.0122819

-0.0431849

0.0705515

xxxxxxxxxx
 
solution=solve_example(triout)

108 ms

xxxxxxxxxx
 
clf(); plot(solution,triout.pointlist, triout.trianglelist);gcf().set_size_inches(4,4);gcf()

72.7 ms

3D Plot ?

7.1 ms

xxxxxxxxxx
 
if do_3d_plot
    clf(); plot3d(solution,triout.pointlist, triout.trianglelist);gcf().set_size_inches(4,4);gcf()
end

44.0 ns