def test_projection_iga_fem(self):
        from sfepy.discrete import FieldVariable
        from sfepy.discrete.fem import FEDomain, Field
        from sfepy.discrete.iga.domain import IGDomain
        from sfepy.mesh.mesh_generators import gen_block_mesh
        from sfepy.discrete.iga.domain_generators import gen_patch_block_domain
        from sfepy.discrete.projections import (make_l2_projection,
                                                make_l2_projection_data)

        shape = [10, 12, 12]
        dims = [5, 6, 6]
        centre = [0, 0, 0]
        degrees = [2, 2, 2]

        nurbs, bmesh, regions = gen_patch_block_domain(dims, shape, centre,
                                                       degrees,
                                                       cp_mode='greville',
                                                       name='iga')
        ig_domain = IGDomain('iga', nurbs, bmesh, regions=regions)

        ig_omega = ig_domain.create_region('Omega', 'all')
        ig_field = Field.from_args('iga', nm.float64, 1, ig_omega,
                                   approx_order='iga', poly_space_base='iga')
        ig_u = FieldVariable('ig_u', 'parameter', ig_field,
                             primary_var_name='(set-to-None)')

        mesh = gen_block_mesh(dims, shape, centre, name='fem')
        fe_domain = FEDomain('fem', mesh)

        fe_omega = fe_domain.create_region('Omega', 'all')
        fe_field = Field.from_args('fem', nm.float64, 1, fe_omega,
                                   approx_order=2)
        fe_u = FieldVariable('fe_u', 'parameter', fe_field,
                             primary_var_name='(set-to-None)')

        def _eval_data(ts, coors, mode, **kwargs):
            return nm.prod(coors**2, axis=1)[:, None, None]

        make_l2_projection_data(ig_u, _eval_data)

        make_l2_projection(fe_u, ig_u) # This calls ig_u.evaluate_at().

        coors = 0.5 * nm.random.rand(20, 3) * dims

        ig_vals = ig_u.evaluate_at(coors)
        fe_vals = fe_u.evaluate_at(coors)

        ok = nm.allclose(ig_vals, fe_vals, rtol=0.0, atol=1e-12)
        if not ok:
            self.report('iga-fem projection failed!')
            self.report('coors:')
            self.report(coors)
            self.report('iga fem diff:')
            self.report(nm.c_[ig_vals, fe_vals, nm.abs(ig_vals - fe_vals)])

        return ok

    def test_interpolation_two_meshes(self):
        from sfepy import data_dir
        from sfepy.discrete import Variables
        from sfepy.discrete.fem import Mesh, FEDomain, Field

        m1 = Mesh.from_file(data_dir + '/meshes/3d/block.mesh')

        m2 = Mesh.from_file(data_dir + '/meshes/3d/cube_medium_tetra.mesh')
        m2.coors[:] *= 2.0

        bbox = m1.get_bounding_box()
        dd = bbox[1,:] - bbox[0,:]
        data = nm.sin(4.0 * nm.pi * m1.coors[:,0:1] / dd[0]) \
               * nm.cos(4.0 * nm.pi * m1.coors[:,1:2] / dd[1])

        variables1 = {
            'u'       : ('unknown field', 'scalar_tp', 0),
            'v'       : ('test field',    'scalar_tp', 'u'),
        }

        variables2 = {
            'u'       : ('unknown field', 'scalar_si', 0),
            'v'       : ('test field',    'scalar_si', 'u'),
        }

        d1 = FEDomain('d1', m1)
        omega1 = d1.create_region('Omega', 'all')
        field1 = Field.from_args('scalar_tp', nm.float64, (1,1), omega1,
                                 approx_order=1)
        ff1 = {field1.name : field1}

        d2 = FEDomain('d2', m2)
        omega2 = d2.create_region('Omega', 'all')
        field2 = Field.from_args('scalar_si', nm.float64, (1,1), omega2,
                                 approx_order=0)
        ff2 = {field2.name : field2}

        vv1 = Variables.from_conf(transform_variables(variables1), ff1)
        u1 = vv1['u']
        u1.set_from_mesh_vertices(data)

        vv2 = Variables.from_conf(transform_variables(variables2), ff2)
        u2 = vv2['u']

        # Performs interpolation, if other field differs from self.field
        # or, in particular, is defined on a different mesh.
        u2.set_from_other(u1, strategy='interpolation', close_limit=0.1)

        fname = in_dir(self.options.out_dir)
        u1.save_as_mesh(fname('test_mesh_interp_block_scalar.vtk'))
        u2.save_as_mesh(fname('test_mesh_interp_cube_scalar.vtk'))

        return True

    def from_conf(conf, options):
        import sfepy
        from sfepy.discrete.fem import Mesh, Domain, Field
        mesh = Mesh.from_file('meshes/2d/rectangle_tri.mesh',
                              prefix_dir=sfepy.data_dir)
        domain = Domain('domain', mesh)
        dim = domain.shape.dim

        min_x, max_x = domain.get_mesh_bounding_box()[:,0]
        eps = 1e-8 * (max_x - min_x)

        omega = domain.create_region('Omega', 'all')
        gamma1 = domain.create_region('Gamma1',
                                      'vertices in x < %.10f' % (min_x + eps),
                                      'facet')
        gamma2 = domain.create_region('Gamma2',
                                      'vertices in x > %.10f' % (max_x - eps),
                                      'facet')

        field = Field.from_args('fu', nm.float64, 'vector', omega,
                                approx_order=2)

        test = Test(conf=conf, options=options, dim=dim,
                    omega=omega, gamma1=gamma1, gamma2=gamma2,
                    field=field)
        return test

def nodal_stress(out, pb, state, extend=False, integrals=None):
    """
    Calculate stresses at nodal points.
    """

    # Point load.
    mat = pb.get_materials()["Load"]
    P = 2.0 * mat.get_data("special", "val")[1]

    # Calculate nodal stress.
    pb.time_update()

    if integrals is None:
        integrals = pb.get_integrals()

    stress = pb.evaluate("ev_cauchy_stress.ivn.Omega(Asphalt.D, u)", mode="qp", integrals=integrals)
    sfield = Field.from_args("stress", nm.float64, (3,), pb.domain.regions["Omega"])
    svar = FieldVariable("sigma", "parameter", sfield, primary_var_name="(set-to-None)")
    svar.set_data_from_qp(stress, integrals["ivn"])

    print "\n=================================================================="
    print "Given load = %.2f N" % -P
    print "\nAnalytical solution"
    print "==================="
    print "Horizontal tensile stress = %.5e MPa/mm" % (-2.0 * P / (nm.pi * 150.0))
    print "Vertical compressive stress = %.5e MPa/mm" % (-6.0 * P / (nm.pi * 150.0))
    print "\nFEM solution"
    print "============"
    print "Horizontal tensile stress = %.5e MPa/mm" % (svar()[0][0])
    print "Vertical compressive stress = %.5e MPa/mm" % (-svar()[0][1])
    print "=================================================================="
    return out

def _get_bqp(geometry, order):
    from sfepy.discrete import Integral
    from sfepy.discrete.fem.geometry_element import GeometryElement
    from sfepy.discrete.fem import Mesh, FEDomain, Field

    gel = GeometryElement(geometry)

    mesh = Mesh.from_data('aux', gel.coors, None,
                          [gel.conn[None, :]], [[0]], [geometry])
    domain = FEDomain('domain', mesh)
    omega = domain.create_region('Omega', 'all')
    surf =  domain.create_region('Surf', 'vertices of surface', 'facet')
    field = Field.from_args('f', nm.float64, shape=1,
                            region=omega, approx_order=1)
    field.setup_surface_data(surf)

    integral = Integral('aux', order=order)
    field.create_bqp('Surf', integral)

    sd = field.surface_data['Surf']
    qp = field.qp_coors[(integral.order, sd.bkey)]

    output('geometry:', geometry, 'order:', order, 'num. points:',
           qp.vals.shape[1], 'true_order:',
           integral.qps[gel.surface_facet_name].order)
    output('min. weight:', qp.weights.min())
    output('max. weight:', qp.weights.max())

    return (gel, qp.vals.reshape((-1, mesh.dim)),
            nm.tile(qp.weights, qp.vals.shape[0]))

def nodal_stress(out, pb, state, extend=False, integrals=None):
    '''
    Calculate stresses at nodal points.
    '''

    # Point load.
    mat = pb.get_materials()['Load']
    P = 2.0 * mat.get_data('special', 'val')[1]

    # Calculate nodal stress.
    pb.time_update()

    if integrals is None: integrals = pb.get_integrals()

    stress = pb.evaluate('ev_cauchy_stress.ivn.Omega(Asphalt.D, u)', mode='qp',
                         integrals=integrals, copy_materials=False)
    sfield = Field.from_args('stress', nm.float64, (3,),
                             pb.domain.regions['Omega'])
    svar = FieldVariable('sigma', 'parameter', sfield,
                         primary_var_name='(set-to-None)')
    svar.set_from_qp(stress, integrals['ivn'])

    print('\n==================================================================')
    print('Given load = %.2f N' % -P)
    print('\nAnalytical solution')
    print('===================')
    print('Horizontal tensile stress = %.5e MPa/mm' % (-2.*P/(nm.pi*150.)))
    print('Vertical compressive stress = %.5e MPa/mm' % (-6.*P/(nm.pi*150.)))
    print('\nFEM solution')
    print('============')
    print('Horizontal tensile stress = %.5e MPa/mm' % (svar()[0]))
    print('Vertical compressive stress = %.5e MPa/mm' % (-svar()[1]))
    print('==================================================================')
    return out

def do_interpolation(m2, m1, data, field_name, force=False):
    """Interpolate data from m1 to m2. """
    from sfepy.discrete import Variables
    from sfepy.discrete.fem import FEDomain, Field

    fields = {
        'scalar_si' : ((1,1), 'Omega', 2),
        'vector_si' : ((3,1), 'Omega', 2),
        'scalar_tp' : ((1,1), 'Omega', 1),
        'vector_tp' : ((3,1), 'Omega', 1),
    }

    d1 = FEDomain('d1', m1)

    omega1 = d1.create_region('Omega', 'all')

    f = fields[field_name]

    field1 = Field.from_args('f', nm.float64, f[0], d1.regions[f[1]],
                             approx_order=f[2])
    ff = {field1.name : field1}

    vv = Variables.from_conf(transform_variables(variables), ff)
    u1 = vv['u']
    u1.set_from_mesh_vertices(data)

    d2 = FEDomain('d2', m2)
    omega2 = d2.create_region('Omega', 'all')

    field2 = Field.from_args('f', nm.float64, f[0], d2.regions[f[1]],
                             approx_order=f[2])
    ff2 = {field2.name : field2}

    vv2 = Variables.from_conf(transform_variables(variables), ff2)
    u2 = vv2['u']

    if not force:
        # Performs interpolation, if other field differs from self.field
        # or, in particular, is defined on a different mesh.
        u2.set_from_other(u1, strategy='interpolation', close_limit=0.5)

    else:
        coors = u2.field.get_coor()
        vals = u1.evaluate_at(coors, close_limit=0.5)
        u2.set_data(vals)

    return u1, u2

def run(domain, order):
    omega = domain.create_region('Omega', 'all')
    bbox = domain.get_mesh_bounding_box()
    min_x, max_x = bbox[:, 0]
    min_y, max_y = bbox[:, 1]
    eps = 1e-8 * (max_x - min_x)
    gamma1 = domain.create_region('Gamma1',
                                  'vertices in (x < %.10f)' % (min_x + eps),
                                  'facet')
    gamma2 = domain.create_region('Gamma2',
                                  'vertices in (x > %.10f)' % (max_x - eps),
                                  'facet')
    gamma3 = domain.create_region('Gamma3',
                                  'vertices in y < %.10f' % (min_y + eps),
                                  'facet')
    gamma4 = domain.create_region('Gamma4',
                                  'vertices in y > %.10f' % (max_y - eps),
                                  'facet')

    field = Field.from_args('fu', nm.float64, 1, omega, approx_order=order)

    u = FieldVariable('u', 'unknown', field)
    v = FieldVariable('v', 'test', field, primary_var_name='u')

    integral = Integral('i', order=2*order)

    t1 = Term.new('dw_laplace(v, u)',
                  integral, omega, v=v, u=u)
    eq = Equation('eq', t1)
    eqs = Equations([eq])

    fix1 = EssentialBC('fix1', gamma1, {'u.0' : 0.4})
    fix2 = EssentialBC('fix2', gamma2, {'u.0' : 0.0})

    def get_shift(ts, coors, region):
        return nm.ones_like(coors[:, 0])

    dof_map_fun = Function('dof_map_fun', per.match_x_line)
    shift_fun = Function('shift_fun', get_shift)

    sper = LinearCombinationBC('sper', [gamma3, gamma4], {'u.0' : 'u.0'},
                               dof_map_fun, 'shifted_periodic',
                               arguments=(shift_fun,))

    ls = ScipyDirect({})

    pb = Problem('laplace', equations=eqs, auto_solvers=None)

    pb.time_update(ebcs=Conditions([fix1, fix2]), lcbcs=Conditions([sper]))

    ev = pb.get_evaluator()
    nls = Newton({}, lin_solver=ls,
                 fun=ev.eval_residual, fun_grad=ev.eval_tangent_matrix)

    pb.set_solver(nls)

    state = pb.solve()

    return pb, state

def main():
    from sfepy import data_dir

    parser = OptionParser(usage=usage, version="%prog")
    parser.add_option("-s", "--show", action="store_true", dest="show", default=False, help=help["show"])
    options, args = parser.parse_args()

    mesh = Mesh.from_file(data_dir + "/meshes/2d/rectangle_tri.mesh")
    domain = Domain("domain", mesh)

    min_x, max_x = domain.get_mesh_bounding_box()[:, 0]
    eps = 1e-8 * (max_x - min_x)
    omega = domain.create_region("Omega", "all")
    gamma1 = domain.create_region("Gamma1", "vertices in x < %.10f" % (min_x + eps), "facet")
    gamma2 = domain.create_region("Gamma2", "vertices in x > %.10f" % (max_x - eps), "facet")

    field = Field.from_args("fu", nm.float64, "vector", omega, approx_order=2)

    u = FieldVariable("u", "unknown", field)
    v = FieldVariable("v", "test", field, primary_var_name="u")

    m = Material("m", lam=1.0, mu=1.0)
    f = Material("f", val=[[0.02], [0.01]])

    integral = Integral("i", order=3)

    t1 = Term.new("dw_lin_elastic_iso(m.lam, m.mu, v, u)", integral, omega, m=m, v=v, u=u)
    t2 = Term.new("dw_volume_lvf(f.val, v)", integral, omega, f=f, v=v)
    eq = Equation("balance", t1 + t2)
    eqs = Equations([eq])

    fix_u = EssentialBC("fix_u", gamma1, {"u.all": 0.0})

    bc_fun = Function("shift_u_fun", shift_u_fun, extra_args={"shift": 0.01})
    shift_u = EssentialBC("shift_u", gamma2, {"u.0": bc_fun})

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = Problem("elasticity", equations=eqs, nls=nls, ls=ls)
    pb.save_regions_as_groups("regions")

    pb.time_update(ebcs=Conditions([fix_u, shift_u]))

    vec = pb.solve()
    print nls_status

    pb.save_state("linear_elasticity.vtk", vec)

    if options.show:
        view = Viewer("linear_elasticity.vtk")
        view(vector_mode="warp_norm", rel_scaling=2, is_scalar_bar=True, is_wireframe=True)

def solve_problem(shape, dims, young, poisson, force, transform=None):
    domain = make_domain(dims[:2], shape, transform=transform)

    omega = domain.regions['Omega']
    gamma1 = domain.regions['Gamma1']
    gamma2 = domain.regions['Gamma2']

    field = Field.from_args('fu', nm.float64, 6, omega, approx_order=1,
                            poly_space_base='shell10x')
    u = FieldVariable('u', 'unknown', field)
    v = FieldVariable('v', 'test', field, primary_var_name='u')

    thickness = dims[2]
    if transform is None:
        pload = [[0.0, 0.0, force / shape[1], 0.0, 0.0, 0.0]] * shape[1]

    elif transform == 'bend':
        pload = [[force / shape[1], 0.0, 0.0, 0.0, 0.0, 0.0]] * shape[1]

    elif transform == 'twist':
        pload = [[0.0, force / shape[1], 0.0, 0.0, 0.0, 0.0]] * shape[1]

    m = Material('m', D=sh.create_elastic_tensor(young=young, poisson=poisson),
                 values={'.drill' : 1e-7})
    load = Material('load', values={'.val' : pload})

    aux = Integral('i', order=3)
    qp_coors, qp_weights = aux.get_qp('3_8')
    qp_coors[:, 2] = thickness * (qp_coors[:, 2] - 0.5)
    qp_weights *= thickness

    integral = Integral('i', coors=qp_coors, weights=qp_weights, order='custom')

    t1 = Term.new('dw_shell10x(m.D, m.drill, v, u)',
                  integral, omega, m=m, v=v, u=u)
    t2 = Term.new('dw_point_load(load.val, v)',
                  integral, gamma2, load=load, v=v)
    eq = Equation('balance', t1 - t2)
    eqs = Equations([eq])

    fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = Problem('elasticity with shell10x', equations=eqs)
    pb.set_bcs(ebcs=Conditions([fix_u]))
    pb.set_solver(nls)

    state = pb.solve()

    return pb, state, u, gamma2

    def from_conf(conf, options):
        from sfepy.discrete import FieldVariable, Variables, Problem
        from sfepy.discrete.fem import Mesh, FEDomain, Field

        mesh = Mesh.from_file(data_dir + '/meshes/2d/square_unit_tri.mesh')
        domain = FEDomain('domain', mesh)

        omega = domain.create_region('Omega', 'all')
        domain.create_region('Left',
                             'vertices in (x < -0.499)',
                             'facet')
        domain.create_region('LeftStrip',
                             'vertices in (x < -0.499)'
                             ' & (y > -0.199) & (y < 0.199)',
                             'facet')
        domain.create_region('LeftFix',
                             'r.Left -v r.LeftStrip',
                             'facet')
        domain.create_region('Right',
                             'vertices in (x > 0.499)',
                             'facet')
        domain.create_region('RightStrip',
                             'vertices in (x > 0.499)'
                             ' & (y > -0.199) & (y < 0.199)',
                             'facet')
        domain.create_region('RightFix',
                             'r.Right -v r.RightStrip',
                             'facet')

        fu = Field.from_args('fu', nm.float64, 'vector', omega, approx_order=2)
        u = FieldVariable('u', 'unknown', fu)

        fp = Field.from_args('fp', nm.float64, 'scalar', omega, approx_order=2)
        p = FieldVariable('p', 'unknown', fp)

        pb = Problem('test', domain=domain, fields=[fu, fp],
                     auto_conf=False, auto_solvers=False)

        test = Test(problem=pb, variables=Variables([u, p]),
                    conf=conf, options=options)
        return test

    def from_conf(conf, options):
        mesh = Mesh.from_file('meshes/2d/square_unit_tri.mesh',
                              prefix_dir=sfepy.data_dir)
        domain = FEDomain('domain', mesh)

        omega = domain.create_region('Omega', 'all')

        field = Field.from_args('linear', nm.float64, 'scalar', omega,
                                approx_order=1)

        test = Test(conf=conf, options=options, omega=omega, field=field)
        return test

    def test_evaluate_at(self):
        from sfepy import data_dir
        from sfepy.discrete.fem import Mesh
        from sfepy.discrete import Variables
        from sfepy.discrete.fem import FEDomain, Field

        meshes = {
            'tp' : Mesh.from_file(data_dir + '/meshes/3d/block.mesh'),
        }
        datas = gen_datas(meshes)

        fields = {
            'scalar_tp' : ((1,1), 'Omega', 1),
            'vector_tp' : ((3,1), 'Omega', 1),
        }

        ok = True
        for field_name in ['scalar_tp', 'vector_tp']:
            d = FEDomain('d', meshes['tp'])
            d.create_region('Omega', 'all')

            f = fields[field_name]
            field = Field.from_args('f', nm.complex128, f[0],
                                    d.regions[f[1]],
                                    approx_order=f[2])
            ff = {field.name : field}

            vv = Variables.from_conf(transform_variables(variables), ff)
            u = vv['u']

            bbox = d.get_mesh_bounding_box()
            t = nm.expand_dims(nm.linspace(0, 1, 100), 1)
            coors = nm.expand_dims(bbox[1] - bbox[0], 0) * t + bbox[0]

            data_r = datas[field_name]
            data_i = 2. / (1 + datas[field_name])

            u.set_from_mesh_vertices(data_r)
            vals_r = u.evaluate_at(coors)
            u.set_from_mesh_vertices(data_i)
            vals_i = u.evaluate_at(coors)
            u.set_from_mesh_vertices(data_r + data_i * 1j)
            vals = u.evaluate_at(coors)

            _ok = nm.allclose(vals_r + vals_i * 1j, vals, rtol=0.0, atol=1e-12)
            _ok = _ok and nm.abs(vals).sum() > 1
            self.report('evaluating complex field %s: %s' % (field_name, _ok))

            ok = ok and _ok

        return ok

def run(domain, order):
    omega = domain.create_region("Omega", "all")
    bbox = domain.get_mesh_bounding_box()
    min_x, max_x = bbox[:, 0]
    min_y, max_y = bbox[:, 1]
    eps = 1e-8 * (max_x - min_x)
    gamma1 = domain.create_region("Gamma1", "vertices in (x < %.10f)" % (min_x + eps), "facet")
    gamma2 = domain.create_region("Gamma2", "vertices in (x > %.10f)" % (max_x - eps), "facet")
    gamma3 = domain.create_region("Gamma3", "vertices in y < %.10f" % (min_y + eps), "facet")
    gamma4 = domain.create_region("Gamma4", "vertices in y > %.10f" % (max_y - eps), "facet")

    field = Field.from_args("fu", nm.float64, 1, omega, approx_order=order)

    u = FieldVariable("u", "unknown", field)
    v = FieldVariable("v", "test", field, primary_var_name="u")

    integral = Integral("i", order=2 * order)

    t1 = Term.new("dw_laplace(v, u)", integral, omega, v=v, u=u)
    eq = Equation("eq", t1)
    eqs = Equations([eq])

    fix1 = EssentialBC("fix1", gamma1, {"u.0": 0.4})
    fix2 = EssentialBC("fix2", gamma2, {"u.0": 0.0})

    def get_shift(ts, coors, region):
        return nm.ones_like(coors[:, 0])

    dof_map_fun = Function("dof_map_fun", per.match_x_line)
    shift_fun = Function("shift_fun", get_shift)

    sper = LinearCombinationBC(
        "sper", [gamma3, gamma4], {"u.0": "u.0"}, dof_map_fun, "shifted_periodic", arguments=(shift_fun,)
    )

    ls = ScipyDirect({})

    pb = Problem("laplace", equations=eqs, auto_solvers=None)

    pb.time_update(ebcs=Conditions([fix1, fix2]), lcbcs=Conditions([sper]))

    ev = pb.get_evaluator()
    nls = Newton({}, lin_solver=ls, fun=ev.eval_residual, fun_grad=ev.eval_tangent_matrix)

    pb.set_solver(nls)

    state = pb.solve()

    return pb, state

    def test_invariance_qp(self):
        from sfepy import data_dir
        from sfepy.discrete import Variables, Integral
        from sfepy.discrete.fem import Mesh, FEDomain, Field
        from sfepy.terms import Term
        from sfepy.discrete.common.mappings import get_physical_qps

        mesh = Mesh.from_file(data_dir + '/meshes/3d/block.mesh')

        bbox = mesh.get_bounding_box()
        dd = bbox[1,:] - bbox[0,:]
        data = nm.sin(4.0 * nm.pi * mesh.coors[:,0:1] / dd[0]) \
               * nm.cos(4.0 * nm.pi * mesh.coors[:,1:2] / dd[1])

        variables = {
            'u'       : ('unknown field', 'scalar_tp', 0),
            'v'       : ('test field',    'scalar_tp', 'u'),
        }

        domain = FEDomain('domain', mesh)
        omega = domain.create_region('Omega', 'all')
        field = Field.from_args('scalar_tp', nm.float64, 1, omega,
                                approx_order=1)
        ff = {field.name : field}

        vv = Variables.from_conf(transform_variables(variables), ff)
        u = vv['u']
        u.set_from_mesh_vertices(data)

        integral = Integral('i', order=2)
        term = Term.new('ev_volume_integrate(u)', integral, omega, u=u)
        term.setup()
        val1 = term.evaluate(mode='qp')
        val1 = val1.ravel()

        qps = get_physical_qps(omega, integral)
        coors = qps.values

        val2 = u.evaluate_at(coors).ravel()

        self.report('max. difference:', nm.abs(val1 - val2).max())
        ok = nm.allclose(val1, val2, rtol=0.0, atol=1e-12)
        self.report('invariance in qp: %s' % ok)

        return ok

def prepare_variable(filename, n_components):
    from sfepy.discrete import FieldVariable
    from sfepy.discrete.fem import Mesh, FEDomain, Field

    mesh = Mesh.from_file(filename)

    bbox = mesh.get_bounding_box()
    dd = bbox[1, :] - bbox[0, :]
    data = nm.sin(4.0 * nm.pi * mesh.coors[:, 0:1] / dd[0]) * nm.cos(4.0 * nm.pi * mesh.coors[:, 1:2] / dd[1])

    domain = FEDomain("domain", mesh)
    omega = domain.create_region("Omega", "all")
    field = Field.from_args("field", nm.float64, n_components, omega, approx_order=2)

    u = FieldVariable("u", "parameter", field, primary_var_name="(set-to-None)")
    u.set_from_mesh_vertices(data * nm.arange(1, n_components + 1)[None, :])

    return u

    def test_projection_tri_quad(self):
        from sfepy.discrete.projections import make_l2_projection

        source = FieldVariable('us', 'unknown', self.field)

        coors = self.field.get_coor()
        vals = nm.sin(2.0 * nm.pi * coors[:,0] * coors[:,1])
        source.set_data(vals)

        name = op.join(self.options.out_dir,
                       'test_projection_tri_quad_source.vtk')
        source.save_as_mesh(name)

        mesh = Mesh.from_file('meshes/2d/square_quad.mesh',
                              prefix_dir=sfepy.data_dir)
        domain = FEDomain('domain', mesh)

        omega = domain.create_region('Omega', 'all')


        field = Field.from_args('bilinear', nm.float64, 'scalar', omega,
                                approx_order=1)

        target = FieldVariable('ut', 'unknown', field)

        make_l2_projection(target, source)

        name = op.join(self.options.out_dir,
                       'test_projection_tri_quad_target.vtk')
        target.save_as_mesh(name)

        bbox = self.field.domain.get_mesh_bounding_box()
        x = nm.linspace(bbox[0, 0] + 0.001, bbox[1, 0] - 0.001, 20)
        y = nm.linspace(bbox[0, 1] + 0.001, bbox[1, 1] - 0.001, 20)

        xx, yy = nm.meshgrid(x, y)
        test_coors = nm.c_[xx.ravel(), yy.ravel()].copy()

        vec1 = source.evaluate_at(test_coors)
        vec2 = target.evaluate_at(test_coors)

        ok = (nm.abs(vec1 - vec2) < 0.01).all()

        return ok

def prepare_variable(filename, n_components):
    from sfepy.discrete import FieldVariable
    from sfepy.discrete.fem import Mesh, FEDomain, Field

    mesh = Mesh.from_file(filename)

    bbox = mesh.get_bounding_box()
    dd = bbox[1,:] - bbox[0,:]
    data = (nm.sin(4.0 * nm.pi * mesh.coors[:,0:1] / dd[0])
            * nm.cos(4.0 * nm.pi * mesh.coors[:,1:2] / dd[1]))

    domain = FEDomain('domain', mesh)
    omega = domain.create_region('Omega', 'all')
    field = Field.from_args('field', nm.float64, n_components, omega,
                            approx_order=2)

    u = FieldVariable('u', 'parameter', field,
                      primary_var_name='(set-to-None)')
    u.set_from_mesh_vertices(nm.c_[tuple([data] * n_components)])

    return u

    def test_project_tensors(self):
        from sfepy.discrete import FieldVariable
        from sfepy.discrete.projections import project_by_component

        ok = True

        u = FieldVariable('u', 'parameter', self.field,
                          primary_var_name='(set-to-None)')
        u.set_constant(1.0)

        component = FieldVariable('component', 'parameter', self.field,
                                  primary_var_name='(set-to-None)')

        nls_options = {'eps_a' : 1e-16, 'i_max' : 1}

        u_qp = u.evaluate()
        u2 = FieldVariable('u2', 'parameter', self.field,
                           primary_var_name='(set-to-None)')
        project_by_component(u2, u_qp, component, self.field.approx_order,
                             nls_options=nls_options)

        _ok = self.compare_vectors(u(), u2())
        ok = ok and _ok

        gu_qp = u.evaluate(mode='grad')

        gfield = Field.from_args('gu', nm.float64, 2, self.field.region,
                                 approx_order=self.field.approx_order)
        gu = FieldVariable('gu', 'parameter', gfield,
                           primary_var_name='(set-to-None)')

        project_by_component(gu, gu_qp, component, gfield.approx_order,
                             nls_options=nls_options)

        _ok = self.compare_vectors(gu(), nm.zeros_like(gu()))
        ok = ok and _ok

        return ok

    def test_linearization(self):
        from sfepy.base.base import Struct
        from sfepy.discrete.fem import Mesh, FEDomain, Field
        from sfepy import data_dir

        geometries = ['2_3', '2_4', '3_4', '3_8']
        approx_orders = [1, 2]
        funs = [nm.cos, nm.sin, lambda x: x]

        ok = True
        for geometry in geometries:
            name = os.path.join(data_dir,
                                'meshes/elements/%s_1.mesh' % geometry)
            mesh = Mesh.from_file(name)

            domain = FEDomain('', mesh)
            domain = domain.refine()

            domain.mesh.write(self.join('linearizer-%s-0.mesh' % geometry))

            omega = domain.create_region('Omega', 'all')

            for approx_order in approx_orders:
                for dpn in [1, mesh.dim]:
                    self.report('geometry: %s, approx. order: %d, dpn: %d' %
                                (geometry, approx_order, dpn))

                    field = Field.from_args('fu', nm.float64, dpn, omega,
                                            approx_order=approx_order)

                    cc = field.get_coor()
                    dofs = nm.zeros((field.n_nod, dpn), dtype=nm.float64)

                    for ic in range(dpn):
                        dofs[:, ic] = funs[ic](3 * (cc[:, 0] * cc[:, 1]))

                    vmesh, vdofs, level = field.linearize(dofs,
                                                          min_level=0,
                                                          max_level=3,
                                                          eps=1e-2)

                    if approx_order == 1:
                        _ok = level == 0

                    else:
                        _ok = level > 0
                    self.report('max. refinement level: %d: %s' % (level, _ok))

                    ok = ok and _ok

                    rdofs = nm.zeros((vmesh.n_nod, dpn), dtype=nm.float64)
                    cc = vmesh.coors
                    for ic in range(dpn):
                        rdofs[:, ic] = funs[ic](3 * (cc[:, 0] * cc[:, 1]))

                    _ok = nm.allclose(rdofs, vdofs, rtol=0.0, atol=0.03)
                    self.report('interpolation: %s' % _ok)
                    ok = ok and _ok

                    out = {
                        'u' : Struct(name='output_data',
                                     mode='vertex', data=vdofs,
                                     var_name='u', dofs=None)
                    }

                    name = self.join('linearizer-%s-%d-%d'
                                     % (geometry, approx_order, dpn))

                    vmesh.write(name + '.mesh')
                    vmesh.write(name + '.vtk', out=out)

        return ok