Expressions.hh

/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/



/** \file Expressions.hh */

namespace AMDiS {
  
namespace detail {
  
  template<typename Term, typename T, typename Enable = void>
  struct GenericOperatorTerm { };
  
  template<typename Term, typename T>
  struct GenericOperatorTerm<Term, T, typename boost::enable_if<traits::is_scalar<T> >::type> { 
    typedef GenericZeroOrderTerm<Term> type;
  };
  
  template<typename Term, typename T>
  struct GenericOperatorTerm<Term, T, typename boost::enable_if<traits::is_vector<T> >::type> { 
    typedef GenericFirstOrderTerm_b<Term> type;
  };
  
  template<typename Term, typename T >
  struct GenericOperatorTerm<Term, T, typename boost::enable_if<traits::is_matrix<T> >::type> { 
    typedef GenericSecondOrderTerm_A<Term, true> type;
  };
  
} // end namespace detail

template<typename Term>
struct GenericOperatorTerm { 
  typedef typename detail::GenericOperatorTerm<Term, typename Term::value_type>::type type;
};


template<typename Term>
inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
integrate(Term term, Mesh* mesh_)
{
  typedef typename Term::value_type TOut;
  
  typename GenericOperatorTerm<Term>::type ot(term);
  std::set<const FiniteElemSpace*> feSpaces = ot.getAuxFeSpaces();
  
  TEST_EXIT(mesh_ || !feSpaces.empty())("The Expression must contain a DOFVector or FeSpace depended value!\n");
  Mesh* mesh = mesh_ ? mesh_ : (*feSpaces.begin())->getMesh();
  
  int deg = term.getDegree();
  int dim = mesh->getDim();
  Quadrature* quad = Quadrature::provideQuadrature(dim, deg);
  
  TOut value; nullify(value);

  Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_DET;
  TraverseStack stack;
  ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
  while (elInfo) {
    term.initElement(&ot, elInfo, NULL, quad);
    TOut tmp; nullify(tmp);
    for (int iq = 0; iq < quad->getNumPoints(); iq++) {
      tmp += quad->getWeight(iq) * term(iq);
    }
    value += tmp * elInfo->getDet();

    elInfo = stack.traverseNext(elInfo);
  }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  Parallel::mpi::globalAdd(value);
#endif

  return value;
}


template<typename Term, typename Functor>
inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
accumulate(Term term, Functor f, typename Term::value_type value0)
{
  typedef typename Term::value_type TOut;
  typename GenericOperatorTerm<Term>::type ot(term);
  std::set<const FiniteElemSpace*> feSpaces = ot.getAuxFeSpaces();
  
  TEST_EXIT(!feSpaces.empty())("The Expression must contain a DOFVector or FeSpace depended value!\n");
  const FiniteElemSpace* feSpace0 = *feSpaces.begin();  
  Mesh* mesh = feSpace0->getMesh();
  
  const BasisFunction *basisFcts = feSpace0->getBasisFcts();  
  int nBasisFcts = basisFcts->getNumber();
  
  DOFVector<bool> assigned(feSpace0, "assigned");
  assigned.set(false);
  
  std::vector<DegreeOfFreedom> localIndices(nBasisFcts);
  TraverseStack stack;
  ElInfo *elInfo = stack.traverseFirst(mesh, -1,
					Mesh::CALL_LEAF_EL | 
					Mesh::FILL_COORDS | Mesh::FILL_GRD_LAMBDA);
  
  while (elInfo) {
    term.initElement(&ot, elInfo, NULL, NULL, basisFcts);
    basisFcts->getLocalIndices(elInfo->getElement(), feSpace0->getAdmin(), localIndices);	
    
    for (int i = 0; i < nBasisFcts; i++) {
      if (!assigned[localIndices[i]]) {
	value0 = f(value0, term(i));
	assigned[localIndices[i]] = true;
      }
    }
    elInfo = stack.traverseNext(elInfo);
  }
  
  return value0;
}



// works only for nodal basis functions!
template<typename T, typename Term>
inline typename boost::enable_if<
  typename boost::mpl::and_<typename traits::is_expr<Term>::type, 
			    typename boost::is_convertible<typename Term::value_type, T>::type
			    >::type
  >::type
transformDOF(Term term, DOFVector<T>* result)
{
  typedef typename Term::value_type TOut;
  TOut tmp; nullify(tmp);
  
  DOFVector<TOut> temp(result->getFeSpace(), "temp");
  DOFVector<int> assigned(result->getFeSpace(), "assigned");
  
  typename GenericOperatorTerm<Term>::type ot(term);
  Mesh* mesh = result->getFeSpace()->getMesh();
  
  const FiniteElemSpace* resultFeSpace = temp.getFeSpace();
  const BasisFunction *basisFcts = resultFeSpace->getBasisFcts();  
  int nBasisFcts = basisFcts->getNumber();
  
  assigned.set(0);
  temp.set(tmp);
  
  std::vector<DegreeOfFreedom> localIndices(nBasisFcts);
  TraverseStack stack;
  ElInfo *elInfo = stack.traverseFirst(mesh, -1,
					Mesh::CALL_LEAF_EL | 
					Mesh::FILL_COORDS | Mesh::FILL_GRD_LAMBDA);
  
  while (elInfo) {
    term.initElement(&ot, elInfo, NULL, NULL, basisFcts);
    basisFcts->getLocalIndices(elInfo->getElement(), resultFeSpace->getAdmin(), localIndices);	
    
    for (int i = 0; i < nBasisFcts; i++) {
      temp[localIndices[i]] += term(i);
      assigned[localIndices[i]]++;
    }
    elInfo = stack.traverseNext(elInfo);
  }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  Parallel::MeshDistributor::globalMeshDistributor->synchAddVector(temp);
  Parallel::MeshDistributor::globalMeshDistributor->synchAddVector(assigned);
#endif
  
  
  DOFIterator<TOut> tempIter(&temp, USED_DOFS);
  DOFIterator<T> resultIter(result, USED_DOFS);
  DOFIterator<int> assignedIter(&assigned, USED_DOFS);
  for (tempIter.reset(), resultIter.reset(), assignedIter.reset(); !resultIter.end(); ++tempIter, ++resultIter, ++assignedIter) {
    *resultIter = (*tempIter) * (1.0/static_cast<double>(*assignedIter));
  }
}


template<typename T, typename Term>
typename boost::enable_if<
  typename boost::mpl::and_<typename traits::is_expr<Term>::type, 
			    typename boost::is_convertible<typename Term::value_type, T>::type
			    >::type,
  DOFVector<T>& >::type
operator<<(DOFVector<T>& result, const Term& term)
{
  transformDOF(term, &result);
  return result;
}


template<typename Term>
typename boost::enable_if<typename traits::is_expr<Term>::type, 
			  std::ostream& >::type
operator<<(std::ostream& result, const Term& term)
{
  result << term.str();
  return result;
}

} // end namespace AMDiS