Operator.inc.hpp

#pragma once

namespace AMDiS
{
  template <class MeshView>
    template <class RowFeSpace, class ColFeSpace>
  void Operator<MeshView>::init(RowFeSpace const& rowFeSpace,
				ColFeSpace const& colFeSpace)
  {
    AMDIS_FUNCNAME("Operator::init()");

    using IdType = typename Operator<MeshView>::IdType;
    lastMatrixId = std::numeric_limits<IdType>::max();
    lastVectorId = std::numeric_limits<IdType>::max();

//     auto const& rowFE = rowView.tree().finiteElement();
//     auto const& colFE = colView.tree().finiteElement();

//     psiDegree = rowFE.localBasis().order();
//     phiDegree = colFE.localBasis().order();

    psiDegree = getPolynomialDegree<RowFeSpace>;
    phiDegree = getPolynomialDegree<ColFeSpace>;

    // TODO: calc quadrature degree here.
  }


  template <class MeshView>
    template <class Element>
  int Operator<MeshView>::getQuadratureDegree(Element const& element, int order, FirstOrderType firstOrderType)
  {
    std::list<OperatorTermType*>* terms = NULL;

    auto const& t = element.type();
    auto geometry = element.geometry();

    switch(order)
    {
    case 0:
      terms = &zeroOrder;
      break;
    case 1:
      if (firstOrderType == GRD_PHI)
        terms = &firstOrderGrdPhi;
      else
        terms = &firstOrderGrdPsi;
      break;
    case 2:
      terms = &secondOrder;
      break;
    }
    int maxTermDegree = 0;

    for (OperatorTermType* term : *terms)
      maxTermDegree = std::max(maxTermDegree, term->getDegree(t));

    int degree = psiDegree + phiDegree + maxTermDegree;
    if (t.isSimplex())
      degree -= order;

    if (!geometry.affine())
      degree += 1; // oder += (order+1)

    return degree;
  }


  template <class MeshView>
    template <class LocalView>
  typename Operator<MeshView>::IdType
  Operator<MeshView>::getElementId(LocalView const& localView)
  {
    auto const& element = localView.element();
    auto const& grid = localView.globalBasis().gridView().grid();

    auto const& idSet = grid.localIdSet();
    return idSet.id(element);
  }


  template <class MeshView>
    template <class RowView, class ColView>
  bool Operator<MeshView>::getElementMatrix(RowView const& rowView,
                                            ColView const& colView,
                                            ElementMatrix& elementMatrix,
                                            double* factor)
  {
    AMDIS_FUNCNAME("Operator::getElementMatrix()");

    double fac = factor ? *factor : 1.0;

    if (fac == 0.0)
      return false;

    const auto nRows = rowView.tree().finiteElement().size();
    const auto nCols = colView.tree().finiteElement().size();

    auto currentId = getElementId(rowView);
    bool useCachedMatrix = (lastMatrixId == currentId
                            && num_rows(cachedElementMatrix) == nRows
                            && num_cols(cachedElementMatrix) == nCols);

    bool assign = true;
    if (!useCachedMatrix) {

      cachedElementMatrix.change_dim(nRows, nCols);
      set_to_zero(cachedElementMatrix);

      if (!zeroOrder.empty())
        assembleZeroOrderTerms(rowView, colView, cachedElementMatrix);
      if (!firstOrderGrdPhi.empty())
        assembleFirstOrderTermsGrdPhi(rowView, colView, cachedElementMatrix);
      if (!firstOrderGrdPsi.empty())
        assembleFirstOrderTermsGrdPsi(rowView, colView, cachedElementMatrix);
      if (!secondOrder.empty())
        assembleSecondOrderTerms(rowView, colView, cachedElementMatrix);

      assign = !zeroOrder.empty() || !firstOrderGrdPhi.empty() ||
               !firstOrderGrdPsi.empty() || !secondOrder.empty();
    }

    if (assign)
      elementMatrix += fac * cachedElementMatrix;

    lastMatrixId = currentId;

    return true;
  }


  template <class MeshView>
    template <class RowView>
  bool Operator<MeshView>::getElementVector(RowView const& rowView,
					                                  ElementVector& elementVector,
                                            double* factor)
  {
    AMDIS_FUNCNAME("Operator::getElementVector()");

    double fac = factor ? *factor : 1.0;

    if (fac == 0.0)
      return false;

    test_exit( firstOrderGrdPhi.empty() && secondOrder.empty(),
               "Derivatives on ansatz-functions not allowed on the vector-side!");

    const auto nRows = rowView.tree().finiteElement().size();

    auto currentId = getElementId(rowView);
    bool useCachedVector = (lastVectorId == currentId && size(cachedElementVector) == nRows);

    bool assign = true;
    if (!useCachedVector) {
      cachedElementVector.change_dim(nRows);
      set_to_zero(cachedElementVector);

      if (!zeroOrder.empty())
        assembleZeroOrderTerms(rowView, cachedElementVector);
      if (!firstOrderGrdPsi.empty())
        assembleFirstOrderTermsGrdPsi(rowView, cachedElementVector);

      assign = !zeroOrder.empty() || !firstOrderGrdPsi.empty();
    }

    if (assign)
      elementVector += fac * cachedElementVector;

    lastVectorId = currentId;

    return true;
  }


  template <class MeshView>
    template <class RowView, class ColView>
  void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
                                                  ColView const& colView,
                                                  ElementMatrix& elementMatrix)
  {
    AMDIS_FUNCNAME("Operator::assembleZeroOrderTerms(elementMatrix)");

    using Element = typename RowView::Element;
    auto element = rowView.element();
    const int dim = Element::dimension;
    auto geometry = element.geometry();

    auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
    auto const& colLocalBasis = colView.tree().finiteElement().localBasis();

    int order = getQuadratureDegree(element, 0);
    auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);

    for (auto* operatorTerm : zeroOrder)
      operatorTerm->init(element, quad);

    for (std::size_t iq = 0; iq < quad.size(); ++iq) {
      // Position of the current quadrature point in the reference element
      Dune::FieldVector<double,dim> const& quadPos = quad[iq].position();

      // The multiplicative factor in the integral transformation formula
      const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();

      std::vector<Dune::FieldVector<double,1> > rowShapeValues, colShapeValues;
      rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);

      if (psiDegree == phiDegree)
        colShapeValues = rowShapeValues;
      else
        colLocalBasis.evaluateFunction(quadPos, colShapeValues);

      for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
        for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
          const int local_i = rowView.tree().localIndex(i);
          const int local_j = colView.tree().localIndex(j);
          for (auto* operatorTerm : zeroOrder)
            elementMatrix[local_i][local_j]
              += operatorTerm->evalZot(iq, rowShapeValues[i], colShapeValues[j]) * factor;
        }
      }
    }
  }



  template <class MeshView>
    template <class RowView>
  void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
						                                      ElementVector& elementvector)
  {
    AMDIS_FUNCNAME("Operator::assembleZeroOrderTerms(elementvector)");

    using Element = typename RowView::Element;
    auto element = rowView.element();
    const int dim = Element::dimension;
    auto geometry = element.geometry();

    auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();

    int order = getQuadratureDegree(element, 0);
    auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);

    for (auto* operatorTerm : zeroOrder)
	    operatorTerm->init(element, quad);

    for (std::size_t iq = 0; iq < quad.size(); ++iq) {
      // Position of the current quadrature point in the reference element
      const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();

      // The multiplicative factor in the integral transformation formula
      const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();

      std::vector<Dune::FieldVector<double,1> > rowShapeValues;
      rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);

      for (std::size_t i = 0; i < size(elementvector); ++i) {
        const int local_i = rowView.tree().localIndex(i);
        for (auto* operatorTerm : zeroOrder)
          elementvector[local_i] += operatorTerm->evalZot(iq, rowShapeValues[i]) * factor;
      }
    }
  }


  template <class MeshView>
    template <class RowView, class ColView>
  void Operator<MeshView>::assembleFirstOrderTermsGrdPhi(RowView const& rowView,
                                                         ColView const& colView,
                                                         ElementMatrix& elementMatrix)
  {
    AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPhi(elementMatrix)");

    auto element = rowView.element();
    auto geometry = element.geometry();
    const int dim = RowView::GridView::dimension;
    const int dow = RowView::GridView::dimensionworld;

    auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
    auto const& colLocalBasis = colView.tree().finiteElement().localBasis();

    int order = getQuadratureDegree(element, 1, GRD_PHI);
    auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);

    for (auto* operatorTerm : firstOrderGrdPhi)
      operatorTerm->init(element, quad);

    for (std::size_t iq = 0; iq < quad.size(); ++iq) {
      // Position of the current quadrature point in the reference element
      const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();

      // The transposed inverse Jacobian of the map from the reference element to the element
      const auto jacobian = geometry.jacobianInverseTransposed(quadPos);

      // The multiplicative factor in the integral transformation formula
      const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();

      std::vector<Dune::FieldVector<double,1> > rowShapeValues;
      rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);

      // The gradients of the shape functions on the reference element
      std::vector<Dune::FieldMatrix<double,1,dim> > colReferenceGradients;
      colLocalBasis.evaluateJacobian(quadPos, colReferenceGradients);

      // Compute the shape function gradients on the real element
      std::vector<Dune::FieldVector<double,dow> > colGradients(colReferenceGradients.size());

      for (std::size_t i = 0; i < colGradients.size(); ++i)
        jacobian.mv(colReferenceGradients[i][0], colGradients[i]);

      for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
        for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
          const int local_i = rowView.tree().localIndex(i);
          const int local_j = colView.tree().localIndex(j);
          for (auto* operatorTerm : firstOrderGrdPhi)
            elementMatrix[local_i][local_j]
              += operatorTerm->evalFot1(iq, rowShapeValues[i], colGradients[j]) * factor;
        }
      }
    }
  }


  template <class MeshView>
    template <class RowView, class ColView>
  void Operator<MeshView>::assembleFirstOrderTermsGrdPsi(RowView const& rowView,
                                                         ColView const& colView,
                                                         ElementMatrix& elementMatrix)
  {
    AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPsi(elementMatrix)");

    auto element = rowView.element();
    auto geometry = element.geometry();
    const int dim = RowView::GridView::dimension;
    const int dow = RowView::GridView::dimensionworld;

    auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
    auto const& colLocalBasis = colView.tree().finiteElement().localBasis();

    int order = getQuadratureDegree(element, 1, GRD_PSI);
    auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);

    for (auto* operatorTerm : firstOrderGrdPsi)
      operatorTerm->init(element, quad);

    for (std::size_t iq = 0; iq < quad.size(); ++iq) {
      // Position of the current quadrature point in the reference element
      const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();

      // The transposed inverse Jacobian of the map from the reference element to the element
      const auto jacobian = geometry.jacobianInverseTransposed(quadPos);

      // The multiplicative factor in the integral transformation formula
      const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();

      // The gradients of the shape functions on the reference element
      std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
      rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);

      // Compute the shape function gradients on the real element
      std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());

      for (std::size_t i = 0; i < rowGradients.size(); ++i)
        jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);

      std::vector<Dune::FieldVector<double,1> > colShapeValues;
      colLocalBasis.evaluateFunction(quadPos, colShapeValues);

      for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
        for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
          const int local_i = rowView.tree().localIndex(i);
          const int local_j = colView.tree().localIndex(j);
          for (auto* operatorTerm : firstOrderGrdPsi)
            elementMatrix[local_i][local_j]
              += operatorTerm->evalFot2(iq, rowGradients[i], colShapeValues[j]) * factor;
        }
      }
    }
  }



  template <class MeshView>
    template <class RowView>
  void Operator<MeshView>::assembleFirstOrderTermsGrdPsi(RowView const& rowView,
                                                         ElementVector& elementvector)
  {
    AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPsi(elementvector)");

    auto element = rowView.element();
    const int dim = RowView::GridView::dimension;
    const int dow = RowView::GridView::dimensionworld;
    auto geometry = element.geometry();

    auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();

    int order = getQuadratureDegree(element, 1, GRD_PSI);
    auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);

    for (auto* operatorTerm : firstOrderGrdPsi)
        operatorTerm->init(element, quad);

    for (std::size_t iq = 0; iq < quad.size(); ++iq) {
      // Position of the current quadrature point in the reference element
      const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();

      // The transposed inverse Jacobian of the map from the reference element to the element
      const auto jacobian = geometry.jacobianInverseTransposed(quadPos);

      // The multiplicative factor in the integral transformation formula
      const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();

      // The gradients of the shape functions on the reference element
      std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
      rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);

      // Compute the shape function gradients on the real element
      std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());

      for (std::size_t i = 0; i < rowGradients.size(); ++i)
        jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);

      for (std::size_t i = 0; i < size(elementvector); ++i) {
        const int local_i = rowView.tree().localIndex(i);
        for (auto* operatorTerm : firstOrderGrdPsi)
          elementvector[local_i] += operatorTerm->evalFot2(iq, rowGradients[i]) * factor;
      }
    }
  }


  template <class MeshView>
    template <class RowView, class ColView>
  void Operator<MeshView>::assembleSecondOrderTerms(RowView const& rowView,
                                                    ColView const& colView,
                                                    ElementMatrix& elementMatrix)
  {
    AMDIS_FUNCNAME("Operator::assembleSecondOrderTerms(elementMatrix)");

    auto element = rowView.element();
    const int dim = RowView::GridView::dimension;
    const int dow = RowView::GridView::dimensionworld;
    auto geometry = element.geometry();

    auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
//     auto const& colLocalBasis = colView.tree().finiteElement().localBasis();

    int order = getQuadratureDegree(element, 2);
    auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);

    for (auto* operatorTerm : secondOrder)
      operatorTerm->init(element, quad);

    // TODO: currently only the implementation for equal fespaces
    assert( psiDegree == phiDegree );

    for (std::size_t iq = 0; iq < quad.size(); ++iq) {
      // Position of the current quadrature point in the reference element
      const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();

      // The transposed inverse Jacobian of the map from the reference element to the element
      const auto jacobian = geometry.jacobianInverseTransposed(quadPos);

      // The multiplicative factor in the integral transformation formula
      const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();

      // The gradients of the shape functions on the reference element
      std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients;
      rowLocalBasis.evaluateJacobian(quadPos, referenceGradients);

      // Compute the shape function gradients on the real element
      std::vector<Dune::FieldVector<double,dow> > gradients(referenceGradients.size());

      for (std::size_t i = 0; i < gradients.size(); ++i)
        jacobian.mv(referenceGradients[i][0], gradients[i]);

      for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
        for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
          const int local_i = rowView.tree().localIndex(i);
          const int local_j = colView.tree().localIndex(j);
          for (auto* operatorTerm : secondOrder)
            elementMatrix[local_i][local_j]
              += operatorTerm->evalSot(iq, gradients[i], gradients[j]) * factor;
        }
      }
    }
  }



  template <class MeshView>
    template <class Term>
  Operator<MeshView>& Operator<MeshView>::addZOTImpl(Term const& term)
  {
    zeroOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
    return *this;
  }


  template <class MeshView>
    template <class Term>
  Operator<MeshView>& Operator<MeshView>::addFOTImpl(Term const& term,
                                                     FirstOrderType firstOrderType)
  {
    using OpTerm = GenericOperatorTerm<MeshView, Term>;
    if (firstOrderType == GRD_PHI)
      firstOrderGrdPhi.push_back(new OpTerm(term));
    else
      firstOrderGrdPsi.push_back(new OpTerm(term));
    return *this;
  }


  template <class MeshView>
    template <class Term>
  Operator<MeshView>& Operator<MeshView>::addFOTImpl(Term const& term,
                                                     std::size_t i,
                                                     FirstOrderType firstOrderType)
  {
    using OpTerm = GenericOperatorTerm<MeshView, Term, VectorComponent>;

    if (firstOrderType == GRD_PHI)
      firstOrderGrdPhi.push_back(new OpTerm(term, {i}));
    else
      firstOrderGrdPsi.push_back(new OpTerm(term, {i}));
    return *this;
  }


  template <class MeshView>
    template <class Term>
  Operator<MeshView>& Operator<MeshView>::addSOTImpl(Term const& term)
  {
    secondOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
    return *this;
  }


  template <class MeshView>
    template <class Term>
  Operator<MeshView>& Operator<MeshView>::addSOTImpl(Term const& term,
                                                     std::size_t i, std::size_t j)
  {
    using OpTerm = GenericOperatorTerm<MeshView, Term, MatrixComponent>;
    secondOrder.push_back(new OpTerm(term, {i,j}));
    return *this;
  }

} // end namespace AMDiS