Skip to content
Snippets Groups Projects
Commit a10f0cd7 authored by Jö Fahlke's avatar Jö Fahlke
Browse files

Eliminate std::vector from local operators 

This uses our own Q1 local basis implementation that hands out std::array
results
parent cb0c1534
No related branches found
No related tags found
1 merge request!18Don't use std::vectors inside local operators
Pipeline #19676 failed
Stage: test
    Hide whitespace changes
    Inline Side-by-side
    src/gridoperator.hh
    + 12
    − 16
    View file @ a10f0cd7
    ......@@ -16,8 +16,7 @@
    #include <dune/geometry/referenceelements.hh>
    #include <dune/localfunctions/lagrange/q1.hh>
    #include "q1localbasis.hh"
    #include "sparsitypattern.hh"
    namespace PPS {
    ......@@ -35,11 +34,11 @@ namespace PPS {
    class GridOperator
    {
    using DF = typename GridView::ctype;
    using LFE = Dune::Q1LocalFiniteElement<DF, RF, GridView::dimension>;
    using LB = Q1LocalBasis<DF, RF, GridView::dimension>;
    GridView gv_;
    LocalOperator lop_;
    LFE lfe_;
    LB lb_;
    std::set<std::size_t> cv_;
    public:
    ......@@ -67,7 +66,6 @@ namespace PPS {
    std::vector<typename Domain::field_type> xl;
    const auto &iset = gv_.indexSet();
    const auto &lb = lfe_.localBasis();
    // Traverse grid view
    for (const auto& element : elements(gv_))
    ......@@ -87,7 +85,7 @@ namespace PPS {
    Dune::Hybrid::ifElse
    (std::integral_constant<bool, LocalOperator::doAlphaVolume>{},
    [&](auto id) {
    id(lop_).jacobian_volume(element, lb, xl, al);
    id(lop_).jacobian_volume(element, lb_, xl, al);
    });
    // Skip if no intersection iterator is needed
    ......@@ -99,7 +97,7 @@ namespace PPS {
    if(intersection.boundary())
    {
    // Boundary integration
    id(lop_).jacobian_boundary(intersection, lb, xl, al);
    id(lop_).jacobian_boundary(intersection, lb_, xl, al);
    }
    });
    ......@@ -144,7 +142,6 @@ namespace PPS {
    std::vector<typename Range::field_type> yl;
    const auto &iset = gv_.indexSet();
    const auto &lb = lfe_.localBasis();
    // Traverse grid view
    for (const auto& element : elements(gv_))
    ......@@ -166,7 +163,7 @@ namespace PPS {
    }
    // Volume integration
    lop_.jacobian_apply_volume(element, lb, lpl, xl, yl);
    lop_.jacobian_apply_volume(element, lb_, lpl, xl, yl);
    Dune::Hybrid::ifElse
    (std::integral_constant<bool, LocalOperator::doAlphaBoundary>{},
    ......@@ -177,7 +174,7 @@ namespace PPS {
    {
    // Boundary integration
    id(lop_).
    jacobian_apply_boundary(intersection, lb, lpl, xl, yl);
    jacobian_apply_boundary(intersection, lb_, lpl, xl, yl);
    }
    });
    ......@@ -201,7 +198,6 @@ namespace PPS {
    std::vector<typename Range::field_type> rl;
    const auto &iset = gv_.indexSet();
    const auto &lb = lfe_.localBasis();
    // Traverse grid view
    for (const auto& element : elements(gv_))
    ......@@ -212,7 +208,7 @@ namespace PPS {
    rl.assign(size,0.0);
    // Volume integration
    lop_.lambda_volume(element, lb, rl);
    lop_.lambda_volume(element, lb_, rl);
    // Notify assembler engine about bind
    xl.resize(size);
    ......@@ -222,20 +218,20 @@ namespace PPS {
    xl[i] = x[iset.subIndex(element, i, element.dimension)];
    // Volume integration
    lop_.alpha_volume(element, lb, xl, rl);
    lop_.alpha_volume(element, lb_, xl, rl);
    // Traverse intersections
    for(const auto& intersection : intersections(gv_,element))
    if(intersection.boundary())
    {
    // Boundary integration
    lop_.lambda_boundary(intersection, lb, rl);
    lop_.lambda_boundary(intersection, lb_, rl);
    // Boundary integration
    Dune::Hybrid::ifElse
    (std::integral_constant<bool, LocalOperator::doAlphaBoundary>{},
    [&](auto id) {
    id(lop_).alpha_boundary(intersection, lb, xl, rl);
    id(lop_).alpha_boundary(intersection, lb_, xl, rl);
    });
    }
    ......@@ -263,7 +259,7 @@ namespace PPS {
    for (const auto& element : elements(gv_))
    {
    // Volume integration
    lop_.pattern_volume(lfe_, localpattern);
    lop_.pattern_volume(lb_, localpattern);
    // Notify assembler engine about unbinds
    for (auto link : localpattern)
    ......
    src/loppattern.hh
    + 3
    − 3
    View file @ a10f0cd7
    ......@@ -13,11 +13,11 @@ namespace PPS {
    public:
    // define sparsity pattern of operator representation
    template<class LocalFiniteElement, typename LocalPattern>
    void pattern_volume (const LocalFiniteElement& lfe,
    template<class LocalBasis, typename LocalPattern>
    void pattern_volume (const LocalBasis& lb,
    LocalPattern& pattern) const
    {
    auto size = lfe.size();
    auto size = lb.size();
    pattern.reserve(pattern.size() + size*size);
    for(auto row : Dune::range(size))
    for(auto col : Dune::range(size))
    ......
    src/nitschenonlinearpoissonfem.hh
    + 9
    − 13
    View file @ a10f0cd7
    #ifndef PACXX_PROJECTSEMINAR_2019_SRC_NITSCHENONLINEARPOISSONFEM_HH
    #define PACXX_PROJECTSEMINAR_2019_SRC_NITSCHENONLINEARPOISSONFEM_HH
    #include <array>
    #include <numeric>
    #include <vector>
    #include <dune/common/fvector.hh>
    #include <dune/common/rangeutilities.hh>
    ......@@ -62,7 +62,7 @@ namespace PPS {
    // types & dimension
    const int dim = Intersection::Entity::dimension;
    using RF = typename LocalBasis::Traits::RangeFieldType;
    using RF = typename LocalBasis::RangeField;
    // select quadrature rule
    auto geo = isect.geometry();
    ......@@ -71,9 +71,7 @@ namespace PPS {
    // geometry of inside cell
    auto geo_inside = isect.inside().geometry();
    std::vector<typename LocalBasis::Traits::RangeType> phihat;
    std::vector<typename LocalBasis::Traits::JacobianType> gradphi(lb.size());
    std::vector<typename LocalBasis::Traits::JacobianType> gradphihat;
    std::array<typename LocalBasis::Jacobian, lb.size()> gradphi;
    // loop over quadrature points
    for (const auto& ip : PPS::quadratureRule(geo,order))
    ......@@ -82,13 +80,13 @@ namespace PPS {
    auto local = localgeo.global(ip.position());
    // evaluate basis functions
    lb.evaluateFunction(local, phihat);
    auto phihat = lb.evaluateFunction(local);
    // evaluate u
    auto u = std::inner_product(begin(x), end(x), begin(phihat), RF(0));
    // evaluate gradient of shape functions
    lb.evaluateJacobian(local, gradphihat);
    auto gradphihat = lb.evaluateJacobian(local);
    // transform gradients of shape functions to real element
    const auto S = geo_inside.jacobianInverseTransposed(local);
    ......@@ -134,7 +132,7 @@ namespace PPS {
    // types & dimension
    const int dim = Intersection::Entity::dimension;
    using RF = typename LocalBasis::Traits::RangeFieldType;
    using RF = typename LocalBasis::RangeField;
    // select quadrature rule
    auto geo = isect.geometry();
    ......@@ -143,9 +141,7 @@ namespace PPS {
    // geometry of inside cell
    auto geo_inside = isect.inside().geometry();
    std::vector<typename LocalBasis::Traits::RangeType> phihat;
    std::vector<typename LocalBasis::Traits::JacobianType> gradphi(lb.size());
    std::vector<typename LocalBasis::Traits::JacobianType> gradphihat;
    std::array<typename LocalBasis::Jacobian, lb.size()> gradphi;
    // loop over quadrature points
    for (const auto& ip : PPS::quadratureRule(geo,order))
    ......@@ -154,13 +150,13 @@ namespace PPS {
    auto local = localgeo.global(ip.position());
    // evaluate basis functions
    lb.evaluateFunction(local, phihat);
    auto phihat = lb.evaluateFunction(local);
    // evaluate u
    auto u = std::inner_product(begin(x), end(x), begin(phihat), RF(0));
    // evaluate gradient of shape functions
    lb.evaluateJacobian(local, gradphihat);
    auto gradphihat = lb.evaluateJacobian(local);
    // transform gradients of shape functions to real element
    const auto S = geo_inside.jacobianInverseTransposed(local);
    ......
    src/nonlinearpoissonfem.hh
    + 11
    − 19
    View file @ a10f0cd7
    #ifndef PACXX_PROJECTSEMINAR_2019_SRC_NONLINEARPOISSONFEM_HH
    #define PACXX_PROJECTSEMINAR_2019_SRC_NONLINEARPOISSONFEM_HH
    #include <array>
    #include <numeric>
    #include <vector>
    #include <dune/common/fvector.hh>
    #include <dune/common/rangeutilities.hh>
    ......@@ -60,13 +60,11 @@ namespace PPS {
    auto geo = elem.geometry();
    const int order = incrementorder() + 2*lb.order();
    std::vector<typename LocalBasis::Traits::RangeType> phihat;
    // loop over quadrature points
    for (const auto& ip : PPS::quadratureRule(geo,order))
    {
    // evaluate basis functions
    lb.evaluateFunction(ip.position(), phihat);
    auto phihat = lb.evaluateFunction(ip.position());
    // integrate -f*phi_i
    decltype(ip.weight()) factor = ip.weight()*
    ......@@ -95,8 +93,6 @@ namespace PPS {
    auto globalgeo = isect.geometry();
    const int order = incrementorder() + 2*lb.order();
    std::vector<typename LocalBasis::Traits::RangeType> phihat;
    // loop over quadrature points and integrate normal flux
    for (const auto& ip : PPS::quadratureRule(globalgeo, order))
    {
    ......@@ -104,7 +100,7 @@ namespace PPS {
    auto local = localgeo.global(ip.position());
    // evaluate shape functions (assume Galerkin method)
    lb.evaluateFunction(local, phihat);
    auto phihat = lb.evaluateFunction(local);
    // integrate j
    decltype(ip.weight()) factor = ip.weight()*
    ......@@ -122,27 +118,25 @@ namespace PPS {
    {
    // types & dimension
    const int dim = Element::dimension;
    using RF = typename LocalBasis::Traits::RangeFieldType;
    using RF = typename LocalBasis::RangeField;
    // select quadrature rule
    auto geo = elem.geometry();
    const int order = incrementorder() + 2*lb.order();
    std::vector<typename LocalBasis::Traits::RangeType> phihat;
    std::vector<typename LocalBasis::Traits::JacobianType> gradphi(lb.size());
    std::vector<typename LocalBasis::Traits::JacobianType> gradphihat;
    std::array<typename LocalBasis::Jacobian, lb.size()> gradphi;
    // loop over quadrature points
    for (const auto& ip : PPS::quadratureRule(geo, order))
    {
    // evaluate basis functions
    lb.evaluateFunction(ip.position(), phihat);
    auto phihat = lb.evaluateFunction(ip.position());
    // evaluate u
    auto u = std::inner_product(begin(x), end(x), begin(phihat), RF(0));
    // evaluate gradient of shape functions
    lb.evaluateJacobian(ip.position(), gradphihat);
    auto gradphihat = lb.evaluateJacobian(ip.position());
    // transform gradients of shape functions to real element
    const auto S = geo.jacobianInverseTransposed(ip.position());
    ......@@ -176,21 +170,19 @@ namespace PPS {
    {
    // types & dimension
    const int dim = Element::dimension;
    using RF = typename LocalBasis::Traits::RangeFieldType;
    using RF = typename LocalBasis::RangeField;
    // select quadrature rule
    auto geo = elem.geometry();
    const int order = incrementorder() + 2*lb.order();
    std::vector<typename LocalBasis::Traits::RangeType> phihat;
    std::vector<typename LocalBasis::Traits::JacobianType> gradphi(lb.size());
    std::vector<typename LocalBasis::Traits::JacobianType> gradphihat;
    std::array<typename LocalBasis::Jacobian, lb.size()> gradphi;
    // loop over quadrature points
    for (const auto& ip : PPS::quadratureRule(geo, order))
    {
    // evaluate basis functions
    lb.evaluateFunction(ip.position(), phihat);
    auto phihat = lb.evaluateFunction(ip.position());
    // evaluate u
    auto u = std::inner_product(begin(x), end(x), begin(phihat), RF(0));
    ......@@ -198,7 +190,7 @@ namespace PPS {
    auto w = std::inner_product(begin(lp), end(lp), begin(phihat), RF(0));
    // evaluate gradient of shape functions
    lb.evaluateJacobian(ip.position(), gradphihat);
    auto gradphihat = lb.evaluateJacobian(ip.position());
    // transform gradients of shape functions to real element
    const auto S = geo.jacobianInverseTransposed(ip.position());
    ......
    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