<commit_before>// dune-multiscale
// Copyright Holders: <NAME>, <NAME>
// License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

#ifndef DUNE_RIGHT_HAND_SIDE_ASSEMBLER_HH
#define DUNE_RIGHT_HAND_SIDE_ASSEMBLER_HH

#ifdef HAVE_CMAKE_CONFIG
  #include "cmake_config.h"
#else
  #include "config.h"
#endif // ifdef HAVE_CMAKE_CONFIG

// - Dune includes
#include <dune/fem/quadrature/quadrature.hh>
#include <dune/multiscale/hmm/cell_problem_numbering.hh>
#include <dune/multiscale/hmm/cell_problem_solver.hh>

//#include <dune/multiscale/tools/solver/MsFEM/msfem_localproblems/subgrid-list.hh>
#include <dune/multiscale/tools/misc.hh>
#include <dune/stuff/fem/functions/checks.hh>
#include <dune/stuff/common/logging.hh>

namespace Dune {
//! Assembler for right rand side
//! We assemble the right hand side in a LSE, i.e. f \cdot \Phi_H + G \cdot \nabala \Phi_H
//! we call f the first Source and G the second Source
template< class DiscreteFunctionImp >
class RightHandSideAssembler
{
private:
  typedef DiscreteFunctionImp DiscreteFunctionType;
  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType
    DiscreteFunctionSpaceType;
  typedef typename DiscreteFunctionType::LocalFunctionType
    LocalFunctionType;
  typedef typename DiscreteFunctionSpaceType::BasisFunctionSetType
    BasisFunctionSetType;
  typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
  typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
  typedef DomainFieldType TimeType;
  typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
  typedef typename GridPartType::GridType GridType;
  typedef typename DiscreteFunctionSpaceType::JacobianRangeType
    JacobianRangeType;
  typedef typename DiscreteFunctionSpaceType::DomainType
    DomainType;
  typedef typename GridType::template Codim< 0 >::Entity EntityType;
  typedef typename EntityType::Geometry GeometryType;
  typedef Fem::CachingQuadrature< GridPartType, 0 > Quadrature;
  enum { dimension = GridType::dimension };

public:
  static void printRHS(const DiscreteFunctionType& rhs) {
    for (auto entity : rhs.space())
    {
      const LocalFunctionType elementOfRHS = rhs.localFunction(entity);
      const int numDofs = elementOfRHS.numDofs();
      for (int i = 0; i < numDofs; ++i)
      {
        DSC_LOG_DEBUG << "Number of Dof: " << i << " ; " << rhs.name() << " : " << elementOfRHS[i] << std::endl;
      }
    }
  }  // end method

  //! --------------- The rhs-assemble()-methods for linear elliptic problems -----------------

private:
  // need a virtual base to work around local classes not being allowed in templated scopes
  struct FunctorBase {
    virtual RangeType operator()(const DomainType& global_quad_point, const JacobianRangeType& gradientPhi) const = 0;
    virtual ~FunctorBase(){}
  };

  template< class FirstSourceType >
  static void assemble_common(const FirstSourceType& f,
                       const FunctorBase& functor,
                       const int polOrd,
                       DiscreteFunctionType& rhsVector)
  {
    // set rhsVector to zero:
    rhsVector.clear();
    for (const auto& entity : rhsVector.space())
    {
      const GeometryType& geometry = entity.geometry(); // Referenz auf Geometrie
      LocalFunctionType elementOfRHS = rhsVector.localFunction(entity);   // entity zeigt auf ein bestimmtes Element der
                                                                       // entity
      // hier wird sozusagen ein Pointer von localFunction auf discreteFunction erzeugt. Befinden wir uns auf einer
      // bestimmten entity, so berechnet localFunction alle noetigen Werte und speichert sie (da Pointer) in
      // discreteFunction(aktuelleEntity)

      const BasisFunctionSetType baseSet // BaseFunctions leben immer auf Refernzelement!!!
        = rhsVector.space().basisFunctionSet(entity);     // entity Referenz auf eine bestimmtes Element der entity. In der
                                                          // ersten Klasse war das Element fest, deshalb konnte man sich
                                                          // dort Pointer sparen. //loeschen: discreteFunctionSpace
                                                          // statt
                                                          // functionSpace

      const Fem::CachingQuadrature< GridPartType, 0 > quadrature(entity, polOrd);   // 0 --> codim 0
      const int numQuadraturePoints = quadrature.nop();
      const int numDofs = elementOfRHS.numDofs();
      for (int quadraturePoint = 0; quadraturePoint < numQuadraturePoints; ++quadraturePoint)
      {
        // the return values:
        RangeType f_x;
        std::vector<RangeType>phi_x;
        std::vector<JacobianRangeType> gradientPhi;
        // to save: A \nabla PHI_H * \nabla phi_h;
        RangeType res = 0;
        const double det = geometry.integrationElement(quadrature.point(quadraturePoint));

        f.evaluate(geometry.global( quadrature.point(quadraturePoint) ), f_x);
        baseSet.evaluateAll(quadrature[quadraturePoint], phi_x);
        baseSet.jacobianAll(quadrature[quadraturePoint], gradientPhi);

        for (int i = 0; i < numDofs; ++i)
        {
          res = functor(geometry.global(quadrature.point(quadraturePoint)), gradientPhi[i]);
          elementOfRHS[i] += det * quadrature.weight(quadraturePoint) * (f_x * phi_x[i]);
          elementOfRHS[i] += det * quadrature.weight(quadraturePoint) * (res);
        }
      }<commit_msg>where jacobianAll does resize the output arg, evaluateAll does not...<commit_after>// dune-multiscale
// Copyright Holders: <NAME>, <NAME>
// License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

#ifndef DUNE_RIGHT_HAND_SIDE_ASSEMBLER_HH
#define DUNE_RIGHT_HAND_SIDE_ASSEMBLER_HH

#ifdef HAVE_CMAKE_CONFIG
  #include "cmake_config.h"
#else
  #include "config.h"
#endif // ifdef HAVE_CMAKE_CONFIG

// - Dune includes
#include <dune/fem/quadrature/quadrature.hh>
#include <dune/multiscale/hmm/cell_problem_numbering.hh>
#include <dune/multiscale/hmm/cell_problem_solver.hh>

//#include <dune/multiscale/tools/solver/MsFEM/msfem_localproblems/subgrid-list.hh>
#include <dune/multiscale/tools/misc.hh>
#include <dune/stuff/fem/functions/checks.hh>
#include <dune/stuff/common/logging.hh>

namespace Dune {
//! Assembler for right rand side
//! We assemble the right hand side in a LSE, i.e. f \cdot \Phi_H + G \cdot \nabala \Phi_H
//! we call f the first Source and G the second Source
template< class DiscreteFunctionImp >
class RightHandSideAssembler
{
private:
  typedef DiscreteFunctionImp DiscreteFunctionType;
  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType
    DiscreteFunctionSpaceType;
  typedef typename DiscreteFunctionType::LocalFunctionType
    LocalFunctionType;
  typedef typename DiscreteFunctionSpaceType::BasisFunctionSetType
    BasisFunctionSetType;
  typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
  typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
  typedef DomainFieldType TimeType;
  typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
  typedef typename GridPartType::GridType GridType;
  typedef typename DiscreteFunctionSpaceType::JacobianRangeType
    JacobianRangeType;
  typedef typename DiscreteFunctionSpaceType::DomainType
    DomainType;
  typedef typename GridType::template Codim< 0 >::Entity EntityType;
  typedef typename EntityType::Geometry GeometryType;
  typedef Fem::CachingQuadrature< GridPartType, 0 > Quadrature;
  enum { dimension = GridType::dimension };

public:
  static void printRHS(const DiscreteFunctionType& rhs) {
    for (auto entity : rhs.space())
    {
      const LocalFunctionType elementOfRHS = rhs.localFunction(entity);
      const int numDofs = elementOfRHS.numDofs();
      for (int i = 0; i < numDofs; ++i)
      {
        DSC_LOG_DEBUG << "Number of Dof: " << i << " ; " << rhs.name() << " : " << elementOfRHS[i] << std::endl;
      }
    }
  }  // end method

  //! --------------- The rhs-assemble()-methods for linear elliptic problems -----------------

private:
  // need a virtual base to work around local classes not being allowed in templated scopes
  struct FunctorBase {
    virtual RangeType operator()(const DomainType& global_quad_point, const JacobianRangeType& gradientPhi) const = 0;
    virtual ~FunctorBase(){}
  };

  template< class FirstSourceType >
  static void assemble_common(const FirstSourceType& f,
                       const FunctorBase& functor,
                       const int polOrd,
                       DiscreteFunctionType& rhsVector)
  {
    // set rhsVector to zero:
    rhsVector.clear();
    for (const auto& entity : rhsVector.space())
    {
      const GeometryType& geometry = entity.geometry(); // Referenz auf Geometrie
      LocalFunctionType elementOfRHS = rhsVector.localFunction(entity);   // entity zeigt auf ein bestimmtes Element der
                                                                       // entity
      // hier wird sozusagen ein Pointer von localFunction auf discreteFunction erzeugt. Befinden wir uns auf einer
      // bestimmten entity, so berechnet localFunction alle noetigen Werte und speichert sie (da Pointer) in
      // discreteFunction(aktuelleEntity)

      const BasisFunctionSetType baseSet // BaseFunctions leben immer auf Refernzelement!!!
        = rhsVector.space().basisFunctionSet(entity);     // entity Referenz auf eine bestimmtes Element der entity. In der
                                                          // ersten Klasse war das Element fest, deshalb konnte man sich
                                                          // dort Pointer sparen. //loeschen: discreteFunctionSpace
                                                          // statt
                                                          // functionSpace

      const Fem::CachingQuadrature< GridPartType, 0 > quadrature(entity, polOrd);   // 0 --> codim 0
      const int numQuadraturePoints = quadrature.nop();
      const int numDofs = elementOfRHS.numDofs();
      for (int quadraturePoint = 0; quadraturePoint < numQuadraturePoints; ++quadraturePoint)
      {
        // the return values:
        RangeType f_x;
        std::vector<RangeType> phi_x(numDofs);
        std::vector<JacobianRangeType> gradientPhi(numDofs);
        // to save: A \nabla PHI_H * \nabla phi_h;
        RangeType res = 0;
        const double det = geometry.integrationElement(quadrature.point(quadraturePoint));

        f.evaluate(geometry.global( quadrature.point(quadraturePoint) ), f_x);
        baseSet.evaluateAll(quadrature[quadraturePoint], phi_x);
        baseSet.jacobianAll(quadrature[quadraturePoint], gradientPhi);

        for (int i = 0; i < numDofs; ++i)
        {
          res = functor(geometry.global(quadrature.point(quadraturePoint)), gradientPhi[i]);
          elementOfRHS[i] += det * quadrature.weight(quadraturePoint) * (f_x * phi_x[i]);
          elementOfRHS[i] += det * quadrature.weight(quadraturePoint) * (res);
        }
      }
{'fraction_non_alphanumeric': 0.09538378475659139, 'fraction_numerical': 0.0021596328624133897, 'mean_word_length': 2.896914446002805, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6264781', 'n_tokens_mistral': 3227, 'n_tokens_neox': 3122, 'n_words': 802}