hammerSoluce.py

# -------------------------------------------------------------------------
#
# PYTHON for FEM DUMMIES 18-19
# Problème 1
#
# Solution détaillée
#  Vincent Legat
#
# -------------------------------------------------------------------------
# 

from numpy import *
import glfem 

#
# -1- Implémentation de la règle de Hammer
#     Ne pas oublier la valeur absolue du jacobien :-)
#

def interpolate(u,xsi,eta) :
  return u[0]*xsi + u[1]*eta + u[2]*(1.0-xsi-eta)

def hammerIntegrate(x,y,f,graphics) :
  I = 0
  xsi    = [0.166666666666667,0.666666666666667,0.166666666666667]
  eta    = [0.166666666666667,0.166666666666667,0.666666666666667]
  weight = [0.166666666666667,0.166666666666667,0.166666666666667]  
  jac = abs((x[0]-x[1]) * (y[0]-y[2]) - (x[0]-x[2]) * (y[0]-y[1]))  
  xLoc = zeros(3)
  yLoc = zeros(3)
  for i in range(3) :
    xsiLoc = xsi[i]
    etaLoc = eta[i]
    xLoc[i] = interpolate(x,xsiLoc,etaLoc);
    yLoc[i] = interpolate(y,xsiLoc,etaLoc);     
    I += f(xLoc[i],yLoc[i]) * weight[i];
  if (graphics) :
    glfem.drawElement(x,y,'-k')
    glfem.drawNodes(x,y,'or')
    glfem.drawNodes(xLoc,yLoc,'ob')
  return I*jac

#
# -2- Implémentation récursive en divisant chaque triangle en 4
#     Sans doute pas le plus efficace, mais un bon compromis
#     entre efficacité et simplicité :-)
#

def hammerIntegrateRecursive(x,y,f,n,graphics) :
  nodes = [[0,3,5],[3,1,4],[5,4,2],[3,4,5]]
  xsi   = [0.0,1.0,0.0,0.5,0.5,0.0]
  eta   = [0.0,0.0,1.0,0.0,0.5,0.5]
  if (n <= 0) :
    return hammerIntegrate(x,y,f,graphics)
  I = 0
  xLoc = zeros(3)
  yLoc = zeros(3)
  for i in range(4) :
    for j in range(3) :
      xsiLoc = xsi[nodes[i][j]]
      etaLoc = eta[nodes[i][j]]
      xLoc[j] = interpolate(x,xsiLoc,etaLoc)
      yLoc[j] = interpolate(y,xsiLoc,etaLoc)
    I += hammerIntegrateRecursive(xLoc,yLoc,f,n-1,graphics)
  return I;