dgSoluce.py

# -------------------------------------------------------------------------
#
# PYTHON for FEM DUMMIES 18-19
# Problème 5
#
# Solution détaillée du calcul de l'équation de Poisson
# avec des triangles linéaires continus dits de Turner
#
#  Vincent Legat
#
# -------------------------------------------------------------------------
# 

import numpy as np
import glfem 
import fem

# -------------------------------------------------------------------------
#
# -1- Calcul du coefficient zeta (appelé aussi facteur IP)
#     
#     Référence : Shahbazi, K., 2005. 
#        An explicit expression for the penalty parameter of the interior penalty method. 
#        Journal of Computational Physics 205, 401–407. doi:10.1016/j.jcp.2004.11.017.
#
#     Pas mal d'étudiants ont trouvé l'expression :-)
#     Un peu trop facile ?
#
# -------------------------------------------------------------------------


def computePenaltyFactor(h):
  zeta = (4.5) /h
  return zeta
  
    
def computeGeneralizedPenaltyFactor(p,dim,nSide,h) :
  zeta= (p + 1) * (p + dim) * nSide  / (2 * dim * h)
  return zeta
  
# -------------------------------------------------------------------------
#
# -2- Solution analytique (fournie dans le canevas)
#
# -------------------------------------------------------------------------

def analytic(x,y) :
  Uref = np.zeros_like(x)
  imax = 200;
  for i in range(1,imax,2):
    for j in range(1,imax,2):
      Uref += (1/((i*i+j*j)*i*j)) * np.sin(i*np.pi*(x+1)/2) * np.sin(j*np.pi*(y+1)/2)		
  Uref *= 64 / np.pi**4
  return Uref

# -------------------------------------------------------------------------
#
# -3- Numérotation des inconnues :-)
#
# -------------------------------------------------------------------------
    
def mapTriangle(theProblem,iElem) :
  return [3*iElem+j for j in range(3)]
  
  
def mapEdge(theProblem,iEdge) :
  theMesh  = theProblem.mesh
  theEdges = theProblem.edges
  myEdge = theEdges.edges[iEdge]
  elementLeft  = myEdge[2]
  nodesLeft    = theMesh.elem[elementLeft]
  mapEdgeLeft  = [3*elementLeft + np.nonzero(nodesLeft == myEdge[j])[0][0]  for j in range(2)]
  mapLeft      = [3*elementLeft + j                                         for j in range(3)]
  elementRight = myEdge[3]
  if (elementRight != -1) :
    nodesRight   = theMesh.elem[elementRight]
    mapEdgeRight = [3*elementRight + np.nonzero(nodesRight == myEdge[j])[0][0] for j in range(2)]
    mapRight     = [3*elementRight + j                                         for j in range(3)]
  else :
    mapEdgeRight = []
    mapRight     = []
  return [mapEdgeLeft,mapLeft,mapEdgeRight,mapRight]


# -------------------------------------------------------------------------
#
# -4- Assemblage de la matrice
#     en utilisant quelques fonctions intermédiaires :-)
#
# -------------------------------------------------------------------------


def computeShapeTriangle(theMesh,theElement) :
  dphidxsi = np.array([ 1.0, 0.0,-1.0])
  dphideta = np.array([ 0.0, 1.0,-1.0])

  nodes = theMesh.elem[theElement]
  x = theMesh.X[nodes]
  y = theMesh.Y[nodes]
  dxdxsi = x @ dphidxsi
  dxdeta = x @ dphideta
  dydxsi = y @ dphidxsi
  dydeta = y @ dphideta
  jac = abs(dxdxsi*dydeta - dxdeta*dydxsi)
  dphidx = (dphidxsi * dydeta - dphideta * dydxsi) / jac
  dphidy = (dphideta * dxdxsi - dphidxsi * dxdeta) / jac
  return [dphidx,dphidy,jac]

# -------------------------------------------------------------------------

  def computeShapeEdge(theEdges,iEdge): 
  phiphi   = np.array([[ 2.0, 1.0],[ 1.0, 2.0]]) / 3.0

  nodes = theEdges.edges[iEdge][0:2]
  x = theEdges.mesh.X[nodes]   
  y = theEdges.mesh.Y[nodes]   
  dx = x[1] - x[0]
  dy = y[1] - y[0]
  jac = np.sqrt(dx*dx+dy*dy)
  nx =  dy / jac
  ny = -dx / jac
  return[phiphi,nx,ny,jac]

# -------------------------------------------------------------------------

  
def assemble(theProblem):
  theMesh  = theProblem.mesh
  theEdges = theProblem.edges
  A        = theProblem.system.A
  B        = theProblem.system.B
 
#
# -4.1- Assemblage des intégrales de surface
#
  
  for iElem in range(theMesh.nElem) :
    [dphidx,dphidy,jac] = computeShapeTriangle(theMesh,iElem)    
    map = mapTriangle(theProblem,iElem)
    for i in range(3):
      for j in range(3):
        A[map[i],map[j]] += (dphidx[i] * dphidx[j] + dphidy[i] * dphidy[j]) * jac / 2.0
      B[map[i]] += jac / 6.0
                
#
# -4.2- Assemblage des intégrales de ligne frontières
#

  
  for iEdge in range(theEdges.nBoundary):
    [phiphi,nx,ny,jac] = computeShapeEdge(theEdges,iEdge)
    
    myEdge = theEdges.edges[iEdge]
    elementLeft  = myEdge[2]
    nodesLeft    = theMesh.elem[elementLeft]
    mapEdgeLeft  = [3*elementLeft + np.nonzero(nodesLeft == myEdge[j])[0][0]  for j in range(2)]
    mapLeft      = [3*elementLeft + j                                         for j in range(3)]
    [dphidx,dphidy,jacTriangle] = computeShapeTriangle(theMesh,elementLeft)
    dphidnLeft   = nx * dphidx + ny *dphidy      
    zeta         = computePenaltyFactor(jacTriangle/(jac*2)) 
    
    for i in range(2):
      for j in range(2):
        A[mapEdgeLeft[i] ,mapEdgeLeft[j] ] += zeta * phiphi[i,j] * jac / 4.0
      for j in range(3):
        A[mapEdgeLeft[i] ,mapLeft[j]     ] -= dphidnLeft[j] * jac / 2.0
        A[mapLeft[j]     ,mapEdgeLeft[i] ] -= dphidnLeft[j] * jac / 2.0

#
# -4.3- Assemblage des intégrales de ligne intérieures
#


  for iEdge in range(theEdges.nBoundary,theEdges.nEdges):
    [phiphi,nx,ny,jac] = computeShapeEdge(theEdges,iEdge)
    
    myEdge = theEdges.edges[iEdge]

    elementLeft  = myEdge[2]
    nodesLeft    = theMesh.elem[elementLeft]
    mapEdgeLeft  = [3*elementLeft + np.nonzero(nodesLeft == myEdge[j])[0][0]  for j in range(2)]
    mapLeft      = [3*elementLeft + j                                         for j in range(3)]
    [dphidx,dphidy,jacTriangle] = computeShapeTriangle(theMesh,elementLeft)
    dphidnLeft   = nx * dphidx + ny *dphidy      
    zeta         = computePenaltyFactor(jacTriangle/(jac*2)) 
 
    elementRight = myEdge[3]
    nodesRight   = theMesh.elem[elementRight]
    mapEdgeRight = [3*elementRight + np.nonzero(nodesRight == myEdge[j])[0][0] for j in range(2)]
    mapRight     = [3*elementRight + j                                         for j in range(3)]
    [dphidx,dphidy,jacTriangle] = computeShapeTriangle(theMesh,elementRight)
    dphidnRight  = nx * dphidx + ny *dphidy      
    zeta         = max(zeta,computePenaltyFactor(jacTriangle/(jac*2)))
          
    for i in range(2):
      for j in range(2):
        A[mapEdgeLeft[i] ,mapEdgeLeft[j] ] += zeta * phiphi[i,j] * jac / 4.0
        A[mapEdgeLeft[i] ,mapEdgeRight[j]] -= zeta * phiphi[i,j] * jac / 4.0 
        A[mapEdgeRight[i],mapEdgeLeft[j] ] -= zeta * phiphi[i,j] * jac / 4.0 
        A[mapEdgeRight[i],mapEdgeRight[j]] += zeta * phiphi[i,j] * jac / 4.0 
      for j in range(3):
        A[mapEdgeLeft[i] ,mapLeft[j]     ] -= dphidnLeft[j]  * jac / 4.0
        A[mapEdgeLeft[i] ,mapRight[j]    ] -= dphidnRight[j] * jac / 4.0
        A[mapEdgeRight[i],mapLeft[j]     ] += dphidnLeft[j]  * jac / 4.0
        A[mapEdgeRight[i],mapRight[j]    ] += dphidnRight[j] * jac / 4.0          
        A[mapLeft[j]     ,mapEdgeLeft[i] ] -= dphidnLeft[j]  * jac / 4.0
        A[mapLeft[j]     ,mapEdgeRight[i]] += dphidnLeft[j]  * jac / 4.0
        A[mapRight[j]    ,mapEdgeLeft[i] ] -= dphidnRight[j] * jac / 4.0
        A[mapRight[j]    ,mapEdgeRight[i]] += dphidnRight[j] * jac / 4.0



# -------------------------------------------------------------------------