# This example implements homogenization of piezoeletric porous media.
# The mathematical model and numerical results are described in: 
#
# <NAME>., <NAME>.
# Homogenization of the fluid-saturated piezoelectric porous media.
# International Journal of Solids and Structures
# Volume 147, 15 August 2018, Pages 110-125
# https://doi.org/10.1016/j.ijsolstr.2018.05.017
#
# Run simulation:
#
#   ./simple.py example_poropiezo-1/poropiezo_macro_dfc.py
#
# The results are stored in `example_poropiezo-1/results` directory.
#

import sys
import numpy as nm
import os.path as osp
from sfepy import data_dir
from sfepy.base.base import Struct
from sfepy.homogenization.micmac import get_homog_coefs_linear
from sfepy.homogenization.recovery import recover_micro_hook_eps

data_dir = 'example_poropiezo-1'


def set_grad(ts, coors, mode=None, problem=None, **kwargs):
    if mode == 'qp':
        out = problem.data.reshape((coors.shape[0], 1, 1))
        return {'cs': out}


# projection of values from integration points into mesh vertices
def linear_projection(pb, cval):
    from sfepy.discrete import (FieldVariable, Material, Integral,
                                Equation, Equations, Problem)
    from sfepy.discrete.fem import Mesh, FEDomain, Field
    from sfepy.terms import Term
    from sfepy.solvers.ls import ScipyDirect
    from sfepy.solvers.nls import Newton
    from sfepy.base.base import IndexedStruct

    mesh = Mesh.from_file(pb.conf.filename_mesh)
    domain = FEDomain('domain', mesh)
    omega = domain.create_region('Omega', 'all')
    field = Field.from_args('scf', nm.float64, 'scalar', omega,
                            approx_order=1)

    g = FieldVariable('g', 'unknown', field)
    f = FieldVariable('f', 'test', field, primary_var_name='g')

    integral = Integral('i', order=2)
    m = Material('m', function=set_grad)

    t1 = Term.new('dw_volume_dot(f, g)', integral, omega, f=f, g=g)
    t2 = Term.new('dw_volume_lvf(m.cs, f)',
                  integral, omega, m=m, f=f)
    eq = Equation('balance', t1 - t2)
    eqs = 
Equations([eq])