test1.cc

#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/common/Literals.hpp>

#define TEST_EQ(expr1, expr2) test_exit( (expr1) == (expr2), "'" , #expr1 , "' != '" , #expr2 , "'" )

#ifndef DIM
#define DIM 2
#endif
#ifndef DOW
#define DOW 2
#endif

using namespace AMDiS;
namespace {

  using TestParam   = ProblemStatTraits<DIM, DOW, 1>;
  using TestProblem = ProblemStat<TestParam>;

  using MeshView = typename TestProblem::MeshView;
  using Geometry = typename MeshView::template Codim<0>::Geometry;
  using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>;



  template <class FeSpace>
  class LocalEvaluation
  {
  public:
    LocalEvaluation(FeSpace const& feSpace)
      : localView(feSpace.localView())
    {}

    template <class Element>
    void bind(Element const& element)
    {
      localView.bind(element);
    }

    /// evaluate DOFVector at barycentric coordinates \p lambda in element that
    /// is bind to instance in \ref bind().
    template <class ValueType, class LocalCoord>
    ValueType eval(DOFVector<FeSpace, ValueType> const& vec, LocalCoord const& lambda)
    {
      auto const& localBasis = localView.tree().finiteElement().localBasis();
      localBasis.evaluateFunction(lambda, shapeValues);

      ValueType data = 0;
      for (std::size_t j = 0; j < shapeValues.size(); ++j) {
        const auto global_idx = localView.index(j);
        data += vec[global_idx] * shapeValues[j];
      }

      return data;
    }

  private:
    typename FeSpace::LocalView localView;

    std::vector<Dune::FieldVector<double,1> > shapeValues;
  };


  class AMDiSTester
  {
  public:
    // create a problemStat
    AMDiSTester(std::string name)
      : prob(make_shared<TestProblem>(name))
    {
      prob->initialize(INIT_ALL);
    }


    void test1() // CreateDOFVector
    {
      using FeSpace = typename TestProblem::template FeSpace<0>;

      // construction of DOFVectors
      DOFVector<FeSpace, double> vec1(prob->getFeSpace(0_c), "vec1");
      auto vec2 = makeDOFVector(prob->getFeSpace(0_c), "vec2");

      // retrieving DOFVectors from problemStat
      auto vec3 = prob->getSolution(0_c);
      auto vec5 = prob->getSolution<0>();
      auto vec6 = prob->getSolution(index_<0>());

      auto vec7 = prob->getSolution()->getDOFVector(0_c);
      auto vec8 = prob->getSolution()->getDOFVector<0>();
      auto vec9 = prob->getSolution()->getDOFVector(index_<0>());

      auto& sys_vec = *prob->getSolution();
      auto vec10 = sys_vec[0_c];
      auto vec11 = sys_vec[index_<0>()];

      // construction of SystemVector
      using FeSpaces = typename TestProblem::FeSpaces;
      SystemVector<FeSpaces, double> sys_vec1(*prob->getFeSpaces(), prob->getComponentNames());
      auto sys_vec2 = makeSystemVector(*prob->getFeSpaces(), prob->getComponentNames());
      auto sys_vec3 = makeSystemVector(*prob->getFeSpaces(), "sys_vec");

      // retrieving systemVector from problemStat
      auto sys_vec4 = *prob->getSolution();
    }


    void test2() // FillDOFVector
    {
      auto& vec = prob->getSolution(0_c);
      auto const& feSpace = vec.getFeSpace();

      // interpolate function to dofvector
      vec.interpol([](auto const& x) { return x*x; });

      // test whether interpolation is fine
      using FeSpace = std::decay_t< decltype(feSpace) >;
      LocalEvaluation<FeSpace> evaluator(feSpace);

      for (auto const& element : elements(feSpace.gridView())) {
        evaluator.bind(element);

        auto const& refElement = RefElements::general(element.type());

        std::size_t nVertices = element.subEntities(DIM);
        for (std::size_t i = 0; i < nVertices; ++i) {
          auto pos = refElement.position(i, DIM);
          auto data = evaluator.eval(vec, pos);

          auto x = element.geometry().global(pos);
          TEST_EQ(data, x*x);
        }
      }
    }

  protected:
    shared_ptr<TestProblem> prob;
  };

} // end namespace

int main(int argc, char** argv)
{
  AMDiS::init(argc, argv);

  auto tester = make_shared<AMDiSTester>("test");
  tester->test1();
  tester->test2();

  AMDiS::finalize();
  return 0;
}