earthSoluce.py

# -------------------------------------------------------------------------
#
# PYTHON for FEM DUMMIES 18-19
# Problème 2
#
# Solution détaillée mais pas totalement optimisée : 
# en d'autres mots, c'est juste mais pas très repide...
# Seule, l'interpolation a été vectorisée.
#
# Pour le moment, je souhaite juste un code correct !
# Pour l'optimisation de l'intégration, 
#    voir le notebook de Louis, Jérome et Gilles,
#    mais ne pas trop regarder ce qu'ils racontent sur le jacobien toutefois !
#
#  Vincent Legat
#
# -------------------------------------------------------------------------
# 

import numpy as np
import glfem 
import fem

#
# -1- Le rayon de la Terre d'après Wikipedia :-)
#

def radius() :
  R = 6371008
  return R
 
#
# -2- Interpolation de la bathymétrie
#
#     La solution de base non optimisée :-)
#     Attention : il faut tester la sortie éventuelle de la boite avant
#     d'avoir arrondi les réels car int(-0.1) = 0 
#     alors que le point est dehors de la grille : merci Célestin qui a trouvé
#     le bug avec 8 années de retard :-)
#     
#

def gridInterpolateBasic(theGrid,x,y) :
  rx = (x - theGrid.Ox)/theGrid.dx
  ry = (y - theGrid.Oy)/theGrid.dx
    
  if (rx < 0 or rx >= theGrid.nx or ry < 0 or ry >= theGrid.ny) : 
    return 0
    
  i = int(rx)
  j = int(ry)
  alpha = rx - i
  beta =  ry - j
  elem = theGrid.elem

  bath = ((1-alpha)*(1-beta) * elem[j][i] 
          + alpha*(1-beta) * elem[j][i+1] 
          + alpha*beta * elem[j+1][i+1] 
          + (1-alpha)*beta * elem[j+1][i])
  return max(0,-bath)

#
# -2- Interpolation de la bathymétrie
#
#     Solution vectorisée pour pouvoir vectoriser efficacement...
#     l'intégration du volume :-)
#      
#     Finalement, vachement proche de la version basique....
#     mais, il faut juste gérer habilement le test !
#     Merci à Louis, Gilles et Jérome au passage : je me suis un peu inspiré 
#     de leur solution...   
#

	
def gridInterpolate(theGrid,x,y):
  rx = (x - theGrid.Ox)/theGrid.dx
  ry = (y - theGrid.Oy)/theGrid.dx
  
  outside = (rx >= 0) * (rx < theGrid.nx) * (ry >= 0) * (ry < theGrid.ny)
  
  i = rx.astype(int)*outside
  j = ry.astype(int)*outside
  alpha = rx-i
  beta  = ry-j
 
  bath = ((1-alpha)*(1-beta) * theGrid.elem[j,i] 
           + alpha*(1-beta) * theGrid.elem[j,i+1] 
           + alpha*beta * theGrid.elem[j+1,i+1] 
           + (1-alpha)*beta * theGrid.elem[j+1,i])
  return np.maximum(0,-bath)*outside

#
# -3- Intégration numérique avec un jacobien courbe sur un élément courbe de la sphère
#     Implémentation de base non optimisée et non vectorisée :-)
#
#     On vous laisse réfléchir à comment améliorer cela :-)
#     Louis, Gilles et Jérome ont de bonnes idées mais un jacobien mal foutu :-)
#
              
def integrateBathymetry(theMesh,theGrid,theRule,R) :
  I = 0;
  nGaussNumber    = theRule.n
  xsi             = theRule.xsi; 
  eta             = theRule.eta
  weight          = theRule.weight 

  for i in range(theMesh.nElem) :
    nodes = theMesh.elem[i]
    lon = theMesh.X[nodes]
    lat = theMesh.Y[nodes]
    for k in range(nGaussNumber) :
      xsiLoc = xsi[k]
      etaLoc = eta[k]
      weightLoc = weight[k]
      phi      = [xsiLoc,etaLoc,1.0-xsiLoc-etaLoc]
      dphidxsi = [ 1.0, 0.0,-1.0]
      dphideta = [ 0.0, 1.0,-1.0]
      lonLoc   = lon @ phi
      dlondxsi = lon @ dphidxsi
      dlondeta = lon @ dphideta
      latLoc   = lat @ phi
      dlatdxsi = lat @ dphidxsi
      dlatdeta = lat @ dphideta
      batLoc   = gridInterpolate(theGrid,lonLoc,latLoc)
      jacobian = abs(R*R*np.cos(latLoc*np.pi/180)*(dlondxsi*dlatdeta - dlondeta*dlatdxsi)*(np.pi/180)*(np.pi/180))
      I += batLoc*weightLoc*jacobian 
  return I

#
# -3- Intégration numérique avec un jacobien constant sur un élément plan 
#     C'est un peu plus rapide et la précision est en général tout-à-fait acceptable :-)
#
   
def integrateBathymetryBis(theMesh,theGrid,theRule,R) :
  I = 0;
  nGaussNumber    = theRule.n
  xsi             = theRule.xsi; 
  eta             = theRule.eta
  weight          = theRule.weight 

  for i in range(theMesh.nElem) :
    nodes = theMesh.elem[i]
    lon = theMesh.X[nodes]
    lat = theMesh.Y[nodes]
    x     = R * np.cos(lon*np.pi/180) * np.cos(lat*np.pi/180)
    y     = R * np.sin(lon*np.pi/180) * np.cos(lat*np.pi/180)
    z     = R * np.sin(lat*np.pi/180)
    dxdxsi = x[1] - x[0]
    dydxsi = y[1] - y[0]
    dzdxsi = z[1] - z[0]
    dxdeta = x[2] - x[0]
    dydeta = y[2] - y[0]
    dzdeta = z[2] - z[0]
    jacx = dydxsi * dzdeta - dzdxsi * dydeta
    jacy = dzdxsi * dxdeta - dxdxsi * dzdeta
    jacz = dxdxsi * dydeta - dydxsi * dxdeta
    jacobian = np.sqrt(jacx * jacx + jacy * jacy + jacz * jacz);   

    for k in range(nGaussNumber) :
      xsiLoc = xsi[k]
      etaLoc = eta[k]
      weightLoc = weight[k]
      phi      = [xsiLoc,etaLoc,1.0-xsiLoc-etaLoc]
      lonLoc   = lon @ phi
      latLoc   = lat @ phi
      batLoc   = gridInterpolate(theGrid,lonLoc,latLoc)
      I += batLoc*weightLoc*jacobian 
  return I