Fixed stupid errors.

apps/tests/monomial_bases_1D.cpp

@@ -80,7 +80,7 @@ int main(void)
     auto ge = 100.0*std::sqrt(grad_error/grad_norm);
 
     bool fepass = false, gepass = false;
-    if (fe < 1.65e-12) fepass = true;
+    if (fe < 1.91e-12) fepass = true;
     if (ge < 8.79e-5) gepass = true;
 
     auto passfail = [](bool pass) {

apps/tests/monomial_bases_2D.cpp

@@ -128,20 +128,22 @@ int test_basis_functions(const Mesh& msh)
     auto dofs_f = π(f);
 
     /* Error on the projected function */
-    auto fun_error = integrate(msh, cl, 2*degree+1, [&](const point_type& pt) {
-        auto s = f(pt) - eval(pt, dofs_f, basis);
+    auto errf = [&](const point_type& pt) {
+        T s = f(pt) - eval(pt, dofs_f, basis);
         return s*s;
-    }); 
+    };
+    auto fun_error = integrate(msh, cl, 2*degree+1, errf); 
 
     /* Norm of the projected function */
-    auto fun_norm = integrate(msh, cl, 2*degree+1, [&](const point_type& pt) {
-        auto s = f(pt);
+    auto normf = [&](const point_type& pt) {
+        T s = f(pt);
         return s*s;
-    });
+    };
+    auto fun_norm = integrate(msh, cl, 2*degree+1, normf);
 
     /* Error on the gradient of the projected function */
     auto grad_error = integrate(msh, cl, 2*degree+1, [&](const point_type& pt) {
-        auto s = grad_f(pt) - eval(pt, dofs_f, grad(basis));
+        Eigen::Matrix<T,2,1> s = grad_f(pt) - eval(pt, dofs_f, grad(basis));
         return s.dot(s);
     });
 
@@ -166,25 +168,38 @@ int test_basis_functions(const Mesh& msh)
     std::cout << "  Function relative error = " << fe << "% " << passfail(fepass) << std::endl;
     std::cout << "  Gradient relative error = " << ge << "% " << passfail(gepass) << std::endl;
 
+    auto fcs = faces(msh, cl);
+    for (auto& fc : fcs)
+    {
+        auto face_basis = scaled_monomial_basis(msh, fc, degree);
+        L2_projector πF(msh, fc, face_basis);
+        auto dofs_f_F = πF(f);
+        T fun_error_f = 0.0;
+        T fun_norm_f = 0.0;
+        auto fqps = integrate(msh, fc, 2*degree+2);
+        for (auto& qp : fqps) {
+            T diff = f(qp.point()) - eval(qp.point(), dofs_f_F, face_basis);
+            fun_error_f += qp.weight() * (diff*diff);
+            fun_norm_f += qp.weight() * normf(qp.point());
+        }
+        std::cout << "    Face Relerr: " << 100*std::sqrt(fun_error_f/fun_norm_f) << std::endl;
+    }
+
+
     return not (fepass and gepass);
 }
 
 int main(void)
 {
     using T = double;
 
-    point<T,2> a{0, 0};
-    point<T,2> b{M_PI, M_PI};
-    point<T,2> c{-M_PI, M_PI};
+    simplicial_mesh<T,2> msh_simp;
+    make_single_element_mesh(msh_simp, {0,0}, {1,0}, {0,1});
+    test_basis_functions(msh_simp);
 
-    simplicial_mesh<T,2> msh;
-    make_single_element_mesh(msh, a, b, c);
-    test_basis_functions(msh);
-
-    //point<T,2> base;
-    //cartesian_mesh<T,2> msh;
-    //make_single_element_mesh(msh, {0.5, 0.5}, 1.0, M_PI);
-    //test_basis_functions(msh);
+    cartesian_mesh<T,2> msh_cart;
+    make_single_element_mesh(msh_cart, {0.5, 0.5}, 1.0, M_PI);
+    test_basis_functions(msh_cart);
 
     return 0;
-}
+}

libdiskpp/include/diskpp/bases/bases_new.hpp

@@ -198,15 +198,15 @@ class scalar_monomial<Mesh<CoordT, 2, Storage>, typename Mesh<CoordT, 2, Storage
         auto pts = points(msh, cl);
         assert(pts.size() >= 3);
         v0_ = (pts[1] - pts[0]).to_vector();
-        v1_ = (pts[2] - pts[1]).to_vector();
+        v1_ = (pts[2] - pts[0]).to_vector();
         v1_ = v1_ - v1_.dot(v0_)*v0_/(v0_.dot(v0_)); // Gram-Schmidt
         tr_.col(0) = v0_;
         tr_.col(1) = v1_;
         tr_ = tr_.inverse().eval();
     }
 
     point_type phys2ref(const point_type& pt) const {
-        auto v = tr_*(pt - bar_).to_vector();
+        Eigen::Matrix<coordinate_type,2,1> v = tr_*(pt - bar_).to_vector();
         return point_type(v(0), v(1));
     }
 
@@ -321,11 +321,11 @@ class scalar_monomial<Mesh<CoordT, 2, Storage>, typename Mesh<CoordT, 2, Storage
     size_t      degree_;
     size_t      size_;
 
-    coordinate_type phys2ref(const point_type& pt) {
+    coordinate_type phys2ref(const point_type& pt) const {
         auto bv = (pb_ - bar_).to_vector();
         auto dv = (pt - bar_).to_vector();
-        auto nbv = h_/2.;
-        return rs_point_type( bv.dot(dv)/(nbv*nbv) );
+        auto nbv = h_;
+        return bv.dot(dv)/(nbv*nbv);
     }
 
 public:
@@ -370,7 +370,7 @@ class scalar_monomial<Mesh<CoordT, 2, Storage>, typename Mesh<CoordT, 2, Storage
     }
 
     size_t integration_degree() const {
-        return 0;
+        return degree_;
     }
 };
 
@@ -501,4 +501,4 @@ class L2_scalar_product
 };
 
 
-} // namespace disk::basis
+} // namespace disk::basis

libdiskpp/include/diskpp/mesh/meshgen.hpp

@@ -416,8 +416,9 @@ void make_single_element_mesh(cartesian_mesh<T,2>& msh, const point<T,2>& base,
 
     storage->points.push_back( point_type(    base.x(),    base.y() ) );
     storage->points.push_back( point_type( base.x()+hx,    base.y() ) );
-    storage->points.push_back( point_type( base.x()+hx, base.y()+hy ) );
     storage->points.push_back( point_type(    base.x(), base.y()+hy ) );
+    storage->points.push_back( point_type( base.x()+hx, base.y()+hy ) );
+
 
     auto p0 = point_identifier<2>(0);
     auto p1 = point_identifier<2>(1);

libdiskpp/include/diskpp/mesh/meshgen_fvca5hex.hpp

@@ -11,6 +11,7 @@
 #pragma once
 
 #include <iostream>
+#include <optional>
 #include <fstream>
 #include <vector>
 #include <map>
@@ -374,4 +375,4 @@ offset(const hexmesh& hm, const point& pt, bool one_based)
 }
 
 
-} //namespace disk::fvca5
+} //namespace disk::fvca5

libdiskpp/include/diskpp/quadratures/quad_hexahedral.hpp

@@ -39,19 +39,29 @@ std::vector<disk::quadrature_point<T, 2>>
 integrate(const disk::cartesian_mesh<T, 2>& msh, const typename disk::cartesian_mesh<T, 2>::cell& cl, const size_t degree)
 {
     const auto pts = points(msh, cl);
-    return disk::quadrature::tensorized_gauss_legendre(degree, pts);
+    std::array<point<T,2>, 4> pp;
+    pp[0] = pts[0];
+    pp[1] = pts[1];
+    pp[2] = pts[3]; // in cartesian_mesh<T,2> the nodes are
+    pp[3] = pts[2]; // *not* stored in ccw order (see issue #46)
+    return disk::quadrature::tensorized_gauss_legendre(degree, pp);
 }
 
 template<typename T>
 std::vector<disk::quadrature_point<T, 2>>
 integrate(const disk::cartesian_mesh<T, 2>& msh, const typename disk::cartesian_mesh<T, 2>::face& fc, const size_t degree)
 {
+    /*
     if (degree == 0)
     {
         return priv::integrate_degree0(msh, fc);
     }
 
     return priv::integrate_2D_face(msh, fc, degree);
+    */
+    const auto pts = points(msh, fc);
+    assert(pts.size() == 2);
+    return disk::quadrature::gauss_legendre(degree, pts[0], pts[1]);
 }
 
 namespace priv {

libdiskpp/include/diskpp/quadratures/quad_raw_gauss.hpp

@@ -161,11 +161,12 @@ gauss_legendre(size_t degree, const point<T,DIM>& a, const point<T,DIM>& b)
     if (rule_num >= num_rules)
         throw std::invalid_argument("gauss_legendre: order too high");
 
-    auto npts = priv::gauss_rules[rule_num].num_entries;
-    std::vector<quadrature_point<T,1>> qps;
+    const auto& qrule = priv::gauss_rules[rule_num];
+    auto npts = qrule.num_entries;
+    std::vector<quadrature_point<T,DIM>> qps;
     qps.reserve(npts);
     for (size_t i = 0; i < npts; i++) {
-        const auto &qp = priv::gauss_rules[rule_num].points[i];
+        const auto &qp = qrule.points[i];
         auto tr_qp = 0.5*(a+b) + 0.5*(b-a)*qp.point;
         auto tr_qw = distance(a,b)*qp.weight;
         qps.push_back({tr_qp, tr_qw});