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

"""
***********************************************************************************
tutorial_dealii_5.py
DAE Tools: pyDAE module, www.daetools.com
***********************************************************************************
DAE Tools is free software; you can redistribute it and/or modify it under the
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__ = """
In this example a simple flow through porous media is solved (deal.II step-20).

.. code-block:: none

K^{-1} u + nabla(p) = 0 in Omega
-div(u) = -f in Omega
p = g on dOmega

The mesh is a simple square:

.. image:: _static/square(-1,1)x(-1,1)-50x50.png
:width: 300 px

The velocity magnitude plot:

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

import os, sys, numpy, json, tempfile, random
from time import localtime, strftime
from daetools.pyDAE import *
from daetools.solvers.deal_II import *
from daetools.solvers.superlu import pySuperLU

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

class permeabilityFunction_2D(TensorFunction_2_2D):
def __init__(self, N):
TensorFunction_2_2D.__init__(self)

numpy.random.seed(1000)
self.centers = 2 * numpy.random.rand(N,2) - 1
# Create a Tensor<rank=2,dim=2> object to serve as a return value (to make the function faster)
self.inv_kappa = Tensor_2_2D()

def value(self, point, component = 0):
# 1) Sinusoidal (a function of the distance to the flowline)
#distance_to_flowline = numpy.fabs(point[1] - 0.2*numpy.sin(10*point[0]))
#permeability = numpy.exp(-(distance_to_flowline*distance_to_flowline)/0.01)
#norm_permeability = max(permeability, 0.001)

# 2) Random permeability field
x2 = numpy.square(point[0]-self.centers[:,0])
y2 = numpy.square(point[1]-self.centers[:,1])
permeability = numpy.sum( numpy.exp(-(x2 + y2) / 0.01) )
norm_permeability = max(permeability, 0.005)

# Set-up the inverse permeability tensor (only the diagonal items)
self.inv_kappa[0][0] = 1.0 / norm_permeability
self.inv_kappa[1][1] = 1.0 / norm_permeability

return self.inv_kappa

class pBoundaryFunction_2D(Function_2D):
def __init__(self, n_components = 1):
Function_2D.__init__(self, n_components)

def value(self, p, component = 0):
alpha = 0.3
beta  = 1.0
return -(alpha*p[0]*p[1]*p[1]/2.0 + beta*p[0] - alpha*p[0]*p[0]*p[0]/6.0)

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

dofs = [dealiiFiniteElementDOF_2D(name='u',
description='Velocity',
fe = FE_RaviartThomas_2D(0),
multiplicity=2),
dealiiFiniteElementDOF_2D(name='p',
description='Pressure',
fe = FE_DGQ_2D(0),
multiplicity=1)]
self.n_components = int(numpy.sum([dof.Multiplicity for dof in dofs]))

meshes_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'meshes')
mesh_file  = os.path.join(meshes_dir, 'square(-1,1)x(-1,1)-50x50.msh')

# Store the object so it does not go out of scope while still in use by daetools
self.fe_system = dealiiFiniteElementSystem_2D(meshFilename    = mesh_file,     # path to mesh
dofs            = dofs)          # degrees of freedom

self.fe_model = daeFiniteElementModel('PorousMedia', self, 'Flow through porous media', self.fe_system)

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

# deal.II Function<dim> wrapper.
self.fun_p_boundary = pBoundaryFunction_2D(self.n_components)
# deal.II TensorFunction<2,dim> wrapper.
self.fun_k_inverse  = permeabilityFunction_2D(40)

# Boundary IDs
left_edge   = 0
top_edge    = 1
right_edge  = 2
bottom_edge = 3

# Create some auxiliary objects for readability
phi_p_i         =  phi_2D('p', fe_i, fe_q)
phi_p_j         =  phi_2D('p', fe_j, fe_q)
dphi_p_i        = dphi_2D('p', fe_i, fe_q)
dphi_p_j        = dphi_2D('p', fe_j, fe_q)
phi_vector_u_i  =  phi_vector_2D('u', fe_i, fe_q)
phi_vector_u_j  =  phi_vector_2D('u', fe_j, fe_q)
dphi_vector_u_i = dphi_vector_2D('u', fe_i, fe_q)
dphi_vector_u_j = dphi_vector_2D('u', fe_j, fe_q)
div_phi_u_i     =     div_phi_2D('u', fe_i, fe_q)
div_phi_u_j     =     div_phi_2D('u', fe_j, fe_q)
normal = normal_2D(fe_q)
xyz    = xyz_2D(fe_q)
JxW    = JxW_2D(fe_q)

dirichletBC = {}

# Function value wrapper
p_boundary = function_value_2D("p_boundary", self.fun_p_boundary, xyz)
faceFi = {
left_edge:   -(phi_vector_u_i * normal) * p_boundary * JxW,
top_edge:    -(phi_vector_u_i * normal) * p_boundary * JxW,
right_edge:  -(phi_vector_u_i * normal) * p_boundary * JxW,
bottom_edge: -(phi_vector_u_i * normal) * p_boundary * JxW
}

# TensorFunction<2,dim>::value wrappers
k_inverse = tensor2_function_value_2D('k_inverse', self.fun_k_inverse, xyz)

# FE weak form terms
accumulation = 0.0 * JxW
velocity     = (phi_vector_u_i * k_inverse * phi_vector_u_j) * JxW
p_gradient   = -(div_phi_u_i * phi_p_j) * JxW
continuity   = -(phi_p_i * div_phi_u_j) * JxW
source       = 0.0 * JxW

weakForm = dealiiFiniteElementWeakForm_2D(Aij = velocity + p_gradient + continuity,
Mij = accumulation,
Fi  = source,
boundaryFaceFi  = faceFi,
functionsDirichletBC = dirichletBC)

# Setting the weak form of the FE system will declare a set of equations:
# [Mij]{dx/dt} + [Aij]{x} = {Fi} and boundary integral equations
self.fe_system.WeakForm = weakForm

class simTutorial(daeSimulation):
def __init__(self):
daeSimulation.__init__(self)
self.m = modTutorial("tutorial_dealii_5")
self.m.Description = __doc__
self.m.fe_model.Description = __doc__

def SetUpParametersAndDomains(self):
pass

def SetUpVariables(self):
pass

def run(**kwargs):
guiRun = kwargs.get('guiRun', False)

simulation = simTutorial()

# Create SuperLU LA solver
lasolver = pySuperLU.daeCreateSuperLUSolver()

# Create and setup two data reporters:
datareporter = daeDelegateDataReporter()
simName = simulation.m.Name + strftime(" [%d.%m.%Y %H:%M:%S]", localtime())
if guiRun:
results_folder = tempfile.mkdtemp(suffix = '-results', prefix = 'tutorial_deal_II_5-')
daeQtMessage("deal.II", "The simulation results will be located in: %s" % results_folder)
else:
results_folder = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'tutorial_deal_II_5-results')
print("The simulation results will be located in: %s" % results_folder)

# 1. deal.II (exports only FE DOFs in .vtk format to the specified directory)
feDataReporter = simulation.m.fe_system.CreateDataReporter()
if not feDataReporter.Connect(results_folder, simName):
sys.exit()

# 2. TCP/IP
tcpipDataReporter = daeTCPIPDataReporter()