Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

GenericOperatorTerm.hh 4.6 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
/******************************************************************************
 *
 * 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 GenericOperatorTerm.hh */

namespace AMDiS {
  
  namespace detail {
    
    template<typename Term, typename T>
    struct GenericOperatorTerm { };
    
    template<typename Term>
    struct GenericOperatorTerm<Term, double> { 
      typedef GenericZeroOrderTerm<Term> type;
    };
    
    template<typename Term>
    struct GenericOperatorTerm<Term, WorldVector<double> > { 
      typedef GenericFirstOrderTerm_b<Term> type;
    };
    
    template<typename Term >
    struct GenericOperatorTerm<Term, WorldMatrix<double> > { 
      typedef GenericSecondOrderTerm_A<Term, true> type;
    };
    
  }
  
  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 boost::is_base_of<LazyOperatorTermBase, Term>::type, typename Term::value_type>::type
56
integrate(Term term, Mesh* mesh_)
Praetorius, Simon's avatar
Praetorius, Simon committed
57 58 59 60 61
{
  typedef typename Term::value_type TOut;
  
  typename GenericOperatorTerm<Term>::type ot(term);
  std::set<const FiniteElemSpace*> feSpaces = ot.getAuxFeSpaces();
62 63 64
  
  TEST_EXIT(mesh_ || !feSpaces.empty())("The Expression must contain a DOFVector or FeSpace depended value!\n");
  Mesh* mesh = mesh_ ? mesh_ : (*feSpaces.begin())->getMesh();
Praetorius, Simon's avatar
Praetorius, Simon committed
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
  
  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;
}


// works only for nodal basis functions!
template<typename Term>
96
inline typename boost::enable_if<typename boost::is_base_of<LazyOperatorTermBase, Term>::type, void>::type
Praetorius, Simon's avatar
Praetorius, Simon committed
97 98 99 100 101
transformDOF(Term term, DOFVector<typename Term::value_type>* result)
{
  typedef typename Term::value_type TOut;
  TOut tmp; nullify(tmp);
  
102 103 104
  DOFVector<TOut> temp(result->getFeSpace(), "temp");
  DOFVector<short int> assigned(result->getFeSpace(), "assigned");
  
Praetorius, Simon's avatar
Praetorius, Simon committed
105
  typename GenericOperatorTerm<Term>::type ot(term);
106
  Mesh* mesh = result->getFeSpace()->getMesh();
Praetorius, Simon's avatar
Praetorius, Simon committed
107
  
108
  const FiniteElemSpace* resultFeSpace = temp.getFeSpace();
Praetorius, Simon's avatar
Praetorius, Simon committed
109 110 111 112
  const BasisFunction *basisFcts = resultFeSpace->getBasisFcts();  
  int nBasisFcts = basisFcts->getNumber();
  
  assigned.set(0);
113
  temp.set(tmp);
Praetorius, Simon's avatar
Praetorius, Simon committed
114 115 116 117 118 119 120 121 122 123 124 125
  
  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++) {
126
      temp[localIndices[i]] += term(i);
Praetorius, Simon's avatar
Praetorius, Simon committed
127 128 129
      assigned[localIndices[i]]++;
    }
    elInfo = stack.traverseNext(elInfo);
130
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
131
  
132 133
  DOFIterator<TOut> tempIter(&temp, USED_DOFS);
  DOFIterator<TOut> resultIter(result, USED_DOFS);
Praetorius, Simon's avatar
Praetorius, Simon committed
134
  DOFIterator<short int> assignedIter(&assigned, USED_DOFS);
135 136
  for (tempIter.reset(), resultIter.reset(), assignedIter.reset(); !resultIter.end(); ++tempIter, ++resultIter, ++assignedIter) {
    *resultIter = (*tempIter)/static_cast<double>(*assignedIter);
Praetorius, Simon's avatar
Praetorius, Simon committed
137 138 139 140 141 142 143 144 145 146 147 148 149 150
  }
}

template<typename Term>
typename boost::enable_if<typename boost::is_base_of<LazyOperatorTermBase, Term>::type, 
			  DOFVector<typename Term::value_type>& >::type
operator<<(DOFVector<typename Term::value_type>& result, const Term& term)
{
  transformDOF(term, &result);
  return result;
}


}