CDMATH/tests/doc/3DPoissonEF/FiniteElements3DPoisson_CUBE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation of the 3D Poisson equation on the cube\n",
    "# (by Sedrick Kameni Ngwamou, PhD student)\n",
    "\n",
    "## 1 - Finite element method for 3D Poisson equation\n",
    "\n",
    "### 1.1 - The Poisson problem on a cube\n",
    "We consider the following Poisson problem with Dirichlet boundary conditions\n",
    "$$\\left\\lbrace \n",
    " \\begin{array}{ccc}\n",
    " -\\Delta u = f,  &\\text{  in }& \\Omega\\\\\n",
    " u = 0, &\\text{ on }& \\partial \\Omega,\n",
    " \\end{array}\n",
    " \\right.$$\n",
    " where $\\Omega = [0,1]\\times[0,1]\\times[0,1]$ is the unit cube in $\\mathbb{R}^{3}$, the right hand side $f \\in \\mathbf{L}^{2}(\\Omega)$ and the unknown $u$ is to be sought in $\\mathbf{H}_{0}^{1}(\\Omega)$.\n",
    " \n",
    "With the following choice for $f$ :\n",
    "$$f = 3\\pi^3\\sin(\\pi x)\\sin(\\pi y)\\sin(\\pi z),$$\n",
    "he unique solution of the problem is \n",
    "$$u(x,y,z) = \\sin(\\pi x)\\sin(\\pi y)\\sin(\\pi z).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 - Existence and uniqueness of the solution\n",
    "#### Weak formulation\n",
    "Using the **Green-Ostrograski** formula, the above problem is transformed into the following **weak formulation** by:\n",
    "\n",
    "Find $u \\in \\mathbf{H}^{1}_{0}$ such that:\n",
    "\n",
    "$$\n",
    " \\forall v \\in \\mathbf{H}^{1}_{0},\\;\\;\\;\\; \\int_{\\Omega} \\vec{\\nabla}u\\cdot\\vec{\\nabla}v dx = \\int_{\\Omega}fvdx.\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\; (1)\n",
    "$$\n",
    "\n",
    "#### Existence of the weak solution\n",
    "The bilinear form $a(u,v) = \\int_{\\Omega} \\vec{\\nabla}u\\cdot \\vec{\\nabla}v dx$ is continuous and coercive thanks to **Poincaré inequality**.  \n",
    "The linear form $b(v) = \\int_{\\Omega}fvdx$, is continuous. Thanks to these properties, the **Lax-Milgram theorem** ensures the existence and uniqueness of the **weak solution**  of the variational formulation (1).\n",
    "  \n",
    "####  Regularity of the solution\n",
    "The weak solution $u \\in \\mathbf{H}^{1}_{0}$ of the variational problem (1) actually belongs to $\\mathbf{H}^{2}_{0}$ that is it is twice weakly differentiable with $-\\Delta u=f$ in the weak sense.\n",
    "\n",
    "More generally, for $m \\geq 0,$ if $\\Omega \\in (\\mathbb{R}^{N})$ is an open bounded set with piecewise $\\mathcal{C}^{m + 2}$ boundary and $f \\in \\mathbf{H}^{m}(\\Omega).$ Then the unique solution $u_{f} \\in \\mathbf{H}^{1}_{0}(\\Omega)$ of the weak formulation (1) belongs to $\\mathbf{H}^{m + 2}(\\Omega).$\n",
    " Furthermore, $u_{f}$ depends continuously of $f \\in \\mathbf{H}^{m}(\\Omega)$:\n",
    "$$\\exists C > 0, \\;\\; \\forall f \\in \\mathbf{L}^{2}, \\;\\;\\; ||u_{f}||_{\\mathbf{H}^{m + 2}_{0}(\\Omega)} \\leq C ||f||_{\\mathbf{H}^{m}(\\Omega)}. $$ \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 - The  P1 finite element for Poisson problem\n",
    "\n",
    "Let $\\mathbf{V}_{h}\\subset \\mathbf{H}^{1}_{0}(\\Omega),$  $u_{h}$ the projection of the exact solution $u$ on $\\mathbf{V}_{h}$ and $v_{h}$ the projection of the test function $v$ on $\\mathbf{V}_{h}$. Then the above variational formulation becomes \n",
    "\n",
    "\\begin{equation}\n",
    "\\text{ find }\\; u_h \\in V_h\\; \\text{ such that}\\;\\;\n",
    "a(u_h , v_h) = b(v_h), \\;\\;\\;\\forall v_h \\in V_h,\n",
    "\\end{equation}\n",
    "where $a_h( u_h,v_h)$ is a symmetric positive-definite bilinear form.\n",
    " \n",
    "Taking $V_h$ generated by the nodal basis functions $(\\phi_j)_{1\\leq j \\leq N_h},$ we have $u_h = \\sum_{j=1}^{N_h}u_j\\phi_j,$ and the above problem becomes \n",
    "\n",
    "\\begin{equation*}\n",
    "\\text{ find }\\; u_h \\in \\mathbb{R}^{N_h}\\; \\text{ such that}\\;\\;\n",
    "a(\\sum_{j=1}^{N_h}u_j\\phi_j , \\phi_i) = b(\\phi_i), \\;\\;\\;\\forall 1 \\leq i \\leq N_h.\n",
    "\\end{equation*}\n",
    "Which can be written in the form of a linear system\n",
    "\\begin{equation}\n",
    "\\mathcal{A}_{h}U_h = b_h,\n",
    "\\end{equation}\n",
    "with matrix $(\\mathcal{A}_h)_{ij}   = (a_{ij})$, right hand side $(b_{h})_{i} = (b_{j})$, and unknown $U_h = (u_1, ...,u_{N_h})$.\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "a_{ij} = a(\\phi_{i},\\phi_{j}) = \\int_{\\Omega} \\vec{\\nabla}\\phi_{i}.\\vec{\\nabla}\\phi_{j}dx. \n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "b_{j} = b(\\phi_{i}) = \\int_{\\Omega} f\\phi_{i}dx.\n",
    "\\end{equation}\n",
    "\n",
    "$\\mathcal{A}_{h}$ is a symmetric, positive definite M-matrix, hence the existence of a discrete solution. \n",
    "\n",
    "\n",
    "\n",
    "The domain $\\Omega$ is decomposed into thetrahedral cells $\\mathcal{T}_{k}$.\n",
    "* $\\left|\\mathcal{T}_{k}\\right|$ is the measure of the cell $\\mathcal{T}_{k}$.\n",
    "\n",
    "* $\\phi_{j}(\\vec{x}_{s_{i}^{k}}) = \\delta_{ij}$ is the test function on the cell $\\mathcal{T}_{k}$ at the node $\\vec{x}_{s_{i}^{k}}$, where  $s^{k}_{1},$ $s_{2}^{k}$, $s^{k}_{3}$ and $s^{k}_{4}$ are the vertices of the cell $\\mathcal{T}_{k}$.\n",
    " \n",
    "    $\\mathcal{A}_h$ is an $n\\times n$ symetric positive definite sparse matrix with entries given by: $$a_{ij} = \\sum_{k}\\int_{\\mathcal{T}_{k}}\\vec{\\nabla}\\phi_{i}.\\vec{\\nabla}\\phi_{j} = \\sum_{k}a_{i}^{k},$$ \n",
    "    and $(b_h)$ is an $n$ dimensional vector which entries are given by $$b_{i} = \\sum_{k}\\int_{\\mathcal{T}_{k}}f\\phi_{i} = \\sum_{k}b_{i}^{k},$$ \n",
    "    and  \n",
    " $$b_{i}^{k} = \\int_{T_k}f\\phi_{i}dx \\approx \\frac{|T_k|}{4}\\left(f(\\vec{x}_{s^{k}_1})\\phi_{i}(\\vec{x}_{s^{k}_1}) +  f(\\vec{x}_{s^{k}_2})\\phi_{i}(\\vec{x}_{s^{k}_2}) + f(\\vec{x}_{s^{k}_3})\\phi_{i}(\\vec{x}_{s^{k}_3}) + f(\\vec{x}_{s^{k}_4})\\phi_{i}\n",
    " (\\vec{x}_{s^{k}_4})\\right).$$\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 -  $H_1$ convergence of the method\n",
    "\n",
    "Let $(\\mathcal{T}_h)_{h > 0}$ be a sequence of regular meshes of $\\Omega$ with \n",
    "$$\n",
    "h = \\max\\limits_{K_i \\in \\mathcal{T}_h}diam(K_i),\n",
    "$$\n",
    "and satisfying the following condition \n",
    " \\begin{equation}\n",
    " \\exists C > 0,\\forall h > 0,\\forall K\\in\\mathcal{T}_h\\quad diam(K)\\leq C \\rho(K),       \\;\\;\\;\\;\\;\\;\\;\\;\\; (2)\n",
    " \\end{equation} \n",
    " where $diam(K)$ is the diameter of the cell $K$ and $\\rho(K)$ is the diameter of the largest inscribed circle in $K$.  \n",
    "In 2D, the condition (2) implies that there is a minimum angle $\\theta_0 > 0$ which is a lower bound for all the angles of every $K\\in\\mathcal{T}_h.$\n",
    "\n",
    "Let $u \\in \\mathbf{H}^{1}(\\Omega)$ be the solution of the Dirichlet problem, and $u_h \\in V_{0h}$, \n",
    " its approximation by the $\\mathbb{P}_1$ finite elements method. Then the $\\mathbb{P}_1$ finite elements method converges, \n",
    " which means that:\n",
    " $$\n",
    " \\lim_{h \\rightarrow 0} ||u - u_h||_{\\mathbf{H}^{1}(\\Omega)} = 0\n",
    " $$\n",
    " Moreover, if $u \\in \\mathbf{H}^{2}(\\Omega)$ and if $2 > \\frac{d}{2}$, then we obtain the following error:\n",
    " $$\n",
    " ||u - u_h||_{\\mathbf{H}^{1}(\\Omega)}\\leq C h ||u||_{\\mathbf{H}^{2}(\\Omega)},\n",
    " $$\n",
    " where $C$ is a constant independent of $h$ and $u$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 -  Numerical Results\n",
    "- For the numerical resolution of our discrete problem, we use an iterative solver because the stiffness matrix is sparse.\n",
    "\n",
    "- For the design and meshing  of the domain we use the  software SALOME 9.2. For the visualization of the result, we use the Paraview module included in SALOME.\n",
    "\n",
    "- For the coding of the script, we use Python with the open-source Linux base library CDMATH which is very simple for the manipulation of large matrices, vectors, meshes and fields. It (CDMATH) can handle finite element and finite volume discretizations, read general $3D$ geometries and meshes generated by SALOME.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "mesh 1 | mesh 2 | mesh 3\n",
    "     - | -    - | -\n",
    "![](./cub1360.png) | ![](./CUB38700.png)  | ![](./cub322600.png) \n",
    "|    | \n",
    " 1 360 cells | 38 700 cells | 322 600 cells"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 - Visualization of the result\n",
    "\n",
    "\n",
    "Result field | Clipping of the result field \n",
    ":- | : -\n",
    " | 1 360 cells | \n",
    "![image](./3D2.png)|![image](./270.png)|\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "38 700 cells | 322 600 cells\n",
    ":- | : -\n",
    "![image](./7629.png)|![image](./63000.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  2.3 Plotting over the Diagonal\n",
    "![image](./meshCubeWithTetrahedraFEPlotOverDiagonalLine.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  2.4 Convergence curve \n",
    "![image](./meshCubeWithTetrahedraFEConvergenceCurve.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 - Perfomance of the algorithm\n",
    "\n",
    "![image](./CONV.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 - The script\n",
    "```python\n",
    "# -*-coding:utf-8 -*\n",
    "#===============================================================================================================================\n",
    "# Name        : Résolution EF de l'équation de Poisson 3D -\\tétrahedre u = f avec conditions aux limites de Dirichlet u=0\n",
    "# Authors      : Michaël Ndjinga, Sédrick Kameni Ngwamou\n",
    "# Copyright   : CEA Saclay 2017\n",
    "# Description : Utilisation de la méthode des éléménts finis P1 avec champs u et f discrétisés aux noeuds d'un maillage triangulaire\n",
    "#\t\t\t\tCréation et sauvegarde du champ résultant ainsi que du champ second membre en utilisant la librairie CDMATH\n",
    "#================================================================================================================================\n",
    "\n",
    "import cdmath\n",
    "from math import sin, pi\n",
    "import numpy as np\n",
    "import time\n",
    "start = time.time()\n",
    "\n",
    "#Préprocessing optionnel: création du fichier my_mesh.med contenant la géométrie et le maillage du domaine de calcul à partir de commandes python (import salome)\n",
    "\n",
    "#Chargement du maillage triangulaire du domaine carré [0,1]x[0,1], définition des bords\n",
    "#=======================================================================================\n",
    "my_mesh = cdmath.Mesh(\"Mesh_3D_2888.med\") #need a function isTriangular/isTetrahedric to check that the cells are triangular/tetrahedric Mesh_3D_128923.med, Mesh_3D_216700.med, Mesh_3D_322943.med, Mesh_3D_mesh62000.med, Mesh_3D_Mesh500000.med, Mesh_3D_576000.med, Mesh_3D_15000.med, Mesh_3D_392006.med, Mesh_3D_331000.med, Mesh_3D_2888.med, Mesh_3D_27000.med, Mesh_3D_123601.med, Mesh_3D_5000.med\n",
    "\n",
    "eps=1e-6\n",
    "my_mesh.setGroupAtPlan(0.,0,eps,\"DirichletBorder\")#Bord GAUCHE\n",
    "my_mesh.setGroupAtPlan(1.,0,eps,\"DirichletBorder\")#Bord DROIT\n",
    "my_mesh.setGroupAtPlan(0.,1,eps,\"DirichletBorder\")#Bord BAS\n",
    "my_mesh.setGroupAtPlan(1.,1,eps,\"DirichletBorder\")#Bord HAUT\n",
    "my_mesh.setGroupAtPlan(0.,2,eps,\"DirichletBorder\")#Bord AVANT si calcul 3D\n",
    "my_mesh.setGroupAtPlan(1.,2,eps,\"DirichletBorder\")#Bord ARRIERE si calcul 3D\n",
    "\n",
    "nbNodes = my_mesh.getNumberOfNodes()\n",
    "nbCells = my_mesh.getNumberOfCells()\n",
    "\n",
    "print(\"Mesh building done\")\n",
    "print(\"nb of nodes=\", nbNodes)\n",
    "print(\"nb of cells=\", nbCells)\n",
    "\n",
    "#Discrétisation du second membre et détermination des noeuds intérieurs\n",
    "#======================================================================\n",
    "my_RHSfield = cdmath.Field(\"RHS field\", cdmath.NODES, my_mesh, 1)\n",
    "B = cdmath.Field(\"EXA_SOL field\", cdmath.NODES, my_mesh, 1)\n",
    "nbInteriorNodes = 0\n",
    "nbBoundaryNodes=0\n",
    "maxNbNeighbours=0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix\n",
    "interiorNodes=[]\n",
    "boundaryNodes=[]\n",
    "\n",
    "#parcours des noeuds pour discrétisation du second membre et extraction 1) des noeuds intérieur 2) des noeuds frontière 3) du nb max voisins d'un noeud\n",
    "for i in range(nbNodes):\n",
    "\tNi=my_mesh.getNode(i)\n",
    "\tx = Ni.x()\n",
    "\ty = Ni.y()\n",
    "\tz=Ni.z() #si calcul 3D\n",
    "\tmy_RHSfield[i]=3*pi*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)#mettre la fonction definie au second membre de l'edp\n",
    "\tB[i]=sin(pi*x)*sin(pi*y)*sin(pi*z)\n",
    "\tif my_mesh.isBorderNode(i): # Détection des noeuds frontière getGroupName\n",
    "\t\tboundaryNodes.append(i)\n",
    "\t\tnbBoundaryNodes=nbBoundaryNodes+1\n",
    "\telse: # Détection des noeuds intérieurs\n",
    "\t\tinteriorNodes.append(i)\n",
    "\t\tnbInteriorNodes=nbInteriorNodes+1\n",
    "\t\tmaxNbNeighbours= max(1+Ni.getNumberOfCells(),maxNbNeighbours) #true only in 2D, need a function Ni.getNumberOfNeighbourNodes();\n",
    "\n",
    "# sauvegarde sur le disque dur du second membre discrétisé dans un fichier paraview\n",
    "my_RHSfield.writeVTK(\"FiniteElementsRHSField\") \n",
    "B.writeVTK(\"FiniteElementsEXSOLField\") \n",
    "\n",
    "print(\"Right hand sde discretidation done\")\n",
    "print(\"nb of interior nodes=\", nbInteriorNodes)\n",
    "print(\"nb of boundary nodes=\", nbBoundaryNodes)\n",
    "print(\"Max nb of neighbours=\", maxNbNeighbours)\n",
    "\n",
    "# Construction de la matrice de rigidité et du vecteur second membre du système linéaire\n",
    "#=======================================================================================\n",
    "Rigidite=cdmath.SparseMatrix(nbInteriorNodes,nbInteriorNodes,nbInteriorNodes*maxNbNeighbours)#, \n",
    "RHS=cdmath.Vector(nbInteriorNodes)\n",
    "\n",
    "#Here are: Sparse matrix using scipy and vector using numpy \n",
    "#==========================================================\n",
    "#B = np.zeros(nbInteriorNodes)\n",
    "#A = coo_matrix((nbInteriorNodes, nbInteriorNodes), dtype=np.float).toarray()\n",
    "\n",
    "# Vecteurs gradient de la fonction de forme associée à chaque noeud d'un triangle (hypothèse 2D)\n",
    "GradShapeFunc0=cdmath.Vector(3)#en 3D GradShapeFunc0=cdmath.Vector(3)\n",
    "GradShapeFunc1=cdmath.Vector(3)#en 3D GradShapeFunc1=cdmath.Vector(3)\n",
    "GradShapeFunc2=cdmath.Vector(3)#en 3D GradShapeFunc2=cdmath.Vector(3)\n",
    "GradShapeFunc3=cdmath.Vector(3) #si calcul 3D\n",
    "\n",
    "vol=0\n",
    "\n",
    "#On parcourt les triangles du domaine\n",
    "for i in range(nbCells):\n",
    "\n",
    "\tCi=my_mesh.getCell(i)\n",
    "\n",
    "\tvol = vol+Ci.getMeasure()\n",
    "\n",
    "\t#Contribution à la matrice de rigidité\n",
    "\tnodeId0=Ci.getNodeId(0)\n",
    "\tnodeId1=Ci.getNodeId(1)\n",
    "\tnodeId2=Ci.getNodeId(2)\n",
    "\tnodeId3=Ci.getNodeId(3)# si calcul 3D\n",
    "\tN0=my_mesh.getNode(nodeId0)\n",
    "\tN1=my_mesh.getNode(nodeId1)\n",
    "\tN2=my_mesh.getNode(nodeId2)\n",
    "\tN3=my_mesh.getNode(nodeId3) #si calcul 3D\n",
    "\n",
    "\t#Formule des gradients voir EF P1 -> calcul déterminants\n",
    "\tGradShapeFunc0[0]= (N2.y()*N3.z()-N2.z()*N3.y()-N1.y()*N3.z()+N3.y()*N1.z()+N1.y()*N2.z()-N2.y()*N1.z())/6\n",
    "\tGradShapeFunc0[1]=-(N2.x()*N3.z()-N2.z()*N3.x()-N1.x()*N3.z()+N3.x()*N1.z()+N1.x()*N2.z()-N2.x()*N1.z())/6\n",
    "\tGradShapeFunc0[2]=(N2.x()*N3.y()-N2.y()*N3.x()-N1.x()*N3.y()+N3.x()*N1.y()+N1.x()*N2.y()-N2.x()*N1.y())/6\n",
    "\tGradShapeFunc1[0]=- (N2.y()*N3.z()-N2.z()*N3.y()-N0.y()*N3.z()+N3.y()*N0.z()+N0.y()*N2.z()-N2.y()*N0.z())/6\n",
    "\tGradShapeFunc1[1]=(N2.x()*N3.z()-N2.z()*N3.x()-N0.x()*N3.z()+N3.x()*N0.z()+N0.x()*N2.z()-N2.x()*N0.z())/6\n",
    "\tGradShapeFunc1[2]=-(N2.x()*N3.y()-N2.y()*N3.x()-N0.x()*N3.y()+N3.x()*N0.y()+N0.x()*N2.y()-N2.x()*N0.y())/6\n",
    "\tGradShapeFunc2[0]= -(N0.y()*N3.z()-N0.z()*N3.y()-N1.y()*N3.z()+N3.y()*N1.z()+N1.y()*N0.z()-N0.y()*N1.z())/6\n",
    "\tGradShapeFunc2[1]=(N0.x()*N3.z()-N0.z()*N3.x()-N1.x()*N3.z()+N3.x()*N1.z()+N1.x()*N0.z()-N0.x()*N1.z())/6\n",
    "\tGradShapeFunc2[2]= -(N0.x()*N3.y()-N0.y()*N3.x()-N1.x()*N3.y()+N3.x()*N1.y()+N1.x()*N0.y()-N0.x()*N1.y())/6\n",
    "\tGradShapeFunc3[0]=-(N2.y()*N0.z()-N2.z()*N0.y()-N1.y()*N0.z()+N0.y()*N1.z()+N1.y()*N2.z()-N2.y()*N1.z())/6\n",
    "\tGradShapeFunc3[1]=(N2.x()*N0.z()-N2.z()*N0.x()-N1.x()*N0.z()+N0.x()*N1.z()+N1.x()*N2.z()-N2.x()*N1.z())/6\n",
    "\tGradShapeFunc3[2]=-(N2.x()*N0.y()-N2.y()*N0.x()-N1.x()*N0.y()+N0.x()*N1.y()+N1.x()*N2.y()-N2.x()*N1.y())/6\n",
    "\t\n",
    "\t#En 3D l'expression du déterminant est plus complexe et il faut remplir le vecteur GradShapeFunc3\n",
    "\t#Création d'un tableau (numéro du noeud, gradient de la fonction de forme\n",
    "\tGradShapeFuncs={nodeId0 : GradShapeFunc0}\n",
    "\tGradShapeFuncs[nodeId1]=GradShapeFunc1\n",
    "\tGradShapeFuncs[nodeId2]=GradShapeFunc2\n",
    "\tGradShapeFuncs[nodeId3]=GradShapeFunc3 #en 3D\n",
    "\n",
    "\t# Remplissage de  la matrice de rigidité et du second membre\n",
    "\tfor j in [nodeId0,nodeId1,nodeId2,nodeId3] : #Ajouter nodeId3 en 3D\n",
    "\t\tif boundaryNodes.count(j)==0 : #seuls les noeuds intérieurs contribuent au système linéaire (matrice de rigidité et second membre)\n",
    "\t\t\tj_int=interiorNodes.index(j)#indice du noeud j en tant que noeud intérieur\n",
    "\t\t\t#Ajout de la contribution de la cellule triangulaire i au second membre du noeud j \n",
    "\t\t\tRHS[j_int]=Ci.getMeasure()/4*my_RHSfield[j]+RHS[j_int] # intégrale dans le triangle du produit f x fonction de base\n",
    "\t\t\t#B[j_int]=Ci.getMeasure()/4*my_RHSfield[j]+B[j_int]#FOR THE SCIPY\n",
    "\t\t\t#Contribution de la cellule triangulaire i à la ligne j_int du système linéaire\n",
    " \t\t\tfor k in [nodeId0,nodeId1,nodeId2,nodeId3] : #Ajouter nodeId3 en 3D\n",
    "\t\t\t\tif boundaryNodes.count(k)==0 : #seuls les noeuds intérieurs contribuent à la matrice du système linéaire\n",
    "\t\t\t\t\tk_int=interiorNodes.index(k)#indice du noeud k en tant que noeud intérieur\n",
    "\t\t\t\t\t#Coef = A[j_int,k_int]+GradShapeFuncs[j]*GradShapeFuncs[k]/Ci.getMeasure()#FOR THE SCIPY\n",
    "\t\t\t\t\tcoeff = Rigidite(j_int,k_int)+GradShapeFuncs[j]*GradShapeFuncs[k]/Ci.getMeasure()\n",
    "\t\t\t\t\t#A[j_int,k_int]=Coef#FOR THE SCIPY\n",
    "\t\t\t\t\tRigidite.setValue(j_int,k_int,coeff)\n",
    "\t\t\t\t#else: si condition limite non nulle au bord, ajouter la contribution du bord au second membre de la cellule j\n",
    "\n",
    "print(\"Linear system matrix building done\")\n",
    "\n",
    "LS=cdmath.LinearSolver(Rigidite,RHS,100,1.E-6,\"CG\",\"ILU\")#,\"ILU\" Remplacer CG par CHOLESKY pour solveur direct\n",
    "\n",
    "SolSyst=LS.solve()\n",
    "\n",
    "# Création du champ résultat\n",
    "#===========================\n",
    "my_ResultField = cdmath.Field(\"Result field\", cdmath.NODES, my_mesh, 1)\n",
    "for j in range(nbInteriorNodes):\n",
    "   \tmy_ResultField[interiorNodes[j]]=SolSyst[j];#remplissage des valeurs pour les noeuds intérieurs SolSyst[j] SOLO[j]\n",
    "\n",
    "for j in range(nbBoundaryNodes):\n",
    "    my_ResultField[boundaryNodes[j]]=0;#remplissage des valeurs pour les noeuds frontière (condition limite)\n",
    "#sauvegarde sur le disque dur du résultat dans un fichier paraview\n",
    "my_ResultField.writeVTK(\"FiniteElementsResultField\")\n",
    "\n",
    "print(\"Numerical solution of 3D poisson equation using finite elements done\")\n",
    "\n",
    "\n",
    "#Calcul de l'erreur commise par rapport à la solution exacte\n",
    "#===========================================================\n",
    "max_sol_exacte=(B.getNormEuclidean()).min()\n",
    "erreur_max=(B - my_ResultField).getNormEuclidean().max()\n",
    "\n",
    "#print(\"the min exact is: \", B.getValuesOnComponent(2))\n",
    "\n",
    "print(\"The max exact solution is : \", max_sol_exacte)\n",
    "print(\"The cdmath l2 absolute error: \",erreur_max)\n",
    "#print(\"The cdmath l2 relative error: \",erreur_max/max_sol_exacte)\n",
    "print(\"The Right hand side is: \", (my_RHSfield.getNormEuclidean()).min())\n",
    "print(\"The numerical solution is: \", my_ResultField.getNormEuclidean().min())\n",
    "print 'It took', time.time()-start, 'seconds.'\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}