Commit cb05fcd5 authored by Andreas Dedner's avatar Andreas Dedner

added a test for the quadrature rules (both python and C++ algorithm)

parent fc030b9f
Pipeline #12739 passed with stage
in 13 minutes and 58 seconds
......@@ -12,6 +12,10 @@ def _duneIntegrate(self,entity,f):
_duneQuadratureRules = {}
def quadratureRule(geometryType, order):
try:
geometryType = geometryType.type
except AttributeError:
pass
try:
rule = _duneQuadratureRules[(geometryType,order)]
except KeyError:
......@@ -20,7 +24,7 @@ def quadratureRule(geometryType, order):
_duneQuadratureRules[(geometryType,order)] = rule
return rule
def quadratureRules(order):
return lambda entity: quadratureRule(entity.type,order)
return lambda entity: quadratureRule(entity,order)
def integrate(rules,entity,f):
return rules(entity).apply(entity,f)
......@@ -35,7 +35,7 @@ try:
self.weights_ = transform.get_vol(self.vertices_)*self.quadrature_.weights
self.quadPoints_ = [ QPQuadPoint(p,w) for p,w in zip(self.points_,self.weights_) ]
self.points_ = self.points_.transpose().copy()
def get():
def get(self):
return self.points_, self.weights_
def apply(self,entity,f):
ie = entity.geometry.integrationElement
......@@ -47,6 +47,10 @@ try:
def order(self):
return self.quadrature_.degree
def rule(gt,quadDescription):
try:
gt = gt.type
except AttributeError:
pass
try:
quad = _cache[(gt,quadDescription)]
return quad
......@@ -94,7 +98,12 @@ try:
return ret
def rules(methods):
return lambda entity: rule(entity.type, methods[entity.type])
def r(entity):
try:
return rule(entity.type, methods[entity.type])
except AttributeError:
return rule(entity, methods[entity])
return r
except ImportError as e:
logger.warning('Unable to import quadpy: ' + " ".join(str(e).splitlines()))
......
from __future__ import print_function, unicode_literals
import time, math, numpy
import dune
from dune.generator import algorithm
from dune.grid import cartesianDomain, yaspGrid
domain = cartesianDomain([0, 0], [1, 0.25], [60, 16])
yaspView = yaspGrid(domain)
@dune.grid.gridFunction(yaspView)
def function(x):
return numpy.cos(2.*numpy.pi/(0.3+abs(x[0]*x[1])))
def interpolate(grid):
mapper = grid.mapper({dune.geometry.vertex: 1})
data = numpy.zeros(mapper.size)
for v in grid.vertices:
data[mapper.index(v)] = function(v.geometry.center)
return mapper, data
mapper, data = interpolate(yaspView)
@dune.grid.gridFunction(yaspView)
def p12dEvaluate(element,x):
indices = mapper(element)
bary = 1-x[0]-x[1], x[0], x[1]
return sum( bary[i] * data[indices[i]] for i in range(3) )
@dune.grid.gridFunction(yaspView)
def error(element,x):
return numpy.abs(p12dEvaluate(element,x)-function(element,x))
rules = dune.geometry.quadratureRules(5)
start = time.time()
l2norm2 = 0
for e in yaspView.elements:
hatxs, hatws = rules(e.type).get()
weights = hatws * e.geometry.integrationElement(hatxs)
l2norm2 += numpy.sum(error(e, hatxs)**2 * weights, axis=-1)
print("Python:",math.sqrt(l2norm2))
print("time used:", round(time.time()-start,2))
algo = algorithm.load('l2norm2', 'test_quad.hh', yaspView, rules, error)
start = time.time()
l2norm2 = algo(yaspView,rules,error)
print("C++:",math.sqrt(l2norm2))
print("time used:", round(time.time()-start,2))
try:
import dune.geometry.quadpy as quadpy
rules = quadpy.rules({dune.geometry.quadrilateral: ("C2 7-2","Stroud")})
start = time.time()
l2norm2 = 0
for e in yaspView.elements:
hatxs, hatws = rules(e.type).get()
weights = hatws * e.geometry.integrationElement(hatxs)
l2norm2 += numpy.sum(error(e, hatxs)**2 * weights, axis=-1)
print("Python:",math.sqrt(l2norm2))
print("time used:", round(time.time()-start,2))
start = time.time()
l2norm2 = algo(yaspView,rules,error)
print("C++:",math.sqrt(l2norm2))
print("time used:", round(time.time()-start,2))
except ImportError:
pass
#include <cstddef>
#include <dune/common/fvector.hh>
// #include <dune/python/grid/entity.hh>
#include <dune/python/pybind11/pybind11.h>
#include <dune/python/pybind11/numpy.h>
template< class GridView, class Rules, class GF >
double l2norm2 ( const GridView &gridView, const Rules &rules, const GF& gf )
{
auto lf = localFunction( gf );
double l2norm2 = 0;
for( const auto &entity : elements( gridView ) )
{
const auto geo = entity.geometry();
typedef typename decltype(geo)::LocalCoordinate LocalCoordinate;
lf.bind( entity );
pybind11::object pyrule = rules( geo.type() );
pybind11::object pyPW = pyrule.attr("get")();
auto pointsWeights = pyPW.template cast<
std::pair<pybind11::array_t<double>,
pybind11::array_t<double>> >();
const auto &valuesArray = lf( pointsWeights.first ).template cast< pybind11::array_t< double > >();
// check shape here...
auto values = valuesArray.template unchecked< 1 >();
for( std::size_t i = 0, sz = pointsWeights.second.size(); i < sz; ++i )
{
LocalCoordinate hatx(0);
for (int c=0;c<2;++c)
hatx[c] = pointsWeights.first.at(c,i);
double weight = pointsWeights.second.at( i ) * geo.integrationElement( hatx );
l2norm2 += (values[ i ] * values[ i ]) * weight;
}
lf.unbind();
}
return l2norm2;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment