# ------------------------------------------------------------------------- # # 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()