#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
***********************************************************************************
                            tutorial16.py
                DAE Tools: pyDAE module, www.daetools.com
                Copyright (C) Dragan Nikolic
***********************************************************************************
DAE Tools is free software; you can redistribute it and/or modify it under the
terms of the GNU General Public License version 3 as published by the Free Software
Foundation. DAE Tools is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with the
DAE Tools software; if not, see <http://www.gnu.org/licenses/>.
************************************************************************************
"""
__doc__ = """
This tutorial shows how to use DAE Tools objects with NumPy arrays to solve a simple
stationary heat conduction in one dimension using the Finite Elements method
with linear elements and two ways of manually assembling a stiffness matrix/load vector:

.. code-block:: none

   d2T(x)/dx2 = F(x);  x in (0, Lx)

Linear finite elements discretisation and simple FE matrix assembly:

.. code-block:: none

                   phi                 phi
                      (k-1)               (k)
                      
                     *                   *             
                   * | *               * | *            
                 *   |   *           *   |   *          
               *     |     *       *     |     *        
             *       |       *   *       |       *      
           *         |         *         |         *        
         *           |       *   *       |           *      
       *             |     *       *     |             *    
     *               |   *           *   |               *  
   *                 | *  element (k)  * |                 *
 *-------------------*+++++++++++++++++++*-------------------*-
                     x                   x
                      (k-i                (k)
                    
                     \_________ _________/
                               |
                               dx

The comparison of the analytical solution and two ways of assembling the system is given
in the following plot:

.. image:: _static/tutorial16-results.png
   :width: 500px
"""

import sys, numpy
from time import localtime, strftime
from daetools.pyDAE import *

# Standard variable types are defined in variable_types.py
from pyUnits import m, kg, s, K, Pa, mol, J, W

class modTutorial(daeModel):
    def __init__(self, Name, Parent = None, Description = ""):
        daeModel.__init__(self, Name, Parent, Description)

        self.x  = daeDomain("x", self, unit(), "x axis domain")

        self.L  = daeParameter("L", unit(), self, "Length")

        self.Ta = daeVariable("Ta", no_t, self, "Temperature - analytical solution", [self.x])
        self.T1 = daeVariable("T1", no_t, self, "Temperature - first way",           [self.x])
        self.T2 = daeVariable("T2", no_t, self, "Temperature - second way",          [self.x])

    def local_dof(self, i):
        return self.mapLocalToGlobalIndices

    def DeclareEquations(self):
        daeModel.DeclareEquations(self)

        N = self.x.NumberOfPoints
        Nelem = N - 1
        Ndofs_per_elem = 2 # since we use linear elements
        
        numpy.set_printoptions(linewidth=1e10)
        ##################################################################################
        # Analytical solution
        ##################################################################################
        dx = self.L.GetValue() / Nelem
        m = [0., 5./2, 4., 9./2]
        for i in range(N):
            eq = self.CreateEquation("Poisson_Analytical(%d)" % i)
            eq.Residual = self.Ta(i) - m[i] * dx**2

        ##################################################################################
        # First way: use constant global stiffness matrix and load array
        ##################################################################################
        print('***************************************************************************')
        print('    First way')
        print('***************************************************************************')
        # Create global stiffness matrix and load vector (dtype = float):
        A = numpy.zeros((N,N))        
        F = numpy.zeros(N)

        # Maps local indices within an element to the indices in the system's matrix
        mapLocalToGlobalIndices = [(0, 1), (1, 2), (2, 3)]   
        # ∆x is equidistant and constant
        dx = self.L.GetValue() / Nelem
        for el in range(Nelem):
            # Get global indices for the current element
            dof_indices = mapLocalToGlobalIndices[el]
            
            # Element stiffness matrix and load vector:
            #       | 1 -1 |          | 1 |
            # Ael = |-1  1 |    Fel = | 1 |
            #       
            Ael = (1/dx) * numpy.array( [[1, -1], [-1, 1]] )
            Fel = (dx/2) * numpy.array( [1, 1] )
            
            # Loop over element DOFs and update the global matrix and vector 
            for i in range(Ndofs_per_elem):
                for j in range(Ndofs_per_elem):
                    A[dof_indices[i], dof_indices[j]] += Ael[i,j]
                F[dof_indices[i]] += Fel[i]
        
        print('The global stiffness matrix (A) before applying boundary conditions:')
        print(A)
        print('The global load vector (F):')
        print(F)
        
        # Boundary conditions:
        # at x = 0: T(0) = 0     (Dirichlet BC)
        # at x = 1: dT(1)/dx = 0 (Neumann BC)
        A[0, 1:-1] = 0
        F[0] = 0
        print('The global stiffness matrix (A) after applying boundary conditions:')
        print(A)
        print('The global load vector (F) after applying boundary conditions:')
        print(F)
        
        # Create a vector of temperatures:
        T = numpy.empty(N, dtype=object)
        T[:] = [self.T1(i) for i in range(N)]
        
        # Generate the system equations
        for i in range(N):
            eq = self.CreateEquation("Poisson_ConstantStiffnessMatrix(%d)" % i)
            eq.Residual = numpy.sum(A[i, :] * T[:]) - F[i]

        ##################################################################################
        # Second way: use global stiffness matrix and load array that depend on DAE Tools
        #             model parameters/variables (not constant in a general case).
        #             Obviously, they are constant here - this is only to show the concept
        ##################################################################################
        print('***************************************************************************')
        print('    Second way')
        print('***************************************************************************')
        # Create global stiffness matrix and load vector (dtype = object). This matrix and
        # load vector will be functions of model parameters/variables.  
        # In this simple example that is not a case; however, the procedure is analogous.
        A = numpy.zeros((N,N), dtype=object)
        #A[:] = adouble(0, 0, True)
        #print A
        
        F = numpy.zeros(N, dtype=object)
        #F[:] = adouble(0, 0, True)
        #print F
        
        # Maps local indices within an element to the indices in the system's matrix
        mapLocalToGlobalIndices = [(0, 1), (1, 2), (2, 3)]
        # ∆x is equidistant but not constant (it depends on the parameter 'L')
        dx = self.L() / Nelem
        for el in range(Nelem):
            # Get global indices for the current element
            dof_indices = mapLocalToGlobalIndices[el]
            
            # Element stiffness matrix and load vector are the same:
            #       | 1 -1 |          | 1 |
            # Ael = |-1  1 |    Fel = | 1 |
            #
            Ael = (1 / dx) * numpy.array( [[1, -1], [-1, 1]] )
            Fel = (dx / 2) * numpy.array([1, 1])
            
            # Loop over element DOFs and update the global matrix and vector 
            for i in range(Ndofs_per_elem):
                for j in range(Ndofs_per_elem):
                    A[dof_indices[i], dof_indices[j]] += Ael[i,j]
                F[dof_indices[i]] += Fel[i]
        
        print('The global stiffness matrix (A) before applying boundary conditions:')
        print(A)
        print('The global load vector (F):')
        print(F)
        
        # Boundary conditions:
        # at x = 0: T(0) = 0     (Dirichlet BC)
        # at x = 1: dT(1)/dx = 0 (Neumann BC)
        A[0, 1:-1] = 0
        F[0] = 0
        print('The global stiffness matrix (A) after applying boundary conditions:')
        print(A)
        print('The global load vector (F) after applying boundary conditions:')
        print(F)
        
        # Create a vector of temperatures:
        T = numpy.empty(N, dtype=object)
        T[:] = [self.T2(i) for i in range(N)]
        
        # Generate the system equations
        for i in range(N):
            eq = self.CreateEquation("Poisson_NonConstantStiffnexMatrix(%d)" % i)
            eq.Residual = numpy.sum(A[i, :] * T[:]) - F[i]
            
class simTutorial(daeSimulation):
    def __init__(self):
        daeSimulation.__init__(self)
        self.m = modTutorial("tutorial16")
        self.m.Description = __doc__
        
    def SetUpParametersAndDomains(self):
        self.m.x.CreateArray(4)
        self.m.L.SetValue(1)

    def SetUpVariables(self):
        pass
    
def run(**kwargs):
    simulation = simTutorial()
    return daeActivity.simulate(simulation, reportingInterval = 10, 
                                            timeHorizon       = 1000,
                                            **kwargs)

if __name__ == "__main__":
    guiRun = False if (len(sys.argv) > 1 and sys.argv[1] == 'console') else True
    run(guiRun = guiRun)