meteoriteSoluce.py

# -------------------------------------------------------------------------
#
# PYTHON for FEM DUMMIES 18-19
# Problème 6
#
# Solution détaillée du calcul du problème de la météorite
# avec des triangles linéaires continus dits de Turner
#
#  Vincent Legat
#
# -------------------------------------------------------------------------
# 

import numpy as np
import glfem 
import fem


# -------------------------------------------------------------------------
#
# -1- Calcul de la contribution de la météorite sur un élément
#     On ne tient pas compte du cas où la météorite sur un sommet et/ou
#     un segment : dans un tel cas, il faudrait adapter le calcul :-)
#
#     Lorsque j'aurai le temps, je ferai une version plus générale pour traiter
#     ce cas particulier dans "assemble" pas dans cette fonction
#
# -------------------------------------------------------------------------


def computeSource(theProblem,theElement) :
  theMesh = theProblem.mesh
  xloc = theProblem.xSource[0]
  yloc = theProblem.xSource[1]
  if (theElement >= theMesh.nElem) :
    return np.zeros(3)

  nodes = theMesh.elem[theElement]
  x = theMesh.X[nodes]
  y = theMesh.Y[nodes]
  jac  = (x[0]-x[1]) * (y[0]-y[2]) - (x[0]-x[2]) * (y[0]-y[1]) 
  phi = np.zeros(3)
  phi[0] = ((xloc-x[1]) * (yloc-y[2]) - (xloc-x[2]) * (yloc-y[1]))/jac  
  phi[1] = ((x[0]-xloc) * (y[0]-y[2]) - (x[0]-x[2]) * (y[0]-yloc))/jac  
  phi[2] = ((x[0]-x[1]) * (y[0]-yloc) - (x[0]-xloc) * (y[0]-y[1]))/jac
  if (phi[0] >= 0.0 and phi[1] >= 0.0 and phi[2] >= 0.0 ) : 
    print(" == %d : %e %e %e %e " % (theElement,jac,*phi[:]))
    return phi
  else :
    return np.zeros(3)

  
def computePenaltyFactor(jac,h):
  return 9.0 * (h/jac)
  

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]


# -------------------------------------------------------------------------
#
# -1- Calcul du problème discret avec les termes de transport :-)
#     C'est une simple adaptation des termes IP en tenant compte du signe
#     produit scalaire beta dot n :-)
#
# -------------------------------------------------------------------------


def assemble(theProblem):
  theMesh  = theProblem.mesh
  theEdges = theProblem.edges
  A        = theProblem.system.A
  B        = theProblem.system.B
  betax     = theProblem.beta[0]
  betay     = theProblem.beta[1]
 
#
# -1- Assemblage des intégrales de surface
#
  
  for iElem in range(theMesh.nElem) :
    [dphidx,dphidy,jac] = computeShapeTriangle(theMesh,iElem)
    
    value = computeSource(theProblem,iElem) 
    map = [3*iElem+j for j in range(3)]
    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
        A[map[i],map[j]] += -(betax*dphidx[i] + betay*dphidy[i]) * jac / 6.0
      B[map[i]] += value[i]
      
#
# -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) 

    beta = nx*betax+ny*betay
    if (beta > 0) :
      alpha = 1.0
    else :
      alpha = 0.0
    
    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] ,mapEdgeLeft[j] ] += beta*alpha * 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
  
       
#
# -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) 
 
    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))
    
    beta = nx*betax+ny*betay
    if (beta > 0) :
      alpha = 1.0
    else :
      alpha = 0.0
 
    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 
        A[mapEdgeLeft[i] ,mapEdgeLeft[j] ] += beta*alpha * phiphi[i,j] * jac / 4.0
        A[mapEdgeLeft[i] ,mapEdgeRight[j]] -= beta*(1-alpha) * phiphi[i,j] * jac / 4.0 
        A[mapEdgeRight[i],mapEdgeLeft[j] ] -= beta*alpha * phiphi[i,j] * jac / 4.0 
        A[mapEdgeRight[i],mapEdgeRight[j]] += beta*(1-alpha) * 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