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

"""
***********************************************************************************
tutorial_dealii_3.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 the Cahn-Hilliard equation is solved using the finite element method.
This equation describes the process of phase separation, where two components of a
binary mixture separate and form domains pure in each component.

.. code-block:: none

dc/dt - D*nabla^2(mu) = 0, in Omega
mu = c^3 - c - gamma*nabla^2(c)

The mesh is a simple square (0-100)x(0-100):

.. image:: _static/square.png
:width: 300 px

The concentration plot at t = 500s:

.. image:: _static/tutorial_dealii_3-results.png
:height: 400 px

Animation:

.. image:: _static/tutorial_dealii_3-animation.gif
:height: 400 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 modTutorial(daeModel):
def __init__(self, Name, Parent = None, Description = ""):
daeModel.__init__(self, Name, Parent, Description)

dofs = [dealiiFiniteElementDOF_2D(name='c',
description='Concentration',
fe = FE_Q_2D(1),
multiplicity=1),
dealiiFiniteElementDOF_2D(name='mu',
description='Chemical potential',
fe = FE_Q_2D(1),
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')
# This mesh is a coarse one (30x30 cells); there are finer meshes available: 50x50, 100x100, 200x200
mesh_file  = os.path.join(meshes_dir, 'square(0,100)x(0,100)-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('CahnHilliard', self, 'Cahn-Hilliard equation', self.fe_system)

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

# Create some auxiliary objects for readability
phi_c_i  =  phi_2D('c', fe_i, fe_q)
phi_c_j  =  phi_2D('c', fe_j, fe_q)
dphi_c_i = dphi_2D('c', fe_i, fe_q)
dphi_c_j = dphi_2D('c', fe_j, fe_q)

phi_mu_i  =  phi_2D('mu', fe_i, fe_q)
phi_mu_j  =  phi_2D('mu', fe_j, fe_q)
dphi_mu_i = dphi_2D('mu', fe_i, fe_q)
dphi_mu_j = dphi_2D('mu', fe_j, fe_q)

# FE approximation of a quantity at the specified quadrature point (adouble object)
c      = dof_approximation_2D('c', fe_q)
normal = normal_2D(fe_q)
xyz    = xyz_2D(fe_q)
JxW    = JxW_2D(fe_q)

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

dirichletBC = {}
surfaceIntegrals = {}

self.useWikipedia_fc = True

# 1) f(c) from the Wikipedia (https://en.wikipedia.org/wiki/Cahn-Hilliard_equation)
if self.useWikipedia_fc:
Diffusivity = 1.0
gamma       = 1.0
def f(c):
return c**3 - c

# 2) f(c) used by Raymond Smith (M.Z.Bazant's group, MIT) for phase-separating battery electrodes
if not self.useWikipedia_fc:
Diffusivity = 1
gamma       = 1
Omg_a       = 3.4
log_fe = feExpression_2D.log
def f(c):
# The original expression is:
#   log_fe(c/(1-c)) + Omg_a*(1-2*c)
# However, the one below is much more computationally efficient and requires less memory,
# since a Finite Element approximation of a DoF is an expensive operation:
#   sum(phi_j * dof(j))
# For vector-valued DoFs it is even more demanding for the approximation is:
#   sum(phi_vector_j * dof(j))
# where phi_vector_j is a rank=1 Tensor.
return log_fe(1 + 1/(1-c)) + Omg_a - (2*Omg_a)*c

# FE weak form terms
c_accumulation    = (phi_c_i * phi_c_j) * JxW
mu_diffusion_c_eq = (dphi_c_i * dphi_mu_j) * Diffusivity * JxW
mu                = phi_mu_i *  phi_mu_j * JxW
c_diffusion_mu_eq = (-dphi_mu_i * dphi_c_j) * gamma * JxW
fun_c             = (phi_mu_i * JxW) * f(c)

cell_Aij = mu_diffusion_c_eq + mu + c_diffusion_mu_eq + c_diffusion_mu_eq
cell_Mij = c_accumulation
cell_Fi  = fun_c

weakForm = dealiiFiniteElementWeakForm_2D(Aij = cell_Aij,
Mij = cell_Mij,
Fi  = cell_Fi,
functionsDirichletBC = dirichletBC,
surfaceIntegrals = surfaceIntegrals)

print('Cahn-Hilliard equation:')
print('    Aij = %s' % str(cell_Aij))
print('    Mij = %s' % str(cell_Mij))
print('    Fi  = %s' % str(cell_Fi))
print('    surfaceIntegrals  = %s' % str([item for item in surfaceIntegrals]))

# 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_3")
self.m.Description = __doc__
self.m.fe_model.Description = __doc__

def SetUpParametersAndDomains(self):
pass

def SetUpVariables(self):
numpy.random.seed(124)

def c_with_noise_wiki(index, overallIndex):
c0     = 0.0
stddev = 0.1
return numpy.random.normal(c0, stddev)

def c_with_noise_ray(index, overallIndex):
c0     = 0.5
stddev = 0.1
return numpy.random.normal(c0, stddev)

if self.m.useWikipedia_fc:
setFEInitialConditions(self.m.fe_model, self.m.fe_system, 'c', c_with_noise_wiki)
else:
setFEInitialConditions(self.m.fe_model, self.m.fe_system, 'c', c_with_noise_ray)

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_3-')
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_3-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()