poissonSoluce.py

# -------------------------------------------------------------------------
#
# PYTHON for FEM DUMMIES 18-19
# Problème 4
#
# 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- Assemblage de la matrice 
#     Il est possible de calculer analytiquement tous les termes !
#     C'est donc totalement inutile de faire appel à une règle d'intégration :-)
#

def assemble(theProblem):
  theMesh         = theProblem.mesh
  theSystem       = theProblem.system
  theEdges        = theProblem.edges
  A = theSystem.A
  B = theSystem.B

  for iElem in range(theMesh.nElem) :
    nodes = theMesh.elem[iElem]
    x = theMesh.X[nodes]
    y = theMesh.Y[nodes]
    dphidxsi = np.array([ 1.0, 0.0,-1.0])
    dphideta = np.array([ 0.0, 1.0,-1.0])
    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

    for i in range(3):
      for j in range(3):
        A[nodes[i],nodes[j]] += (dphidx[i] * dphidx[j] + dphidy[i] * dphidy[j]) *jac / 2.0
      B[nodes[i]] += jac / 6.0
    
  for myEdge in theEdges.edges:
    if (myEdge[3] == -1) :
      theSystem.constrain(myEdge[0],0.0)
      theSystem.constrain(myEdge[1],0.0)

#
# -2- Optimisation de l'élimination gaussienne...
#     A faire :-)
#     D'ici peu, je ferai la synthèse des meilleures contributions.
#     En ne faisant pas cette partie, j'obtiens 8/10 :-)
#     Be patient :-)
#

      
def eliminate(theProblem):
  theProblem.system.eliminate()
       

  

• Download
• Back to MECA1120