Expressions.hh 6.62 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
136
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;
}



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

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

193
194
195
196
197
198
199
200

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
201
202
203
204
205
206
{
  transformDOF(term, &result);
  return result;
}


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