1
$\begingroup$

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: Mapping

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
$\endgroup$
  • $\begingroup$ There is something wrong with your code. Try a simple case where B(x,y)=1, isosceles triangle: x=[0,1/2,0] ; y=[0,sqrt(3)/2,0] Answer should be 1, but I get 0.5 So you are missing a factor of 2 somewhere. Similarly fi we choose: x=[0,0,1] ; y=[0,1,0] for which the mapping is unitary, we still get 0.5 $\endgroup$ – user21387 Jul 19 '17 at 21:39
  • $\begingroup$ Hi, thank you for your answer :) There was a little typo in the code, I fix it now thanks! I add a similar testing at the end of the post (edit of the 20.07.17), and in fact it works... So I really don't have a clue why it's not working on more complex function $B(x,y)$... $\endgroup$ – Electrode Jul 20 '17 at 8:17
3
$\begingroup$

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]$.

$\endgroup$
  • $\begingroup$ Awesome answer! Of course, I completely switched the fact that the mapping is in fact physical edge -> isoparametrical edge -> 1D segment (I forget this last one). I will accept this answer, since it solved my problem :D In order to be useful for other people, I edited my post (at the end) where I add a working code that compute $\int_{P_{1}}^{P_{2}} 4 x^{3} ds \ = -21.21$ (the result is the correct, as expected) with a 3 node quadrature on $[-1;1]$ on this particular edge (specific implementation). $\endgroup$ – Electrode Jul 20 '17 at 11:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.