Line integral along the edge of an isoparametrically mapped triangle

Question

I need to integrate the following function on the line segment from $P_{1} = \begin{bmatrix} -2\\-1 \end{bmatrix}$ to $P_{2} = \begin{bmatrix} 1\\2 \end{bmatrix}$: $$\int_{P_{1}}^{P_{2}} 4x + y \ ds$$

This question take part into the implementation of a 2D finite element solver.

How I plan to do it

Suppose that the following transformation is used to transform a general triangular element K to the standard triangular element $T_{st}$ :

$$ x = P(\xi,\eta) $$ $$y = Q(\xi,\eta)$$ This corresponds to this physical mapping in fact:

Then we have: $$dx = \frac{\delta x}{\delta \xi}d\xi + \frac{\delta x}{\delta \eta}d\eta = J_{11}d\xi + J_{21}d\eta $$

$$dy = \frac{\delta y}{\delta \xi}d\xi + \frac{\delta y}{\delta \eta}d\eta = J_{12}d\xi + J_{22}d\eta $$

Along the side $P_{1} - P_{2}$, the coordinate $\xi$ is in fact always fixed. So we can write $d\eta = 0$ . Using this it comes:

$$ dx = (\frac{\delta x}{\delta \xi})_{\eta = 0} \ d\xi = J_{11}(\xi,0) d\xi $$ $$ dy = (\frac{\delta y}{\delta \xi})_{\eta = 0} \ d\xi = J_{12}(\xi,0) d\xi $$

Therefore: $$ ds = \sqrt{J_{11}^{2}(\xi,0)+ J_{12}^{2}(\xi,0)} d\xi $$

So we can rewrite the integral with a variable change accordingly: $$\int_{P_{1}}^{P_{2}} 4x + y \ ds = \int_{0}^{1} B(P(\xi,0),Q(\xi,0)) \ \sqrt{J_{11}^{2}(\xi,0)+ J_{12}^{2}(\xi,0)} d\xi $$

Using this isoparametrical formulation, we can finally compute the numerical integral thanks to a 4 node quadrature :

$$ I = \sum_{i=1}^{gauss \ point} w_{i} B(P(\xi_{i},0),Q(\xi_{i},0)) \ \sqrt{J_{11}^{2}(\xi_{i},0)+ J_{12}^{2}(\xi_{i},0)} d\xi $$

This can be done with this really simple code in Python (at least I was thinking):

Edit of the code : 20.07.2017 => removed an error (dividing by two without reason)

# -*- coding: utf-8 -*-
from __future__ import division # avoid integer problem of division

#Import zone
import numpy as np
import math


#Define shape function
def P1shapes(r,s):
    S = np.array([1-r-s,r,s])
    dSdr = np.array([-1,1,0])
    dSds = np.array([-1,0,1])
    return S,dSdr,dSds


#Jacobian function
def isopmap(x,y,r,s,shapefcn):
    # x = vector of x coordinate of the element's point
    #shapefcn = P1shapes
    S,dSdr,dSds = shapefcn(r,s);
    j11=np.dot(dSdr,x)
    j12=np.dot(dSdr,y)
    j21=np.dot(dSds,x)
    j22=np.dot(dSds,y)
    detJ=j11*j22-j12*j21
    dSdx=( j22*dSdr-j12*dSds)/detJ
    dSdy=(-j21*dSdr+j11*dSds)/detJ
    return S,dSdx,dSdy,detJ,j11,j12,j21,j22

#Gauss point quadrature
qwgts=np.array([-27/48,25/48,25/48,25/48])
rspts=np.array([[1/3,1/3],
[0.2,0.2],
[0.6,0.2],
[0.2,0.6]])


def B(x,y):
    z = 4 *x + y
    return z

#Point of the 3 nodes triangles (test case)
x = [-2, 1 ,-1]
y = [-1 ,2 ,3]

#P1-P2
int_total = 0

#Begin integral on segment
for q in range(len(qwgts)):
    r = rspts[q,0] # r coordinate of the q_th quadrature point
    s = rspts[q,1] # s coordinate of the q_th quadrature point
    #Define ds and Map x_physical, y_physical
    S,dSdx,dSdy,detJ,j11,j12,j21,j22 = isopmap(x,y,r,0,P1shapes)
    ds = math.sqrt(math.pow(j11,2) + math.pow(j12,2))

    x_physical = np.dot(x,S)
    y_physical = np.dot(y,S)

    B_gausspoint = B(x_physical,y_physical)
    wxarea = B_gausspoint*qwgts[q]*ds
    int_total += wxarea

print int_total

The result gives: $$ I = -16.97 $$

Of course if I change the third point of the triangle, it should change nothing. With this:

x = [-2, 1 ,-1]
y = [-1 ,2 ,7]

We again find for the path $P_{1} - P_{2}$: $$ I = -16.97 $$

What we should expect

From an analytical point of view we have for the parametrisation $x = -2 + 3t$ and $y = -1 + 3t$. This allow to have:

$$ds = \sqrt{ (\frac{dx}{dt})^{2} + (\frac{dy}{dt})^{2}} \ dt = 3\sqrt{2} \ dt$$

So we have: $$\int_{P_{1}}^{P_{2}} 4 x + y \ ds = \int_{0}^{1} 4(-2+3t) + (-1+3t)) \ 3\sqrt{2} \ dt = -6.36$$

Actually, we can see that the numerical method is not working...

My question

Did I make a mystake in the main concept? Is it an implementation error? I think it could be an interesting question for all the people that try to implement a path integral on a FEM mesh.

I really did'nt find litterature example to check my implementation.

Extension of the question

Is there an other way/more elegant way to compute a path integral on the boundary of my FEM model? It seems that I have everything (non linear solver, quadrature integral for assembly of stiffness matrix, shape function, isoparametric formulation), but I'm really stuck at this point. Without this I never would be able to compute for example:

$$ R_{ij} = \int \kappa(x,y) \ \phi_{i} \ \phi_{j} \ dS \ for \ i \ = \ 1,2 $$

EDIT 20.07.17 : Additional test

Actually, if we take for the integral $B(x,y) = 1$, we are integrating the length of the path :

def B(x,y):
    z = 1
    return z

If we run the code, we end up with : $$ I = 4.24$$

Which is the correct value since the length of the path: $$\sqrt{(P_{2}^{x} -P_{1}^{x})^{2} + (P_{2}^{y} -P_{1}^{y})^{2} } = 4.24$$

Thank you in advance.

Correct code (see explanation on accepted answer)

# -*- coding: utf-8 -*-
from __future__ import division # avoid integer problem of division

#Import zone
import numpy as np
import math


#Define shape function
def P1shapes(r,s):
    S = np.array([1-r-s,r,s])
    dSdr = np.array([-1,1,0])
    dSds = np.array([-1,0,1])
    return S,dSdr,dSds


#Jacobian function
def isopmap(x,y,r,s,shapefcn):
    # x = vector of x coordinate of the element's point
    #shapefcn = P1shapes
    S,dSdr,dSds = shapefcn(r,s);
    j11=np.dot(dSdr,x)
    j12=np.dot(dSdr,y)
    j21=np.dot(dSds,x)
    j22=np.dot(dSds,y)
    detJ=j11*j22-j12*j21
    dSdx=( j22*dSdr-j12*dSds)/detJ
    dSdy=(-j21*dSdr+j11*dSds)/detJ
    return S,dSdx,dSdy,detJ,j11,j12,j21,j22

#Gauss point quadrature    
def transfo1D(a,b,ti):
    #Inverse mapping to find r coordinate from 2D [r,0] space corresponding to 1D space [-1;1] defined by t
    epsylon_i = ((b-a)/2.0)*ti + ((b+a)/2.0)
    return epsylon_i

def B(x,y):
    z = 4*math.pow(x,3)
    return z


coordinate=np.array([-0.774596669,0.000000000,0.774596669])
weight=np.array([0.555555556,0.888888889,0.555555556])


#Point of the 3 nodes triangles (test case)
x = [-2, 1 ,-1]
y = [-1 ,2 ,3]

#P1-P2
int_total = 0

#Begin integral on segment
for q in range(len(coordinate)):
    ti = coordinate[q] # r coordinate of the q_th quadrature point
    #Transform this to r space [r,0]
    a = 0
    b = 1
    r = transfo1D(a,b,ti)
    dettransform = (b-a)/2.0

    #Define ds = fct(r) = fct(r(t))
    S,dSdx,dSdy,detJ,j11,j12,j21,j22 = isopmap(x,y,r,0,P1shapes)
    ds = math.sqrt(math.pow(j11,2) + math.pow(j12,2))

    #Define B(P(r,0),Q(r,0)) = B(P(r(t),0),Q(r(t),0))
    x_physical = np.dot(x,S)
    y_physical = np.dot(y,S)
    B_gausspoint = B(x_physical,y_physical)

    wxarea = weight[q] * B_gausspoint* ds * dettransform
    int_total += wxarea

print int_total

Answer 1

Interestingly enough, you are using quadrature rule for a master triangle in order to integrate over a segment. You should use quadrature rule for a master segment (e.g. $[-1, 1]$) instead. You should also provide proper transformations $\mathbf r_i : [-1, 1] \rightarrow \partial K_i$, $i = 1, 2, 3$ for each edge of a physical triangle.

The best way of doing this is via introducing extra mappings $\mathbf T_i : [-1, 1] \rightarrow \partial \hat K_i$, $i = 1, 2, 3$ for each edge of the master triangle $\hat K$. If $\hat K$ is spanned by vertices $\{(0,0), (1, 0), (0, 1)\}$, then $$ \mathbf T_1(t) = \langle \, .5(1-t), .5(1+t) \, \rangle, \\ \mathbf T_2(t) = \langle \, 0, .5(1-t) \, \rangle, \\ \mathbf T_3(t) = \langle \, .5(1+t), 0 \, \rangle. $$ (I assume that the $i$th edge is the edge against the $i$th vertex.) Then you can set $\mathbf r_i := \mathbf T \circ \mathbf T_i$, where $\mathbf T : \hat K \rightarrow K$ is a usual mapping from the master triangle to the physical one.

So your element Robin matrix can be computed as $$ \mathbf R^K_{ij} := \int_{\partial K_m} \kappa \, s^K_j \, s^K_i \, \text{d}s = \int_{-1}^{1} \left( \kappa \, s^K_j \, s^K_i \right) \circ \mathbf r_m \, || \mathbf r_m' || \, \text{d}t = \\ \int_{-1}^{1} \left( \kappa \circ \mathbf T \circ \mathbf T_m \right) \, \left( \hat s_j \circ \mathbf T_m \right) \, \left( \hat s_i \circ \mathbf T_m \right) \, || (\nabla \mathbf T \circ \mathbf T_m) \, \mathbf T_m' || \, \text{d}t. $$ Note that this choice of transformations $\mathbf r_i$ allows you to use master shape functions (so you can compute their images of quadrature nodes only once and then use them for every physical edge). Note also that the “length element” is constant unless you are using curvilinear mapping $\mathbf T$ (defined e.g. via $P2$ or $P3$ Lagrange shape functions). Finally, note that you have to use quadrature rule for the master segment $[-1, 1]$.