Commit aa554992 authored by Andreas Dedner's avatar Andreas Dedner

fix an issue with the lagrange interpolation on quads

parent 9c6d3846
Pipeline #12857 passed with stage
in 7 minutes and 2 seconds
......@@ -4,12 +4,12 @@ import dune
from dune.generator import algorithm
from dune.grid import cartesianDomain, yaspGrid
domain = cartesianDomain([0, 0], [1, 0.25], [60, 16])
domain = cartesianDomain([0, 0], [1, 0.25], [12, 3])
yaspView = yaspGrid(domain)
@dune.grid.gridFunction(yaspView)
def function(x):
return numpy.cos(2.*numpy.pi/(0.3+abs(x[0]*x[1])))
return numpy.cos(2.*numpy.pi/(0.3+x[0]*x[1]))
def interpolate(grid):
mapper = grid.mapper({dune.geometry.vertex: 1})
data = numpy.zeros(mapper.size)
......@@ -20,46 +20,50 @@ 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) )
bary = (1-x[0])*(1-x[1]), x[0]*(1-x[1]), (1-x[0])*x[1], x[0]*x[1]
return sum( bary[i] * data[indices[i]] for i in range(4) )
@dune.grid.gridFunction(yaspView)
def error(element,x):
return numpy.abs(p12dEvaluate(element,x)-function(element,x))
return 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")})
if True:
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))
print("Python:",math.sqrt(l2norm2),flush=True)
print("time used:", round(time.time()-start,2),flush=True)
if True:
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))
print("C++:",math.sqrt(l2norm2),flush=True)
print("time used:", round(time.time()-start,2),flush=True)
try:
import dune.geometry.quadpy as quadpy
rules = quadpy.rules({dune.geometry.quadrilateral: ("C2 7-2","Stroud")})
if True:
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),flush=True)
print("time used:", round(time.time()-start,2),flush=True)
if True:
start = time.time()
l2norm2 = algo(yaspView,rules,error)
print("C++:",math.sqrt(l2norm2),flush=True)
print("time used:", round(time.time()-start,2),flush=True)
except ImportError:
pass
......@@ -31,6 +31,11 @@ double l2norm2 ( const GridView &gridView, const Rules &rules, const GF& gf )
for (std::size_t c=0;c<LocalCoordinate::size();++c)
hatx[c] = pointsWeights.first.at(c,i);
double weight = pointsWeights.second.at( i ) * geo.integrationElement( hatx );
#if 0
std::cout << hatx << " " << pointsWeights.second.at( i )
<< " " << weight << " " << geo.global(hatx)
<< " " << values[i] << std::endl;
#endif
l2norm2 += (values[ i ] * values[ i ]) * weight;
}
lf.unbind();
......
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