Source code for sfepy.discrete.integrals
"""
Classes for accessing quadrature points and weights for various reference
element geometries.
"""
import numpy as nm
from sfepy.base.base import OneTypeList, Container, Struct, basestr
from quadratures import QuadraturePoints
[docs]class Integrals(Container):
"""
Container for instances of :class:`Integral`.
"""
[docs] @staticmethod
def from_conf(conf):
objs = OneTypeList(Integral)
for desc in conf.itervalues():
if hasattr(desc, 'vals'):
aux = Integral(desc.name,
coors=desc.vals,
weights=desc.weights)
else:
aux = Integral(desc.name,
order=desc.order)
objs.append(aux)
obj = Integrals(objs)
return obj
[docs] def get(self, name):
"""
Return existing or new integral.
Parameters
----------
name : str
The name can either be a non-negative integer, a string
representation of a non-negative integer (the integral
order) or 'a' (automatic order) or a string beginning with
'i' (existing custom integral name).
"""
if name == 'a':
raise NotImplementedError
elif isinstance(name, basestr) and (name[0] == 'i'):
try:
obj = self[name]
except IndexError:
raise ValueError('integral %s is not defined!' % name)
else:
try:
order = int(name)
except:
raise ValueError('unsupported integral reference! (%s)' % name)
name = '__o%d' % order
if self.has_key(name):
obj = self[name]
else:
# Create new integral, and add it to self.
obj = Integral(name, order=order)
self.append(obj)
return obj
[docs]class Integral(Struct):
"""
Wrapper class around quadratures.
"""
def __init__(self, name, order=1, coors=None, weights=None,
bounds=None, tp_fix=1.0, weight_fix=1.0, symmetric=False):
self.name = name
self.qps = {}
if coors is None:
self.mode = 'builtin'
else:
self.mode = 'custom'
self.coors = coors
self.weights = weights
self.bounds = bounds
self.tp_fix = tp_fix
self.weight_fix = weight_fix
self.symmetric = symmetric
self.order = 0
if order in ('auto', 'custom', 'a', 'c'):
self.order = -1
else:
self.order = int(order)
[docs] def get_qp(self, geometry):
"""
Get quadrature point coordinates and corresponding weights for given
geometry. For built-in quadratures, the integration order is given by
`self.order`.
Parameters
----------
geometry : str
The geometry key describing the integration domain,
see the keys of `sfepy.discrete.quadratures.quadrature_tables`.
Returns
-------
coors : array
The coordinates of quadrature points.
weights: array
The quadrature weights.
"""
if geometry in self.qps:
qp = self.qps[geometry]
else:
if self.mode == 'builtin':
qp = QuadraturePoints.from_table(geometry, self.order)
else:
qp = QuadraturePoints(None,
coors=self.coors, weights=self.weights,
bounds=self.bounds, tp_fix=self.tp_fix,
weight_fix=self.weight_fix,
symmetric=self.symmetric)
self.qps[geometry] = qp
return qp.coors, qp.weights
[docs] def integrate(self, function, order=1, geometry='1_2'):
"""
Integrate numerically a given scalar function.
Parameters
----------
function : callable(coors)
The function of space coordinates to integrate.
order : int, optional
The integration order. For tensor product geometries, this is the
1D (line) order.
geometry : str
The geometry key describing the integration domain. Default
is `'1_2'`, i.e. a line integral in [0, 1]. For other values
see the keys of `sfepy.discrete.quadratures.quadrature_tables`.
Returns
-------
val : float
The value of the integral.
"""
qp = QuadraturePoints.from_table(geometry, order)
fvals = function(qp.coors)
val = nm.sum(fvals * qp.weights)
return val