diff --git a/hysop/constants.py b/hysop/constants.py
index 7ccba265279bc7af6811ae875b33c589ba2c39b6..ca123a9eb14a81c0fc0bbf038ce658c169c938f6 100644
--- a/hysop/constants.py
+++ b/hysop/constants.py
@@ -147,6 +147,27 @@ BoundaryCondition = EnumFactory.create('BoundaryCondition',
'NEUMANN', 'DIRICHLET' ])
"""Boundary conditions enum"""
+def boundary2str(b):
+ """Helper function to convert a BoundaryCondition to a short string."""
+ sstr = {
+ BoundaryCondition.NONE: 'NONE',
+ BoundaryCondition.MIXED: 'MIXED',
+ BoundaryCondition.PERIODIC: 'PER',
+ BoundaryCondition.HOMOGENEOUS_NEUMANN: 'HNEU',
+ BoundaryCondition.HOMOGENEOUS_DIRICHLET: 'HDIR',
+ BoundaryCondition.NEUMANN: 'NEU',
+ BoundaryCondition.DIRICHLET: 'DIR',
+ }
+ if b in sstr:
+ return sstr[b]
+ else:
+ return str(b)
+
+def format_boundaries(lboundaries, rboundaries):
+ """Helper function format boundary conditions."""
+ return ','.join(('{}/{}'.format(boundary2str(lb), boundary2str(rb))
+ for (lb,rb) in zip(lboundaries, rboundaries)))
+
DomainExtension = EnumFactory.create('DomainExtension',
['PERIODIC', 'EVEN', 'ODD'])
diff --git a/hysop/constants.py.in b/hysop/constants.py.in
index 1149c4c46fd35e0da863e052473e272f2d109baa..35d1b4745174086d0fc3170849287f5e70e78700 100644
--- a/hysop/constants.py.in
+++ b/hysop/constants.py.in
@@ -146,6 +146,27 @@ BoundaryCondition = EnumFactory.create('BoundaryCondition',
'HOMOGENEOUS_NEUMANN', 'HOMOGENEOUS_DIRICHLET',
'NEUMANN', 'DIRICHLET' ])
"""Boundary conditions enum"""
+
+def boundary2str(b):
+ """Helper function to convert a BoundaryCondition to a short string."""
+ sstr = {
+ BoundaryCondition.NONE: 'NONE',
+ BoundaryCondition.MIXED: 'MIXED',
+ BoundaryCondition.PERIODIC: 'PER',
+ BoundaryCondition.HOMOGENEOUS_NEUMANN: 'HNEU',
+ BoundaryCondition.HOMOGENEOUS_DIRICHLET: 'HDIR',
+ BoundaryCondition.NEUMANN: 'NEU',
+ BoundaryCondition.DIRICHLET: 'DIR',
+ }
+ if b in sstr:
+ return sstr[b]
+ else:
+ return str(b)
+
+def format_boundaries(lboundaries, rboundaries):
+ """Helper function format boundary conditions."""
+ return ','.join(('{}/{}'.format(boundary2str(lb), boundary2str(rb))
+ for (lb,rb) in zip(lboundaries, rboundaries)))
DomainExtension = EnumFactory.create('DomainExtension',
['PERIODIC', 'EVEN', 'ODD'])
diff --git a/hysop/fields/continuous_field.py b/hysop/fields/continuous_field.py
index e5e0c82e1d60f5b79bd44bb941a6027716127599..7ba53d130a0d917d8df4089d99293af3b5ff8414 100644
--- a/hysop/fields/continuous_field.py
+++ b/hysop/fields/continuous_field.py
@@ -244,7 +244,9 @@ class FieldContainerI(TaggedObject):
pass
@classmethod
- def from_sympy_expressions(cls, name, exprs, space_symbols, pretty_name=None, dtype=None, **kwds):
+ def from_sympy_expressions(cls, name, exprs, space_symbols,
+ scalar_name_prefix=None, scalar_pretty_name_prefix=None,
+ pretty_name=None, **kwds):
"""
Create a field wich has the same shape as exprs, with optional names.
Expressions should be of kind sympy.Expr and are converted to FieldExpression: this
@@ -256,33 +258,66 @@ class FieldContainerI(TaggedObject):
if isinstance(exprs, sm.Expr):
raise NotImplementedError('Call self.from_sympy_expression instead.')
check_instance(exprs, npw.ndarray, values=sm.Expr)
+ check_instance(name, str)
+ check_instance(pretty_name, (str, unicode), allow_none=True)
+ check_instance(scalar_name_prefix, str, allow_none=True)
+ check_instance(scalar_pretty_name_prefix, (str, unicode), allow_none=True)
+ if isinstance(pretty_name, unicode):
+ pretty_name = pretty_name.encode('utf-8')
+ if isinstance(scalar_pretty_name_prefix, unicode):
+ scalar_pretty_name_prefix = scalar_pretty_name_prefix.encode('utf-8')
+
fields = npw.empty(shape=exprs.shape, dtype=object)
fields[...] = None
- for idx in npw.nditer(*exprs.shape):
- _name = name
- _pretty_name = pretty_name
+ for idx in npw.ndindex(*exprs.shape):
+ if (scalar_name_prefix is not None):
+ sname = TensorField.default_name_formatter(scalar_name_prefix, idx)
+ if (scalar_pretty_name_prefix is not None):
+ spname = TensorField.default_pretty_name_formatter(scalar_pretty_name_prefix, idx)
+ else:
+ spname = TensorField.default_pretty_name_formatter(scalar_name_prefix, idx)
+ else:
+ # names will be autogenerated from sympy expression
+ sname = None
+ spname = None
+
fields[idx] = cls.from_sympy_expression(expr=exprs[idx], space_symbols=space_symbols,
- name=_name, pretty_name=_pname, dtype=dtype, **kwds)
+ name=sname, pretty_name=spname, **kwds)
return TensorField.from_field_array(name=name, pretty_name=pretty_name,
fields=fields)
@classmethod
def from_sympy_expression(cls, expr, space_symbols, **kwds):
from hysop.symbolic.field import FieldExpressionBuilder
+ from hysop.tools.field_utils import print_all_names
assert 'lboundaries' not in kwds
assert 'rboundaries' not in kwds
assert 'domain' not in kwds
+
+ # determine names if not given
+ if ('name' not in kwds) or (kwds['name'] is None):
+ (name, pretty_name, var_name, latex_name) = print_all_names(expr)
+ kwds['name'] = name
+ kwds['pretty_name'] = pretty_name
+ kwds['var_name'] = var_name
+ kwds['latex_name'] = latex_name
+
+ # determine domain and boundary conditions
fe = FieldExpressionBuilder.to_field_expression(
expr=expr, space_symbols=space_symbols, strict=True)
kwds['domain'] = fe.domain
kwds['lboundaries'] = fe.lboundaries
kwds['rboundaries'] = fe.rboundaries
- kwds['name'] = first_not_None(kwds.get('name', None), 'lol')
- kwds['pretty_name'] = first_not_None(kwds.get('pretty_name', None), 'lol')
+
+ # deduce data type from field expression if not specified
+ kwds['dtype'] = first_not_None(kwds.get('dtype', None), fe.dtype)
+
+ # finally return create and return the ScalarField
return ScalarField(**kwds)
def gradient(self, name=None, pretty_name=None,
+ scalar_name_prefix=None, scalar_pretty_name_prefix=None,
directions=None, axis=-1,
space_symbols=None,
dtype=None, **kwds):
@@ -290,11 +325,12 @@ class FieldContainerI(TaggedObject):
Create a field capable of storing the gradient of self,
possibly altered.
"""
- dim = self.dim
- ndim = self.ndim
+ dim = self.dim # dimension of the domain
+ ndim = self.ndim # number of dimension of the np.ndarray
+ frame = self.domain.frame
directions = to_tuple(first_not_None(directions, range(dim)))
- space_symbols = to_tuple(first_not_None(space_symbols, self.domain.frame.coords))
+ space_symbols = to_tuple(first_not_None(space_symbols, frame.coords))
check_instance(directions, tuple, minval=0, maxval=self.dim-1, minsize=1, unique=True)
check_instance(axis, int, minval=-ndim, maxval=ndim-1)
check_instance(space_symbols, tuple, values=SpaceSymbol, size=dim, unique=True)
@@ -306,50 +342,115 @@ class FieldContainerI(TaggedObject):
else:
shape = (ndirs,)
+ name = first_not_None(name, 'grad_{}'.format(self.name))
+ pretty_name = first_not_None(pretty_name, '{}{}'.format(nabla.encode('utf8'),
+ self.pretty_name))
+
if shape==(1,):
- expr = sm.Derivative(self.symbol, space_symbols[directions[0]])
+ expr = self.symbol(frame.time, *space_symbols).diff(space_symbols[directions[0]])
return self.from_sympy_expression(expr=expr, space_symbols=space_symbols,
name=name, pretty_name=pretty_name,
dtype=dtype, **kwds)
else:
exprs = npw.empty(shape=shape, dtype=object)
- for idx in npw.nditer(*shape):
+ for idx in npw.ndindex(*shape):
i = idx[:axis+1] + idx[axis+2:]
d = directions[idx[axis+1]]
- exprs[idx] = sm.Derivative(self[i].symbol, space_symbols[d])
+ exprs[idx] = self[i].symbol(frame.time, *space_symbols).diff(space_symbols[d])
return self.from_sympy_expressions(exprs=exprs, space_symbols=space_symbols,
name=name, pretty_name=pretty_name,
+ scalar_name_prefix=scalar_name_prefix,
+ scalar_pretty_name_prefix=scalar_pretty_name_prefix,
dtype=dtype, **kwds)
-
- #if shape==(1,):
- #name = first_not_None(name, '{d}{}_dx{}'.format(
- #F.name, directions[0], d='d'))
- #pname = first_not_None(pretty_name, u'{d}{}/{d}{x}{}'.format(
- #F.pretty_name.decode('utf-8'),
- #subscript(directions[0]), d=partial, x=xsymbol))
- #return self.field_like(name=name, pretty_name=pretty_name)
- #else:
- #name = first_not_None(name, 'grad_{}'.format(self.name))
- #pretty_name = first_not_None(pretty_name, '{}{}'.format(nabla.encode('utf8'),
- #self.pretty_name))
- #def make_field(idx, **fkwds):
- #i = idx[:axis+1] + idx[axis+2:]
- #d = directions[idx[axis+1]]
-
- #if (ndim==0):
- #Fi = self
- #else:
- #Fi = self[i]
- #fkwds['dtype'] = first_not_None(dtype, Fi.dtype)
- #fkwds['name'] = '{d}{}_{d}x{}'.format(Fi.name, d, d='d')
- #fkwds['pretty_name'] = u'{d}{}/{d}{x}{}'.format(Fi.pretty_name.decode('utf-8'),
- #subscript(d), d=partial, x=xsymbol)
+ def laplacian(self, name=None, pretty_name=None,
+ scalar_name_prefix=None, scalar_pretty_name_prefix=None,
+ dtype=None, **kwds):
+ from hysop.symbolic.field import laplacian
+ frame = self.domain.frame
+ exprs = laplacian(self.symbol(*frame.vars), frame)
+
+ name = first_not_None(name, 'laplacian_{}'.format(self.name))
+ pretty_name = first_not_None(pretty_name, u'\u0394{}'.format(
+ self.pretty_name.decode('utf-8')))
+
+ if (exprs.size == 1):
+ expr = npw.asscalar(exprs)
+ return self.from_sympy_expression(expr=expr, space_symbols=frame.coords,
+ name=name, pretty_name=pretty_name,
+ dtype=dtype, **kwds)
+ else:
+ return self.from_sympy_expressions(exprs=exprs, space_symbols=frame.coords,
+ name=name, pretty_name=pretty_name,
+ scalar_name_prefix=scalar_name_prefix,
+ scalar_pretty_name_prefix=scalar_pretty_name_prefix,
+ dtype=dtype, **kwds)
- #return Fi.field_like(**fkwds)
- #return TensorField(name=name, pretty_name=pretty_name, domain=domain, shape=shape,
- #make_field=make_field, **kwds)
+ def div(self, name=None, pretty_name=None,
+ scalar_name_prefix=None, scalar_pretty_name_prefix=None,
+ axis=-1, dtype=None, **kwds):
+ """
+ Create a field capable of storing the divergence of self,
+ on chosen axis.
+ """
+ from hysop.symbolic.field import div
+ frame = self.domain.frame
+ exprs = npw.asarray(div(self.symbol(*frame.vars), frame))
+
+ name = first_not_None(name, 'div_{}'.format(self.name))
+ pretty_name = first_not_None(pretty_name, u'{}\u22c5{}'.format(nabla,
+ self.pretty_name.decode('utf-8')))
+
+ if exprs.size in (0,1):
+ expr = npw.asscalar(exprs)
+ return self.from_sympy_expression(expr=expr, space_symbols=frame.coords,
+ name=name, pretty_name=pretty_name,
+ dtype=dtype, **kwds)
+ else:
+ return self.from_sympy_expressions(exprs=exprs, space_symbols=frame.coords,
+ name=name, pretty_name=pretty_name,
+ scalar_name_prefix=scalar_name_prefix,
+ scalar_pretty_name_prefix=scalar_pretty_name_prefix,
+ dtype=dtype, **kwds)
+
+ def curl(self, name=None, pretty_name=None,
+ scalar_name_prefix=None, scalar_pretty_name_prefix=None,
+ dtype=None, **kwds):
+ """
+ Create a field capable of storing the curl of self,
+
+ Only 2D and 3D fields are supported as the curl brings
+ a 1-vector to a 2-vector (ie. a vector to a pseudoscalar
+ in 2D and vector to a pseudovector in 3D).
+
+ In 1D the curl is 0, and in 4D the curl would be a 6D 'field'.
+ """
+ from hysop.symbolic.field import curl
+ msg='Can only take curl for a 2D or 3D vector field.'
+ assert (self.dim in (2,3)) and (self.nb_components == self.dim), msg
+ frame = self.domain.frame
+ exprs = curl(self.symbol(*frame.vars), frame)
+
+ name = first_not_None(name, 'curl_{}'.format(self.name))
+ pretty_name = first_not_None(pretty_name, u'{}\u2227{}'.format(nabla,
+ self.pretty_name.decode('utf-8')))
+
+ if (self.dim==2):
+ return self.from_sympy_expression(expr=exprs, space_symbols=frame.coords,
+ name=name, pretty_name=pretty_name,
+ dtype=dtype, **kwds)
+ else:
+ return self.from_sympy_expressions(exprs=exprs, space_symbols=frame.coords,
+ name=name, pretty_name=pretty_name,
+ scalar_name_prefix=scalar_name_prefix,
+ scalar_pretty_name_prefix=scalar_pretty_name_prefix,
+ dtype=dtype, **kwds)
+
+ def rot(self, *args, **kwds):
+ """See curl."""
+ return self.curl(*args, **kwds)
+
def get_attributes(self, *attrs):
"""
@@ -459,6 +560,7 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
@debug
def __new__(cls, domain, name, pretty_name=None,
+ var_name=None, latex_name=None,
initial_values=0, dtype=HYSOP_REAL,
lboundaries=None, rboundaries=None,
is_tmp=False, **kwds):
@@ -474,6 +576,12 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
pretty_name: string or unicode, optional.
A pretty name used for display whenever possible.
Defaults to name.
+ var_name: string, optional.
+ A variable name used for code generation.
+ This will be passed to the symbolic representation of this field.
+ latex_name: string, optional.
+ A variable name used for latex generation.
+ This will be passed to the symbolic representation of this field.
dtype: npw.dtype, optional, defaults to HYSOP_REAL
Underlying data type of this field
initial_values: numeric value, or tuple of numeric values, optional
@@ -510,8 +618,11 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
periodicity: numpy.ndarray of bool
Numpy array mask, True is axis is periodic, else False.
"""
- check_instance(is_tmp, bool)
check_instance(name, str)
+ check_instance(pretty_name, (str, unicode), allow_none=True)
+ check_instance(latex_name, str, allow_none=True)
+ check_instance(var_name, str, allow_none=True)
+ check_instance(is_tmp, bool)
# Data type of the field
if (dtype==npw.bool) or (dtype==bool):
@@ -569,7 +680,9 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
# Symbolic representation of this field
from hysop.symbolic.field import SymbolicField
- obj._symbol = SymbolicField(field=obj)
+ obj._symbol = SymbolicField(field=obj,
+ var_name=var_name,
+ latex_name=latex_name)
# Dictionnary of all the discretizations of this field.
# keys are hysop.topology.topology.Topology,
@@ -649,31 +762,32 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
return s
def format_boundaries(self):
- return ','.join(('{}/{}'.format(
- str(lb).replace('HOMOGENEOUS_','')[:3],
- str(rb).replace('HOMOGENEOUS_','')[:3])
- for (lb,rb) in zip(*self.boundaries)))
+ from hysop.constants import format_boundaries as fb
+ return fb(*self.boundaries)
def long_description(self):
"""Long description of this field."""
s = textwrap.dedent(
- '''{}
+ '''
+ {}
*name: {}
- *pname: {}
+ *pretty_name: {}
+ *var_name: {}
+ *latex_name: {}
*dim: {}
*dtype: {}
- *symbolic repr.: {}
*left boundary: {}
*right boundary: {}
*initial values: {}
*topology tags: [{}]
''').format(self.full_tag,
- self.name, self.pretty_name, self.dim,
- self.dtype, self.symbol,
- self.lboundaries.tolist(), self.rboundaries.tolist(),
+ self.name, self.pretty_name,
+ self.symbol._var_name, self.symbol._latex_name,
+ self.dim, self.dtype,
+ self.lboundaries.tolist(), self.rboundaries.tolist(),
self.initial_values,
','.join([k.full_tag for k in self.discrete_fields.keys()]))
- return s
+ return s[1:]
@debug
@@ -1004,7 +1118,7 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
'''
{}
*name: {}
- *pname: {}
+ *pretty_name: {}
*dim: {}
*shape: {}
*nb_components: {}
diff --git a/hysop/fields/default_fields.py b/hysop/fields/default_fields.py
index a9ae41167e84ad4e64234892bc9f61ae9fa8da33..462f00248b05fe55d319dc255483b57cecbaf537 100644
--- a/hysop/fields/default_fields.py
+++ b/hysop/fields/default_fields.py
@@ -53,7 +53,9 @@ def VorticityField(velocity, name=None, pretty_name=None, **kwds):
assert velocity.nb_components == domain.dim, 'Invalid velocity Field.'
name = first_not_None(name, 'W')
pretty_name = first_not_None(pretty_name, greak[24])
- return velocity.curl(name=name, pretty_name=pretty_name, **kwds)
+ return velocity.curl(name=name, pretty_name=pretty_name,
+ scalar_name_prefix=name, scalar_pretty_name_prefix=pretty_name,
+ **kwds)
def DensityField(domain, name=None, pretty_name=None, **kwds):
diff --git a/hysop/fields/tests/test_fields.py b/hysop/fields/tests/test_fields.py
index b14a279eaede44f90a6c9f0466031960d4adf817..0284189b9f7ccb3e6b214ef1d06f6646d9d2e055 100644
--- a/hysop/fields/tests/test_fields.py
+++ b/hysop/fields/tests/test_fields.py
@@ -339,26 +339,38 @@ def domain_boundary_iterator(dim):
def test_boundaries():
- for dim in (1,2):
+ """This test checks that all boundaries are compatible for velocity and vorticity."""
+ for dim in (1,2,): #3,):
i=0
for (lbd, rbd) in domain_boundary_iterator(dim):
domain = Box(dim=dim, lboundaries=lbd,
rboundaries=rbd)
V = VelocityField(domain)
- print V.gradient()
- W = VorticityField(V)
S = ScalarField(name='S0', domain=domain)
+ divV = V.div()
+ gradV = V.gradient()
+ lapV = V.laplacian()
print
print 'DOMAIN BOUNDARIES:'
print ' *boundaries=[{}]'.format(domain.format_boundaries())
+ print 'SCALAR BOUNDARIES:'
+ print ' *{} boundaries=[{}]'.format(S.pretty_name, S.format_boundaries())
print 'VELOCITY BOUNDARIES:'
for Vi in V.fields:
- print ' *{} boundaries=[{}]'.format(Vi.name, Vi.format_boundaries())
- print 'VORTICITY BOUNDARIES:'
- for Wi in W.fields:
- print ' *{} boundaries=[{}]'.format(Wi.name, Wi.format_boundaries())
- print 'SCALAR BOUNDARIES:'
- print ' *{} boundaries=[{}]'.format(S.name, S.format_boundaries())
+ print ' *{} boundaries=[{}]'.format(Vi.pretty_name, Vi.format_boundaries())
+ print '{} BOUNDARIES:'.format(divV.pretty_name)
+ print ' *{} boundaries=[{}]'.format(divV.pretty_name, divV.format_boundaries())
+ print '{} BOUNDARIES:'.format(gradV.pretty_name)
+ for gVi in gradV.fields:
+ print ' *{} boundaries=[{}]'.format(gVi.pretty_name, gVi.format_boundaries())
+ print '{} BOUNDARIES:'.format(lapV.pretty_name)
+ for lVi in lapV.fields:
+ print ' *{} boundaries=[{}]'.format(lVi.pretty_name, lVi.format_boundaries())
+ if (dim>1):
+ rotV = V.curl()
+ print '{} (VORTICITY) BOUNDARIES:'.format(rotV.pretty_name)
+ for Wi in rotV.fields:
+ print ' *{} boundaries=[{}]'.format(Wi.pretty_name, Wi.format_boundaries())
diff --git a/hysop/symbolic/field.py b/hysop/symbolic/field.py
index d4d2215eccf447961802d985e6877f418decb731..5475b1ee36d80bdc391dbb1b8bd66a41a81c52d3 100644
--- a/hysop/symbolic/field.py
+++ b/hysop/symbolic/field.py
@@ -1,8 +1,11 @@
from abc import abstractmethod
from hysop.deps import sm
+from hysop.constants import BoundaryCondition
+
from hysop.tools.numpywrappers import npw
from hysop.tools.types import check_instance, first_not_None
from hysop.tools.sympy_utils import get_derivative_variables
+from hysop.tools.numerics import find_common_dtype
from hysop.fields.continuous_field import Field, TensorField
from hysop.fields.discrete_field import DiscreteField, DiscreteTensorField
@@ -10,7 +13,6 @@ from hysop.symbolic import Symbol
from hysop.symbolic.base import TensorBase, SymbolicTensor
from hysop.symbolic.func import UndefinedFunction, AppliedSymbolicFunction, FunctionBase, \
SymbolicFunctionTensor
-from hysop.constants import BoundaryCondition
from hysop.domain.domain import Domain
class FieldExpressionI(object):
@@ -26,20 +28,23 @@ class FieldExpressionI(object):
def domain(self):
pass
+ @abstractmethod
+ def dtype(self):
+ pass
+
@property
def boundaries(self):
return (self.lboundaries, self.rboundaries)
def format_boundaries(self):
- return ','.join(('{}/{}'.format(
- str(lb).replace('HOMOGENEOUS_','')[:3],
- str(rb).replace('HOMOGENEOUS_','')[:3])
- for (lb,rb) in zip(*self.boundaries)))
+ from hysop.constants import format_boundaries as fb
+ return fb(*self.boundaries)
class FieldExpression(FieldExpressionI):
def __init__(self, *args, **kwds):
self._domain = kwds.pop('domain', None)
+ self._dtype = kwds.pop('dtype', None)
self._lboundaries = kwds.pop('lboundaries', None)
self._rboundaries = kwds.pop('rboundaries', None)
super(FieldExpression, self).__init__(*args, **kwds)
@@ -74,6 +79,16 @@ class FieldExpression(FieldExpressionI):
check_instance(dom, Domain)
self._domain = dom
+ @property
+ def dtype(self):
+ assert (self._dtype is not None)
+ return self._dtype
+ @dtype.setter
+ def dtype(self, dt):
+ assert (self._dtype is None)
+ check_instance(dt, npw.dtype)
+ self._dtype = dt
+
class FieldExpressionBuilder(object):
class BoundaryIncompatibilityError(ValueError):
@@ -112,7 +127,8 @@ class FieldExpressionBuilder(object):
pass
e = _to_field_expression_impl(expr.args[0])
if cls.is_field_expr(e):
- lb, rb, domain = e.lboundaries.copy(), e.rboundaries.copy(), e.domain
+ dtype, domain = e.dtype, e.domain
+ lb, rb = e.lboundaries.copy(), e.rboundaries.copy(),
assert len(space_symbols)==lb.size==rb.size
for xi in get_derivative_variables(expr):
assert xi in space_symbols, xi
@@ -121,6 +137,7 @@ class FieldExpressionBuilder(object):
rb[i] = cls.update_boundaries(rb[i], +1)
expr = DerivativeFieldExpr(e, *expr.args[1:])
expr.domain = domain
+ expr.dtype = dtype
expr.lboundaries = lb
expr.rboundaries = rb
return expr
@@ -129,10 +146,10 @@ class FieldExpressionBuilder(object):
else:
func = expr.func
args = tuple(_to_field_expression_impl(a) for a in expr.args)
- field_expression_args = filter(lambda x: cls.is_field_expr(x), args)
+ field_expression_args = tuple(filter(lambda x: cls.is_field_expr(x), args))
if field_expression_args:
try:
- return FieldExpression.make_expr(func, args)
+ return cls.make_expr(func, *args)
except cls.BoundaryIncompatibilityError:
msg='\nError during the handling of expression {}.'.format(expr)
msg+='\nSome boundaries were not compatible:'
@@ -144,32 +161,34 @@ class FieldExpressionBuilder(object):
fexpr = _to_field_expression_impl(expr)
if strict and (not cls.is_field_expr(fexpr)):
msg='\nError during the handling of expression {}.'.format(expr)
- msg+='Could not determine boundaries because no FieldExpression '
+ msg+='\nCould not determine boundaries because no FieldExpression '
msg+='was present in expression.'
raise cls.InvalidExpression(msg)
return fexpr
@classmethod
- def mk_expr(cls, func, *args):
+ def make_expr(cls, func, *args):
check_instance(func, type)
- field_expression_args = filter(lambda x: cls.is_field_expr(x), args)
- if field_expression_args:
- msg='No FieldExpression arguments.'
+ field_expression_args = tuple(filter(lambda x: cls.is_field_expr(x), args))
+ if not field_expression_args:
+ msg='No FieldExpression arguments present in args.'
raise ValueError(msg)
if not cls.check_boundary_compatibility(*field_expression_args):
raise cls.BoundaryIncompatibilityError
+ fea0 = field_expression_args[0]
new_func = type(func.__name__+'FieldExpr', (FieldExpression, func), {})
new_expr = new_func(*args)
- new_expr.domain = args[0].domain
- new_expr.lboundaries = args[0].lboundaries.copy()
- new_expr.rboundaries = args[0].rboundaries.copy()
+ new_expr.dtype = npw.dtype(find_common_dtype(*tuple(a.dtype for a in field_expression_args)))
+ new_expr.domain = fea0.domain
+ new_expr.lboundaries = fea0.lboundaries.copy()
+ new_expr.rboundaries = fea0.rboundaries.copy()
return new_expr
@classmethod
- def check_boundary_compatibility(arg0, *args):
+ def check_boundary_compatibility(cls, arg0, *args):
check_instance(args, tuple, values=FieldExpressionI)
- lb, rb = arg0.lboundaries, arg0.rboundaries
+ domain, lb, rb = arg0.domain, arg0.lboundaries, arg0.rboundaries
if args:
match = all((domain == a.domain) for a in args)
match &= all(all(lb==a.lboundaries) for a in args)
@@ -180,7 +199,7 @@ class FieldExpressionBuilder(object):
-class FieldBase(FieldExpressionI, FunctionBase):
+class FieldBase(FunctionBase):
def __new__(cls, field, idx=None, name=None, pretty_name=None, **kwds):
check_instance(field, (Field, DiscreteField))
assert (field.nb_components == 1) or (idx is not None), (field.nb_components, idx)
@@ -199,18 +218,6 @@ class FieldBase(FieldExpressionI, FunctionBase):
hc += (self._field, self._index,)
return hc
- @property
- def lboundaries(self):
- return self._field.lboundaries
-
- @property
- def rboundaries(self):
- return self._field.rboundaries
-
- @property
- def domain(self):
- return self._field.domain
-
@property
def field(self):
"""Get associated field."""
@@ -277,7 +284,7 @@ class SymbolicField(FieldBase, UndefinedFunction):
return not (self==other)
-class AppliedSymbolicField(AppliedSymbolicFunction):
+class AppliedSymbolicField(FieldExpressionI, AppliedSymbolicFunction):
"""Applied scalar fields, hold a reference to a continuous field."""
def __new__(cls, *args, **kwds):
args = args if args else cls.field.domain.frame.vars
@@ -302,6 +309,23 @@ class AppliedSymbolicField(AppliedSymbolicFunction):
def indexed_field(self):
"""Get a unique identifier for an indexed field component."""
return (self.field, self.index)
+
+ @property
+ def lboundaries(self):
+ return self.field.lboundaries
+
+ @property
+ def rboundaries(self):
+ return self.field.rboundaries
+
+ @property
+ def domain(self):
+ return self.field.domain
+
+ @property
+ def dtype(self):
+ return self.field.dtype
+
class SymbolicFieldTensor(SymbolicFunctionTensor):
"""Symbolic tensor symbol."""
@@ -358,17 +382,20 @@ def grad(F, frame, axis=-1):
return gradF.view(TensorBase)
def div(F, frame, axis=-1):
- assert isinstance(F, npw.ndarray)
- assert F.shape[axis] == frame.dim
- shape = F.shape
- ndim = F.ndim
- axis = (axis+ndim)%ndim
-
- divF = npw.empty_like(F)
- for idx in npw.ndindex(*shape):
- divF[idx] = diff(F[idx], frame.coords[idx[axis]])
- divF = divF.sum(axis=axis)
- return divF.view(TensorBase)
+ if isinstance(F, npw.ndarray):
+ assert F.shape[axis] == frame.dim
+ shape = F.shape
+ ndim = F.ndim
+ axis = (axis+ndim)%ndim
+
+ divF = npw.empty_like(F)
+ for idx in npw.ndindex(*shape):
+ divF[idx] = diff(F[idx], frame.coords[idx[axis]])
+ divF = divF.sum(axis=axis)
+ return divF.view(TensorBase)
+ else:
+ assert frame.dim==1
+ return F.diff(frame.coords[0])
def rot(F, frame):
assert (frame.dim in (2,3))
diff --git a/hysop/symbolic/frame.py b/hysop/symbolic/frame.py
index a8743ce3889e78f663d9b665f2a17062d061f255..888e9303d2e541fe360960a92ac8c78009b03d43 100644
--- a/hysop/symbolic/frame.py
+++ b/hysop/symbolic/frame.py
@@ -31,6 +31,7 @@ class SymbolicFrame(object):
"""Get the time variable for conveniance."""
return time_symbol
+
@property
def dtime(self):
"""Get the infinitesimal time variable for conveniance."""
diff --git a/hysop/tools/field_utils.py b/hysop/tools/field_utils.py
index fe09283d7f5fa74f3768708149bca7057451507a..72004242c7d8e4ab5fb05ce49168583c75b206a2 100644
--- a/hysop/tools/field_utils.py
+++ b/hysop/tools/field_utils.py
@@ -1,6 +1,146 @@
from hysop.tools.types import first_not_None, to_tuple
from hysop.tools.sympy_utils import nabla, partial, subscript, subscripts, \
- exponent, exponents, xsymbol
+ exponent, exponents, xsymbol, get_derivative_variables
+
+from sympy.printing.str import StrPrinter, StrReprPrinter
+from sympy.printing.ccode import C99CodePrinter
+from sympy.printing.latex import LatexPrinter
+
+class BasePrinter(object):
+ def print_Derivative(self, expr):
+ (bvar, pvar, vvar, lvar) = print_all_names(expr.args[0])
+ pvar = pvar.decode('utf-8')
+ all_xvars = get_derivative_variables(expr)
+ xvars = tuple(set(all_xvars))
+ varpows = tuple(all_xvars.count(x) for x in xvars)
+ bxvars = tuple(print_name(x) for x in xvars)
+ pxvars = tuple(print_pretty_name(x).decode('utf-8') for x in xvars)
+ vxvars = tuple(print_var_name(x) for x in xvars)
+ lxvars = tuple(print_latex_name(x) for x in xvars)
+ return DifferentialStringFormatter.format_pd(bvar, pvar, vvar, lvar,
+ bxvars, pxvars, vxvars, lxvars,
+ varpows=varpows)
+
+ def _print(self, expr, **kwds):
+ try:
+ return super(BasePrinter, self)._print(expr, **kwds)
+ except:
+ print
+ msg='FATAL ERROR: {} failed to print expression {}.'
+ msg=msg.format(type(self).__name__, expr)
+ print msg
+ print
+ raise
+
+
+class NamePrinter(BasePrinter, StrReprPrinter):
+ def _print(self, expr, **kwds):
+ if hasattr(expr, 'name'):
+ return expr.name
+ elif hasattr(expr, '_name'):
+ return expr._name
+ return super(NamePrinter, self)._print(expr, **kwds)
+ def _print_Derivative(self, expr):
+ return super(NamePrinter, self).print_Derivative(expr)[0]
+ def _print_Add(self, expr):
+ return super(NamePrinter, self)._print_Add(expr).replace(' ', '')
+ def _print_Mul(self, expr):
+ return super(NamePrinter, self)._print_Mul(expr).replace(' ', '')
+ def emptyPrinter(self, expr):
+ msg='\n{} does not implement _print_{}(self, expr).'
+ msg+='\nExpression is {}.'.format(expr)
+ msg+='\nExpression type MRO is:'
+ msg+='\n *'+'\n *'.join(t.__name__ for t in type(expr).__mro__)
+ msg=msg.format(self.__class__.__name__, expr.__class__.__name__)
+ raise NotImplementedError(msg)
+
+
+class PrettyNamePrinter(BasePrinter, StrPrinter):
+ def _print(self, expr, **kwds):
+ if hasattr(expr, 'pretty_name'):
+ return expr.pretty_name
+ elif hasattr(expr, '_pretty_name'):
+ return expr._pretty_name
+ return super(PrettyNamePrinter, self)._print(expr, **kwds)
+ def _print_Derivative(self, expr):
+ return super(PrettyNamePrinter, self).print_Derivative(expr)[1]
+ def emptyPrinter(self, expr):
+ msg='\n{} does not implement _print_{}(self, expr).'
+ msg+='\nExpression is {}.'.format(expr)
+ msg+='\nExpression type MRO is:'
+ msg+='\n *'+'\n *'.join(t.__name__ for t in type(expr).__mro__)
+ msg=msg.format(self.__class__.__name__, expr.__class__.__name__)
+ raise NotImplementedError(msg)
+
+
+class VarNamePrinter(BasePrinter, C99CodePrinter):
+ def _print(self, expr, **kwds):
+ if hasattr(expr, 'var_name'):
+ return expr.var_name
+ elif hasattr(expr, '_var_name'):
+ return expr._var_name
+ return super(VarNamePrinter, self)._print(expr, **kwds).replace(' ', '')
+ def _print_Derivative(self, expr):
+ return super(VarNamePrinter, self).print_Derivative(expr)[2]
+ def _print_Add(self, expr):
+ s = super(VarNamePrinter, self)._print_Add(expr)
+ s = s.replace(' + ', '_plus_').replace(' - ', '_minus_')
+ s = s.replace('+', 'plus_').replace('-', 'minus_')
+ return s
+ def _print_Mul(self, expr):
+ s = super(VarNamePrinter, self)._print_Mul(expr)
+ s = s.replace(' * ', '_times_')
+ return s
+ def emptyPrinter(self, expr):
+ msg='\n{} does not implement _print_{}(self, expr).'
+ msg+='\nExpression is {}.'.format(expr)
+ msg+='\nExpression type MRO is:'
+ msg+='\n *'+'\n *'.join(t.__name__ for t in type(expr).__mro__)
+ msg=msg.format(self.__class__.__name__, expr.__class__.__name__)
+ raise NotImplementedError(msg)
+
+
+class LatexNamePrinter(BasePrinter, LatexPrinter):
+ def _print(self, expr, **kwds):
+ if hasattr(expr, 'latex_name'):
+ return expr.latex_name
+ elif hasattr(expr, '_latex_name'):
+ return expr._latex_name
+ return super(LatexNamePrinter, self)._print(expr, **kwds)
+ def _print_Derivative(self, expr):
+ return super(LatexNamePrinter, self).print_Derivative(expr)[3]
+ def emptyPrinter(self, expr):
+ msg='\n{} does not implement _print_{}(self, expr).'
+ msg+='\nExpression is {}.'.format(expr)
+ msg+='\nExpression type MRO is:'
+ msg+='\n *'+'\n *'.join(t.__name__ for t in type(expr).__mro__)
+ msg=msg.format(self.__class__.__name__, expr.__class__.__name__)
+ raise NotImplementedError(msg)
+
+pbn = NamePrinter()
+ppn = PrettyNamePrinter()
+#pvn = VarNamePrinter()
+pln = LatexNamePrinter()
+
+def print_name(expr):
+ return pbn.doprint(expr)
+
+def print_pretty_name(expr):
+ return ppn.doprint(expr)
+
+def print_var_name(expr):
+ return VarNamePrinter().doprint(expr)
+
+def print_latex_name(expr):
+ return pln.doprint(expr)
+
+def print_all_names(expr):
+ name = print_name(expr)
+ pretty_name = print_pretty_name(expr)
+ var_name = print_var_name(expr)
+ latex_name = print_latex_name(expr)
+ return (name, pretty_name, var_name, latex_name)
+
def to_str(*args):
if len(args)==1:
@@ -181,6 +321,7 @@ class DifferentialStringFormatter(object):
varpows=1, var_components=None, xvars_components=None,
bdivide_fn=bdivide_fn, pdivide_fn=pdivide_fn, vdivide_fn=vdivide_fn, ldivide_fn=ldivide_fn,
fsc=True, **kwds):
+
for k in ('dpow', 'components', 'bvars', 'pvars', 'vvars', 'lvars', 'varpow'):
assert k not in kwds, 'Cannot specify reserved keyword {}.'.format(k)
@@ -278,5 +419,5 @@ if __name__ == '__main__':
_printDifferentialStringFormatter.format_pd(bvar, pvar, vvar, lvar,
bxvars, pxvars, vxvars, lxvars,
varpows=varpows, xvars_components=xvars_components,
- var_components=var_components), multiline=True)
+ var_components=var_components, multiline=True)
diff --git a/hysop/tools/sympy_utils.py b/hysop/tools/sympy_utils.py
index 91c59e3e124cfac7d60825306e6f0a17d08f23dd..d54ede72fb4754042ff79ebbe46a096aa8eac5a9 100644
--- a/hysop/tools/sympy_utils.py
+++ b/hysop/tools/sympy_utils.py
@@ -237,7 +237,6 @@ def tensor_xreplace(tensor,vars):
T[idx] = vars[symbol.name]
return T
-
def non_eval_xreplace(expr, rule):
"""
Duplicate of sympy's xreplace but with non-evaluate statement included.