"""
Classes for accessing quadrature points and weights for various reference
element geometries.
"""
from __future__ import absolute_import
import numpy as nm
from sfepy.base.base import OneTypeList, Container, Struct, basestr
from .quadratures import QuadraturePoints
import six
[docs]
class Integrals(Container):
"""
Container for instances of :class:`Integral`.
"""
[docs]
@staticmethod
def from_conf(conf):
objs = OneTypeList(Integral)
for desc in six.itervalues(conf):
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 name in self.names:
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