r"""
Laplace equation with a field-dependent material parameter.

Find :math:`T(t)` for :math:`t \in [0, t_{\rm final}]` such that:

.. math::
   \int_{\Omega} c(T) \nabla s \cdot \nabla T
    = \int_{\Omega} f s
    \;, \quad \forall s \;.

where :math:`c(T)` is the :math:`T` dependent diffusion coefficient.
Each iteration calculates :math:`T` and adjusts :math:`c(T)`.
"""
from sfepy.base.base import output
from sfepy.mesh.mesh_generators import gen_block_mesh
from sfepy.discrete.fem.meshio import UserMeshIO
import numpy as np


def mesh_hook(mesh, mode):
    """
    Generate the block mesh.
    """
    if mode == 'read':
        mesh = gen_block_mesh(dims, shape, center, name='myPoissonBlock',
                              verbose=False)
        return mesh

    elif mode == 'write':
        pass

# Mesh dimensions
dims = [1, 1]
# Mesh resolution: increase to improve accuracy
shape = [21, 21]
# Center of mesh
center = [0, 0.5]

filename_mesh = UserMeshIO(mesh_hook) 
output_dir = './solutions/'

t0 = 0.0
t1 = 0.1
n_step = 2

materials = {
    'cond' : ({'val' :   2.0},), 
    'flux' : ({'val' : -10.0},),
    'G'    : ({'val' :  20.0},),
    'insulated' : ({'val' : 0.0},),
}

fields = {
    'temperature' : ('real', 1, 'Omega', 1),
}

variables = {
    'T' : ('unknown field', 'temperature', 0),
    's' : ('test field',    'temperature', 'T'),
}

regions = {
    'Omega'  : 'all',
    'Left'   : ('vertices in (x < -0.499999)', 'facet'),
    'Right'  : ('vertices in (x >  0.499999)', 'facet'),
    'Bottom' : ('vertices in (y <  0.000001)', 'facet'),
    'Top'    : ('vertices in (y >  0.999999)', 'facet'),
}

ebcs = {
    'bc_l' : ('Left',   {'T.0' :  2.0}),
    'bc_r' : ('Right',  {'T.0' :  2.0}),
}

integrals = {
    'i' : 1,
}

equations = {
    'Temperature' : """
           dw_laplace.i.Omega(cond.val, s, T)
         - dw_surface_integrate.2.Top(insulated.val, s)
         - dw_volume_lvf.2.Omega(G.val, s)
"""
}

solvers = {
    'ls' : ('ls.scipy_direct', {}),
    'newton' : ('nls.newton', {
        'i_max' : 1,
        'eps_a' : 1e-10,
        'eps_r' : 1.0,
    }),
    'ts' : ('ts.simple', {
        't0' : t0,
        't1' : t1,
        'dt' : None,
        'n_step' : n_step, # has precedence over dt!
        'quasistatic' : True,
    }),
}

options = {
    'nls' : 'newton',
    'ls' : 'ls',
    'ts' : 'ts',
    'save_steps' : -1,
    'output_dir' : output_dir,
}