boundarySoluce.py

# -------------------------------------------------------------------------
#
# PYTHON for FEM DUMMIES 18-19
# Problème 2
#
# Solution détaillée du calcul des frontières
#
#  Vincent Legat
#
# -------------------------------------------------------------------------
# 

import numpy as np
import glfem 
import fem

#
# -1- Remplir la stucture de données à partir du maillage
#     On pourrait sans doute vectoriser tout cela....
#     Mais ce n'est pas l'essentiel du devoir !
#


def expand(theEdges):
  theMesh = theEdges.mesh
  edges = theEdges.edges
  for i in range (theMesh.nElem) :
    for j in range(3) :
      id = i*3 + j
      edges[id][0] = theMesh.elem[i][j]
      edges[id][1] = theMesh.elem[i][(j+1)%3]
      edges[id][2] = i
      edges[id][3] = -1

#
# -2- Le fonction de tri...
#     Ajouter le signe - final permet de commencer par le plus grand des plus petits :-)
#

      
def sort(theEdges):
   theEdges.edges.sort(key = lambda item : -(min(item[0:2])*theEdges.nEdges)-max(item[0:2]))

#
# -3- Détection des segments frontières et réduction en place du tableau
#     A nouveau, la vectorisation dans un nouveau tableau peut être plus efficace....
#

  
def shrink(theEdges):
  edges = theEdges.edges
  index = 0
  for i in range(theEdges.nEdges) :
     if (edges[i][0:2] != edges[i-1][1::-1]) :
         edges[index] = edges[i]
         index += 1
     else :
         edges[index-1][3] = edges[i][2]
  del edges[index:]
  theEdges.nBoundary = 2*index - theEdges.nEdges
  theEdges.nEdges = index   

#
# -4- Calcul de la frontière.
#     Chaque segment est considéré comme une droite reliant deux points de la sphère !
#     comme spécifié dans l'énoncé : si, si !
#

def lengthBoundary(theEdges):
  L = 0
  R = 6371000
  for myEdge in theEdges.edges:
    if (myEdge[3] == -1) :
      nodes = myEdge[:2]
      lon = theEdges.mesh.X[nodes] * (np.pi/180)
      lat = theEdges.mesh.Y[nodes] * (np.pi/180)
      x = R * np.cos(lat) * np.cos(lon)
      y = R * np.cos(lat) * np.sin(lon)
      z = R * np.sin(lat)
      dx = x[1]-x[0]; dy = y[1]-y[0]; dz = z[1]-z[0]; 
      L += np.sqrt(dx*dx + dy*dy + dz*dz)
  return L