Python discrete.Variables类代码示例

    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 standalone_setup(self):
        from sfepy.discrete import create_adof_conns, Variables

        conn_info = {'aux' : self.get_conn_info()}
        adcs = create_adof_conns(conn_info, None)

        variables = Variables(self.get_variables())
        variables.set_adof_conns(adcs)

        materials = self.get_materials(join=True)

        for mat in materials:
            mat.time_update(None, [Struct(terms=[self])])

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 __init__(self, equations):
        Container.__init__(self, equations)

        self.variables = Variables(self.collect_variables())
        self.materials = Materials(self.collect_materials())

        self.domain = self.get_domain()

        self.active_bcs = set()

        self.collect_conn_info()

    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 test_pbc( self ):
        from sfepy.discrete import Variables, Conditions

        problem  = self.problem
        conf = self.conf

        ebcs = Conditions.from_conf(conf.ebcs, problem.domain.regions)
        epbcs = Conditions.from_conf(conf.epbcs, problem.domain.regions)

        variables = Variables.from_conf(conf.variables, problem.fields)
        variables.equation_mapping(ebcs, epbcs, None, problem.functions)
        state = variables.create_state_vector()
        variables.apply_ebc(state)
        return variables.has_ebc(state)

    def test_consistency_d_dw(self):
        from sfepy.discrete import Variables

        ok = True
        pb = self.problem
        for aux in test_terms:
            term_template, (prefix, par_name, d_vars, dw_vars) = aux
            print term_template, prefix, par_name, d_vars, dw_vars

            term1 = term_template % ((prefix,) + d_vars)

            variables = Variables.from_conf(self.conf.variables, pb.fields)

            for var_name in d_vars:
                var = variables[var_name]
                n_dof = var.field.n_nod * var.field.shape[0]
                aux = nm.arange(n_dof, dtype=nm.float64)
                var.set_data(aux)

            if prefix == 'd':
                val1 = pb.evaluate(term1, var_dict=variables.as_dict())

            else:
                val1 = pb.evaluate(term1, call_mode='d_eval',
                                   var_dict=variables.as_dict())

            self.report('%s: %s' % (term1, val1))

            term2 = term_template % (('dw',) + dw_vars[:2])

            vec, vv = pb.evaluate(term2, mode='weak',
                                  var_dict=variables.as_dict(),
                                  ret_variables=True)

            pvec = vv.get_state_part_view(vec, dw_vars[2])
            val2 = nm.dot(variables[par_name](), pvec)
            self.report('%s: %s' % (term2, val2))

            err = nm.abs(val1 - val2) / nm.abs(val1)
            _ok = err < 1e-12
            self.report('relative difference: %e -> %s' % (err, _ok))

            ok = ok and _ok

        return ok

    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

class Equations(Container):

    @staticmethod
    def from_conf(conf, variables, regions, materials, integrals,
                  user=None, verbose=True):

        objs = OneTypeList(Equation)

        conf = copy(conf)

        ii = 0
        for name, desc in conf.iteritems():
            if verbose:
                output('equation "%s":' %  name)
                output(desc)
            eq = Equation.from_desc(name, desc, variables, regions,
                                    materials, integrals, user=user)
            objs.append(eq)
            ii += 1

        obj = Equations(objs)

        return obj

    def __init__(self, equations):
        Container.__init__(self, equations)

        self.variables = Variables(self.collect_variables())
        self.materials = Materials(self.collect_materials())

        self.domain = self.get_domain()

        self.active_bcs = set()

        self.collect_conn_info()

    def create_subequations(self, var_names, known_var_names=None):
        """
        Create sub-equations containing only terms with the given virtual
        variables.

        Parameters
        ----------
        var_names : list
            The list of names of virtual variables.
        known_var_names : list
            The list of  names of (already) known state variables.

        Returns
        -------
        subequations : Equations instance
            The sub-equations.
        """
        from sfepy.discrete import FieldVariable

        known_var_names = get_default(known_var_names, [])

        objs = []
        for iv, var_name in enumerate(var_names):
            terms = [term.copy(name=term.name)
                     for eq in self for term in eq.terms
                     if term.get_virtual_name() == var_name]

            # Make parameter variables from known state variables in terms
            # arguments.
            for known_name in known_var_names:
                for term in terms:
                    if known_name in term.arg_names:
                        ii = term.arg_names.index(known_name)
                        state = self.variables[known_name]
                        par = FieldVariable(known_name, 'parameter',
                                            state.field,
                                            primary_var_name='(set-to-None)')
                        term.args[ii] = par
                        term._kwargs[known_name] = par
                        par.set_data(state())

            new_terms = Terms(terms)
            objs.append(Equation('eq_%d' % iv, new_terms))

        subequations = Equations(objs)

        return subequations

    def get_domain(self):
        domain = None

        for eq in self:
            for term in eq.terms:
                if term.has_region:
                    domain = term.region.domain

        return domain

    def collect_materials(self):
        """
        Collect materials present in the terms of all equations.
        """
        materials = []
        for eq in self:
#.........这里部分代码省略.........