Commit 889c02bb authored by Lukas Riedel's avatar Lukas Riedel

implemented easier quadrature rule syntax in error estimator

parent 2b8bcc39
#include<dune/geometry/referenceelements.hh>
#include<dune/pdelab/common/referenceelements.hh>
#include<dune/pdelab/common/quadraturerules.hh>
#include<dune/pdelab/localoperator/pattern.hh>
#include<dune/pdelab/localoperator/flags.hh>
#include<dune/pdelab/localoperator/idefault.hh>
#include<dune/pdelab/localoperator/defaultimp.hh>
#include<dune/pdelab/finiteelement/localbasiscache.hh>
namespace Dune {
......@@ -85,23 +92,22 @@ namespace Dune {
// geometric factor of penalty
const RF h_F = std::min(ig.inside().geometry().volume(),ig.outside().geometry().volume())/ig.geometry().volume(); // Houston!
// select quadrature rule
Dune::GeometryType gtface = ig.geometry().type();
const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
// get element geometry
auto gtface = ig.geometry();
// container for sum of errors
RF sum(0.0);
int counter = 0;
// loop over quadrature points and integrate normal flux
for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it = rule.begin(); it != rule.end(); ++it)
const auto quadrature_rule = quadratureRule(gtface,intorder);
for (const auto& it : quadrature_rule)
{
const Dune::FieldVector<DF,dim> n_F = ig.unitOuterNormal(it->position());
// retrieve unit normal vector
const Dune::FieldVector<DF,dim> n_F = ig.unitOuterNormal(it.position());
// position of quadrature point in local coordinates of elements
const Domain local_s = ig.geometryInInside().global(it->position());
const Domain local_n = ig.geometryInOutside().global(it->position());
const Domain local_s = ig.geometryInInside().global(it.position());
const Domain local_n = ig.geometryInOutside().global(it.position());
// evaluate basis functions
std::vector<Scalar> phi_s(lfsu_s.size());
......@@ -189,14 +195,14 @@ namespace Dune {
const RF jump = omega_s * condFactor_s * numFlux_s - omega_n * condFactor_n * numFlux_n;
// integration factor
const RF factor = it->weight() * ig.geometry().integrationElement(it->position());
const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
// add square of error
sum += jump*jump*factor*factor;
counter++;
}
// correct for the number of neighbors
sum *= 1. / counter;
sum *= 1. / quadrature_rule.size();
// what is this variable for?
DF h_T = 1;//std::max( diameter(ig.inside().geometry()),
......
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