Expressions.hh 6.58 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
/******************************************************************************
 *
 * 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.
 * 
 ******************************************************************************/



23
/** \file Expressions.hh */
Praetorius, Simon's avatar
Praetorius, Simon committed
24
25
26

namespace AMDiS {
  
27
28
29
30
namespace detail {
  
  template<typename Term, typename T, typename Enable = void>
  struct GenericOperatorTerm { };
Praetorius, Simon's avatar
Praetorius, Simon committed
31
  
32
33
34
  template<typename Term, typename T>
  struct GenericOperatorTerm<Term, T, typename boost::enable_if<traits::is_scalar<T> >::type> { 
    typedef GenericZeroOrderTerm<Term> type;
Praetorius, Simon's avatar
Praetorius, Simon committed
35
  };
36
37
38
39
40
41
42
43
44
45
46
47
  
  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
Praetorius, Simon's avatar
Praetorius, Simon committed
48
49

template<typename Term>
50
51
52
53
54
55
56
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
57
integrate(Term term, Mesh* mesh_)
Praetorius, Simon's avatar
Praetorius, Simon committed
58
59
60
61
62
{
  typedef typename Term::value_type TOut;
  
  typename GenericOperatorTerm<Term>::type ot(term);
  std::set<const FiniteElemSpace*> feSpaces = ot.getAuxFeSpaces();
63
64
65
  
  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
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
  
  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;
}


95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
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)
{
  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;
}



Praetorius, Simon's avatar
Praetorius, Simon committed
136
// works only for nodal basis functions!
137
138
139
140
141
142
143
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)
Praetorius, Simon's avatar
Praetorius, Simon committed
144
145
146
147
{
  typedef typename Term::value_type TOut;
  TOut tmp; nullify(tmp);
  
148
  DOFVector<TOut> temp(result->getFeSpace(), "temp");
149
  DOFVector<int> assigned(result->getFeSpace(), "assigned");
150
  
Praetorius, Simon's avatar
Praetorius, Simon committed
151
  typename GenericOperatorTerm<Term>::type ot(term);
152
  Mesh* mesh = result->getFeSpace()->getMesh();
Praetorius, Simon's avatar
Praetorius, Simon committed
153
  
154
  const FiniteElemSpace* resultFeSpace = temp.getFeSpace();
Praetorius, Simon's avatar
Praetorius, Simon committed
155
156
157
158
  const BasisFunction *basisFcts = resultFeSpace->getBasisFcts();  
  int nBasisFcts = basisFcts->getNumber();
  
  assigned.set(0);
159
  temp.set(tmp);
Praetorius, Simon's avatar
Praetorius, Simon committed
160
161
162
163
164
165
166
167
  
  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) {
168
    term.initElement(&ot, elInfo, NULL, NULL, basisFcts);
Praetorius, Simon's avatar
Praetorius, Simon committed
169
170
171
    basisFcts->getLocalIndices(elInfo->getElement(), resultFeSpace->getAdmin(), localIndices);	
    
    for (int i = 0; i < nBasisFcts; i++) {
172
      temp[localIndices[i]] += term(i);
Praetorius, Simon's avatar
Praetorius, Simon committed
173
174
175
      assigned[localIndices[i]]++;
    }
    elInfo = stack.traverseNext(elInfo);
176
  }
177
178
179
180
181
182

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  Parallel::MeshDistributor::globalMeshDistributor->synchAddVector(temp);
  Parallel::MeshDistributor::globalMeshDistributor->synchAddVector(assigned);
#endif
  
Praetorius, Simon's avatar
Praetorius, Simon committed
183
  
184
  DOFIterator<TOut> tempIter(&temp, USED_DOFS);
185
186
  DOFIterator<T> resultIter(result, USED_DOFS);
  DOFIterator<int> assignedIter(&assigned, USED_DOFS);
187
  for (tempIter.reset(), resultIter.reset(), assignedIter.reset(); !resultIter.end(); ++tempIter, ++resultIter, ++assignedIter) {
188
    *resultIter = (*tempIter) * (1.0/static_cast<double>(*assignedIter));
Praetorius, Simon's avatar
Praetorius, Simon committed
189
190
191
  }
}

192
193
194
195
196
197
198
199

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)
Praetorius, Simon's avatar
Praetorius, Simon committed
200
201
202
203
204
205
{
  transformDOF(term, &result);
  return result;
}


206
207
208
209
210
211
212
213
214
215
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