Source code for sfepy.discrete.evaluate

from __future__ import absolute_import
from copy import copy

import numpy as nm

from sfepy.base.base import output, get_default, OneTypeList, Struct, basestr
from sfepy.discrete import Equations, Variables, Region, Integral, Integrals
from sfepy.discrete.common.fields import setup_extra_data
import six

[docs] def apply_ebc_to_matrix(mtx, ebc_rows, epbc_rows=None): """ Apply E(P)BC to matrix rows: put 1 to the diagonal for EBC DOFs, 1 to the diagonal for master EPBC DOFs, -1 to the [master, slave] entries. It is assumed, that the matrix contains zeros in EBC and master EPBC DOFs rows and columns. When used within a nonlinear solver, the actual values on the EBC DOFs diagonal positions do not matter, as the residual is zero at those positions. """ data, prows, cols = mtx.data, mtx.indptr, mtx.indices # Does not change the sparsity pattern. for ir in ebc_rows: for ic in range(prows[ir], prows[ir + 1]): if (cols[ic] == ir): data[ic] = 1.0 if epbc_rows is not None: import warnings from scipy.sparse import SparseEfficiencyWarning master, slave = epbc_rows with warnings.catch_warnings(): warnings.simplefilter('ignore', SparseEfficiencyWarning) # Changes sparsity pattern in-place - allocates new entries! The # master DOFs are not allocated by Equations.create_matrix_graph(), # see create_adof_conns(). mtx[master, master] = 1.0 mtx[master, slave] = -1.0
## # 02.10.2007, c
[docs] class Evaluator(Struct): """ This class provides the functions required by a nonlinear solver for a given problem. """ def __init__(self, problem, matrix_hook=None): Struct.__init__(self, problem=problem, matrix_hook=matrix_hook)
[docs] @staticmethod def new_ulf_iteration(problem, nls, vec, it, err, err0): vec = problem.equations.make_full_vec(vec) problem.equations.set_state(vec) upd_vars = problem.conf.options.get('mesh_update_variables', None) for varname in upd_vars: try: state = problem.equations.variables[varname] except IndexError: msg = 'variable "%s" does not exist!' % varname raise KeyError( msg ) nods = state.field.get_dofs_in_region(state.field.region, merge=True) coors = problem.domain.get_mesh_coors().copy() coors[nods, :] += state().reshape(len(nods), state.n_components) if len(state.field.mappings0) == 0: state.field.save_mappings() state.field.clear_mappings() problem.set_mesh_coors(coors, update_fields=False, actual=True, clear_all=False)
[docs] def eval_residual(self, vec, is_full=False): if not is_full and self.problem.active_only: vec = self.make_full_vec(vec) else: self.problem.equations.variables.apply_ebc(vec) vec_r = self.problem.equations.eval_residuals(vec) if self.matrix_hook is not None: vec_r = self.matrix_hook(vec_r, self.problem, call_mode='residual') if self.problem.equations.variables.has_lcbc: mtx_lcbc = self.problem.equations.get_lcbc_operator() vec_rr = mtx_lcbc.T * vec_r if self.matrix_hook is not None: vec_rr = self.matrix_hook(vec_rr, self.problem, call_mode='lcbc_residual') vec_r = vec_rr return vec_r
[docs] def eval_tangent_matrix(self, vec, mtx=None, is_full=False): if isinstance(vec, basestr) and vec == 'linear': return get_default(mtx, self.problem.mtx_a) if not is_full and self.problem.active_only: vec = self.make_full_vec(vec) elif vec is not None: self.problem.equations.variables.apply_ebc(vec) pb = self.problem if mtx is None: mtx = pb.mtx_a mtx = pb.equations.eval_tangent_matrices(vec, mtx) if not pb.active_only: apply_ebc_to_matrix(mtx, *pb.get_ebc_indices()) if self.matrix_hook is not None: mtx = self.matrix_hook(mtx, pb, call_mode='basic') if self.problem.equations.variables.has_lcbc: mtx_lcbc = self.problem.equations.get_lcbc_operator() mtx_r = mtx_lcbc.T * mtx * mtx_lcbc mtx_r = mtx_r.tocsr() mtx_r.sort_indices() if self.matrix_hook is not None: mtx_r = self.matrix_hook(mtx_r, self.problem, call_mode='lcbc') mtx = mtx_r return mtx
[docs] def make_full_vec(self, vec): return self.problem.equations.make_full_vec(vec)
[docs] def create_evaluable(expression, fields, materials, variables, integrals, regions=None, ebcs=None, epbcs=None, lcbcs=None, ts=None, functions=None, auto_init=False, mode='eval', extra_args=None, active_only=True, eterm_options=None, verbose=True, kwargs=None): """ Create evaluable object (equations and corresponding variables) from the `expression` string. Parameters ---------- expression : str The expression to evaluate. fields : dict The dictionary of fields used in `variables`. materials : Materials instance The materials used in the expression. variables : Variables instance The variables used in the expression. integrals : Integrals instance The integrals to be used. regions : Region instance or list of Region instances The region(s) to be used. If not given, the regions defined within the fields domain are used. ebcs : Conditions instance, optional The essential (Dirichlet) boundary conditions for 'weak' mode. epbcs : Conditions instance, optional The periodic boundary conditions for 'weak' mode. lcbcs : Conditions instance, optional The linear combination boundary conditions for 'weak' mode. ts : TimeStepper instance, optional The time stepper. functions : Functions instance, optional The user functions for boundary conditions, materials etc. auto_init : bool Set values of all variables to all zeros. mode : one of 'eval', 'el_avg', 'qp', 'weak' The evaluation mode - 'weak' means the finite element assembling, 'qp' requests the values in quadrature points, 'el_avg' element averages and 'eval' means integration over each term region. extra_args : dict, optional Extra arguments to be passed to terms in the expression. active_only : bool If True, in 'weak' mode, the (tangent) matrices and residual vectors (right-hand sides) contain only active DOFs. eterm_options : dict, optional The einsum-based terms evaluation options. verbose : bool If False, reduce verbosity. kwargs : dict, optional The variables (dictionary of (variable name) : (Variable instance)) to be used in the expression. Returns ------- equation : Equation instance The equation that is ready to be evaluated. variables : Variables instance The variables used in the equation. """ if kwargs is None: kwargs = {} if regions is not None: if isinstance(regions, Region): regions = [regions] regions = OneTypeList(Region, regions) else: regions = fields[list(fields.keys())[0]].domain.regions # Create temporary variables. aux_vars = Variables(variables) if extra_args is None: extra_args = kwargs else: extra_args = copy(extra_args) extra_args.update(kwargs) if ts is not None: extra_args.update({'ts' : ts}) equations = Equations.from_conf({'tmp' : expression}, aux_vars, regions, materials, integrals, user=extra_args, eterm_options=eterm_options, verbose=verbose) equations.collect_conn_info() # The true variables used in the expression. variables = equations.variables if auto_init: for var in variables: var.init_data(step=0) if mode == 'weak': equations.time_update(ts, ebcs, epbcs, lcbcs, functions, active_only=active_only, verbose=verbose) else: for eq in equations: for term in eq.terms: term.time_update(ts) setup_extra_data(equations.conn_info) return equations, variables
[docs] def eval_equations(equations, variables, names=None, preserve_caches=False, mode='eval', dw_mode='vector', term_mode=None, active_only=True, any_dof_conn=False, verbose=True): """ Evaluate the equations. Parameters ---------- equations : Equations instance The equations returned by :func:`create_evaluable()`. variables : Variables instance The variables returned by :func:`create_evaluable()`. names : str or sequence of str, optional Evaluate only equations of the given name(s). preserve_caches : bool If True, do not invalidate evaluate caches of variables. mode : one of 'eval', 'el_avg', 'qp', 'weak' The evaluation mode - 'weak' means the finite element assembling, 'qp' requests the values in quadrature points, 'el_avg' element averages and 'eval' means integration over each term region. dw_mode : 'vector' or 'matrix' The assembling mode for 'weak' evaluation mode. term_mode : str The term call mode - some terms support different call modes and depending on the call mode different values are returned. active_only : bool If True, in 'weak' mode, the (tangent) matrices and residual vectors (right-hand sides) contain only active DOFs. any_dof_conn : bool If True, in 'weak' `mode` and 'matrix' `dw_mode`, all DOF connectivities are used to pre-allocate the matrix graph. If False, only cell region connectivities are used. verbose : bool If False, reduce verbosity. Returns ------- out : dict or result The evaluation result. In 'weak' mode it is the vector or sparse matrix, depending on `dw_mode`. Otherwise, it is a dict of results with equation names as keys or a single result for a single equation. """ asm_obj = None if mode == 'weak': if dw_mode == 'vector': asm_obj = equations.create_reduced_vec() else: asm_obj = equations.create_matrix_graph(any_dof_conn=any_dof_conn, active_only=active_only, verbose=verbose) if not preserve_caches: equations.invalidate_term_caches() out = equations.evaluate(names=names, mode=mode, dw_mode=dw_mode, term_mode=term_mode, asm_obj=asm_obj) if variables.has_lcbc and mode == 'weak': mtx_lcbc = variables.mtx_lcbc if dw_mode == 'vector': out = mtx_lcbc.T * out elif dw_mode == 'matrix': out = mtx_lcbc.T * out * mtx_lcbc out = out.tocsr() out.sort_indices() return out
[docs] def eval_in_els_and_qp(expression, iels, coors, fields, materials, variables, functions=None, mode='eval', term_mode=None, extra_args=None, active_only=True, verbose=True, kwargs=None): """ Evaluate an expression in given elements and points. Parameters ---------- expression : str The expression to evaluate. fields : dict The dictionary of fields used in `variables`. materials : Materials instance The materials used in the expression. variables : Variables instance The variables used in the expression. functions : Functions instance, optional The user functions for materials etc. mode : one of 'eval', 'el_avg', 'qp' The evaluation mode - 'qp' requests the values in quadrature points, 'el_avg' element averages and 'eval' means integration over each term region. term_mode : str The term call mode - some terms support different call modes and depending on the call mode different values are returned. extra_args : dict, optional Extra arguments to be passed to terms in the expression. active_only : bool If True, in 'weak' mode, the (tangent) matrices and residual vectors (right-hand sides) contain only active DOFs. verbose : bool If False, reduce verbosity. kwargs : dict, optional The variables (dictionary of (variable name) : (Variable instance)) to be used in the expression. Returns ------- out : array The result of the evaluation. """ weights = nm.ones_like(coors[:, 0]) integral = Integral('ie', coors=coors, weights=weights) domain = list(fields.values())[0].domain region = Region('Elements', 'given elements', domain, '') region.cells = iels region.update_shape() domain.regions.append(region) for field in six.itervalues(fields): field.clear_mappings(clear_all=True) field.clear_qp_basis() aux = create_evaluable(expression, fields, materials, variables.itervalues(), Integrals([integral]), functions=functions, mode=mode, extra_args=extra_args, active_only=active_only, verbose=verbose, kwargs=kwargs) equations, variables = aux out = eval_equations(equations, variables, preserve_caches=False, mode=mode, term_mode=term_mode, active_only=active_only) domain.regions.pop() return out
[docs] def assemble_by_blocks(conf_equations, problem, ebcs=None, epbcs=None, dw_mode='matrix', active_only=True): """Instead of a global matrix, return its building blocks as defined in `conf_equations`. The name and row/column variables of each block have to be encoded in the equation's name, as in:: conf_equations = { 'A,v,u' : "dw_lin_elastic.i1.Y2( inclusion.D, v, u )", } Notes ----- `ebcs`, `epbcs` must be either lists of BC names, or BC configuration dictionaries. """ if isinstance( ebcs, list ) and isinstance( epbcs, list ): bc_mode = 0 elif isinstance( ebcs, dict ) and isinstance( epbcs, dict ): bc_mode = 1 else: raise TypeError('bad BC!') matrices = {} for key, mtx_term in six.iteritems(conf_equations): ks = key.split( ',' ) mtx_name, var_names = ks[0], ks[1:] output( mtx_name, var_names ) problem.set_equations({'eq': mtx_term}) variables = problem.get_variables() indx = variables.get_indx if bc_mode == 0: problem.select_bcs( ebc_names = ebcs, epbc_names = epbcs ) else: problem.time_update(ebcs=ebcs, epbcs=epbcs) ir = indx(var_names[0], reduced=True, allow_dual=True) ic = indx(var_names[1], reduced=True, allow_dual=True) problem.update_materials() mtx = problem.evaluate(mtx_term, auto_init=True, mode='weak', dw_mode='matrix', copy_materials=False, active_only=active_only) matrices[mtx_name] = mtx[ir,ic] return matrices