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
##
# 02.10.2007, c
[docs]class Evaluator( Struct ):
pass
##
# 02.10.2007, c
[docs]class BasicEvaluator( Evaluator ):
def __init__(self, problem, matrix_hook=None):
Evaluator.__init__(self, problem=problem,
matrix_hook=matrix_hook)
[docs] def new_ulf_iteration(self, nls, vec, it, err, err0):
pb = self.problem
vec = self.make_full_vec(vec)
pb.equations.set_variables_from_state(vec)
upd_vars = pb.conf.options.get('mesh_update_variables', None)
for varname in upd_vars:
try:
state = pb.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 = pb.domain.get_mesh_coors().copy()
vs = state()
coors[nods,:] = coors[nods,:] + vs.reshape(len(nods), state.n_components)
if pb.ts.step == 1 and it == 0:
state.field.save_mappings()
state.field.clear_mappings()
pb.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:
vec = self.make_full_vec( vec )
vec_r = self.problem.equations.eval_residuals(vec)
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:
vec = self.make_full_vec( vec )
pb = self.problem
if mtx is None:
mtx = pb.mtx_a
mtx = pb.equations.eval_tangent_matrices(vec, mtx)
if self.matrix_hook is not None:
mtx = self.matrix_hook(mtx, pb, call_mode='basic')
return mtx
[docs] def make_full_vec( self, vec ):
return self.problem.equations.make_full_vec(vec)
##
# 04.10.2007, c
[docs]class LCBCEvaluator( BasicEvaluator ):
##
# 04.10.2007, c
def __init__(self, problem, matrix_hook=None):
BasicEvaluator.__init__(self, problem, matrix_hook=matrix_hook)
self.mtx_lcbc = problem.equations.get_lcbc_operator()
##
# 04.10.2007, c
[docs] def eval_residual( self, vec, is_full = False ):
if not is_full:
vec = self.make_full_vec( vec )
vec_r = BasicEvaluator.eval_residual( self, vec, is_full = True )
vec_rr = self.mtx_lcbc.T * vec_r
return vec_rr
##
# 04.10.2007, c
[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:
vec = self.make_full_vec( vec )
mtx = BasicEvaluator.eval_tangent_matrix( self, vec, mtx = mtx,
is_full = True )
mtx_r = self.mtx_lcbc.T * mtx * self.mtx_lcbc
mtx_r = mtx_r.tocsr()
mtx_r.sort_indices()
## import pylab
## from sfepy.base.plotutils import spy
## spy( mtx_r )
## pylab.show()
## print mtx_r.__repr__()
if self.matrix_hook is not None:
mtx_r = self.matrix_hook(mtx_r, self.problem, call_mode='lcbc')
return mtx_r
[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,
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.
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[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, 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,
verbose=verbose)
else:
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,
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.
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_stripped_state_vector()
else:
asm_obj = equations.create_matrix_graph(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, ig, iels, coors,
fields, materials, variables,
functions=None, mode='eval', term_mode=None,
extra_args=None, 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.
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 = fields.values()[0].domain
region = Region('Elements', 'given elements', domain, '')
region.cells = iels + domain.mesh.el_offsets[ig]
region.update_shape()
domain.regions.append(region)
for field in fields.itervalues():
field.clear_mappings(clear_all=True)
for ap in field.aps.itervalues():
ap.clear_qp_base()
aux = create_evaluable(expression, fields, materials,
variables.itervalues(), Integrals([integral]),
functions=functions,
mode=mode, extra_args=extra_args, verbose=verbose,
kwargs=kwargs)
equations, variables = aux
out = eval_equations(equations, variables,
preserve_caches=False,
mode=mode, term_mode=term_mode)
domain.regions.pop()
return out
[docs]def assemble_by_blocks(conf_equations, problem, ebcs=None, epbcs=None,
dw_mode='matrix'):
"""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_iso.i1.Y2( inclusion.lame, 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 conf_equations.iteritems():
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], stripped = True, allow_dual = True )
ic = indx( var_names[1], stripped = True, allow_dual = True )
problem.update_materials()
mtx = problem.evaluate(mtx_term, auto_init=True,
mode='weak', dw_mode='matrix',
copy_materials=False)
matrices[mtx_name] = mtx[ir,ic]
return matrices