valueOf.hpp 15 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * 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.
18
 *
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
 ******************************************************************************/



/** \file valueOf.hpp */

#ifndef AMDIS_VALUE_OF_HPP
#define AMDIS_VALUE_OF_HPP

#include "AMDiS_fwd.h"
#include "LazyOperatorTerm.h"
#include "DOFVector.h"

#include "traits/category.hpp"
#include "traits/at.hpp"
34
#include "traits/resize.hpp"
35
36
37
38
39

#include <boost/static_assert.hpp>
#include <boost/mpl/bool.hpp>


40
namespace AMDiS
41
42
43
{
  struct _unknown {};

44
  namespace expressions
45
46
47
48
  {
    /// Expressions that extracts the values of a DOFVector at QPs
    template<typename Vector, typename Name, typename Enable = void>
    struct ValueOf : public LazyOperatorTermBase {};
49

50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
    template<typename T, typename Name>
    struct ValueOf<DOFVector<T>, Name> : public LazyOperatorTermBase
    {
      typedef T value_type;
      typedef Name  id;

      DOFVector<T>* vecDV;
      mutable mtl::dense_vector<T> vec;
      mutable mtl::dense_vector<T> coeff;

      ValueOf(DOFVector<T>& vector) : vecDV(&vector) {}
      ValueOf(DOFVector<T>* vector) : vecDV(vector) {}

      template<typename List>
      inline void insertFeSpaces(List& feSpaces) const
      {
	feSpaces.insert(vecDV->getFeSpace());
      }
68

69
70
71
72
73
74
75
      inline int getDegree() const
      {
	return vecDV->getFeSpace()->getBasisFcts()->getDegree();
      }

      template<typename OT>
      inline void initElement(OT* ot, const ElInfo* elInfo,
76
			      SubAssembler* subAssembler, Quadrature *quad,
77
78
79
80
81
82
83
84
			      const BasisFunction *basisFct = NULL)
      {
	if (ot && subAssembler)
	  ot->getVectorAtQPs(vecDV, elInfo, subAssembler, quad, vec);
	else if (quad)
	  vecDV->getVecAtQPs(elInfo, quad, NULL, vec);
	else if (basisFct) {
	  const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
85

86
87
88
	  // get coefficients of DOFVector
	  coeff.change_dim(localBasisFct->getNumber());
	  vecDV->getLocalVector(elInfo->getElement(), coeff);
89

90
91
	  // eval basisfunctions of DOFVector at coords of given basisFct
	  size_t nBasisFct = basisFct->getNumber();
92
	  vec.change_dim(nBasisFct);
93
94
95
96
97
98
99
100
	  for (size_t i = 0; i < nBasisFct; i++)
	    vec[i] = localBasisFct->evalUh(*basisFct->getCoords(i), coeff);
	}
      }


      template<typename OT>
      inline void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
101
			      SubAssembler* subAssembler, Quadrature *quad,
102
103
104
105
106
107
108
109
			      const BasisFunction *basisFct = NULL)
      {
	if (ot && subAssembler)
	  ot->getVectorAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
	else if (quad)
	  vecDV->getVecAtQPs(smallElInfo, largeElInfo, quad, NULL, vec);
	else if (basisFct) {
	  const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
110

111
112
113
	  // get coefficients of DOFVector
	  coeff.change_dim(localBasisFct->getNumber());
	  vecDV->getLocalVector(smallElInfo->getElement(), coeff);
114

115
116
	  // eval basisfunctions of DOFVector at coords of given basisFct
	  size_t nBasisFct = basisFct->getNumber();
117
	  vec.change_dim(nBasisFct);
118
119
120
121
122
123
	  for (size_t i = 0; i < nBasisFct; i++)
	    vec[i] = localBasisFct->evalUh(*basisFct->getCoords(i), coeff);
	}
      }

      inline value_type operator()(const int& iq) const { return vec[iq]; }
124

125
126
      std::string str() const { return std::string("value(") + vecDV->getName() + ")"; }
    };
127
128


129
130
    /// Expressions that extracts the matrix-value of a Matrix<DOFVector> at QPs
    template<template<class> class Matrix, typename T, typename Name>
131
132
    struct ValueOf<Matrix<DOFVector<T>*>, Name,
		  typename enable_if< traits::is_matrix<Matrix<T> > >::type >
133
134
135
136
137
138
139
140
141
      : public LazyOperatorTermBase
    {
      typedef Matrix<T> value_type;
      typedef Name  id;

      Matrix<DOFVector<T>*> vecDV;
      mutable mtl::dense_vector<value_type> vec;
      mutable Matrix<mtl::dense_vector<T> > coeff;

142
      ValueOf(Matrix<DOFVector<T>*> const& vector) : vecDV(vector)
143
144
145
146
147
148
149
150
151
152
153
      {
	resize(coeff, num_rows(vecDV), num_cols(vecDV));
      }

      template<typename List>
      inline void insertFeSpaces(List& feSpaces) const
      {
	for (size_t i = 0; i < num_rows(vecDV); i++)
	for (size_t j = 0; j < num_cols(vecDV); j++)
	  feSpaces.insert(at(vecDV, i, j)->getFeSpace());
      }
154

155
156
157
158
159
160
161
      inline int getDegree() const
      {
	return at(vecDV, 0, 0)->getFeSpace()->getBasisFcts()->getDegree();
      }

      template<typename OT>
      inline void initElement(OT* ot, const ElInfo* elInfo,
162
			      SubAssembler* subAssembler, Quadrature *quad,
163
164
165
166
167
168
169
170
			      const BasisFunction *basisFct = NULL)
      {
	Matrix<mtl::dense_vector<T> > helper; resize(helper, num_rows(vecDV), num_cols(vecDV));
	for (size_t i = 0; i < num_rows(vecDV); i++) {
	for (size_t j = 0; j < num_cols(vecDV); j++) {
	  if (ot && subAssembler)
	    ot->getVectorAtQPs(at(vecDV, i, j), elInfo, subAssembler, quad, at(helper, i, j));
	  else if (quad)
171
	    at(vecDV, i, j)->getVecAtQPs(elInfo, quad, NULL, at(helper, i, j));
172
173
	  else if (basisFct) {
	    const BasisFunction *localBasisFct = at(vecDV, i, j)->getFeSpace()->getBasisFcts();
174

175
176
177
	    // get coefficients of DOFVector
	    at(coeff, i, j).change_dim(localBasisFct->getNumber());
	    at(vecDV, i, j)->getLocalVector(elInfo->getElement(), at(coeff, i, j));
178

179
180
	    // eval basisfunctions of DOFVector at coords of given basisFct
	    size_t nBasisFct = basisFct->getNumber();
181
	    at(helper, i, j).change_dim(nBasisFct);
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
	    for (size_t k = 0; k < nBasisFct; k++)
	      at(helper, i, j)[k] = localBasisFct->evalUh(*basisFct->getCoords(k), at(coeff, i, j));
	  }
	}
	}
	vec.change_dim(num_rows(at(helper, 0, 0)));
	for (size_t iq = 0; iq < num_rows(at(helper, 0, 0)); iq++) {
	  value_type tmp; resize(tmp, num_rows(vecDV), num_cols(vecDV));
	  for (size_t i = 0; i < num_rows(vecDV); i++) {
	  for (size_t j = 0; j < num_cols(vecDV); j++) {
	    at(tmp, i, j) = at(helper, i, j)[iq];
	  }
	  }
	  vec[iq] = tmp;
	}
      }


      template<typename OT>
      inline void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
202
			      SubAssembler* subAssembler, Quadrature *quad,
203
			      const BasisFunction *basisFct = NULL)
204
      {
205
206
207
208
	initElement(ot, smallElInfo, subAssembler, quad, basisFct);
      }

      inline value_type operator()(const int& iq) const { return vec[iq]; }
209

210
211
      std::string str() const { return std::string("value_(") + at(vecDV, 0, 0)->getName() + ")"; }
    };
212
213


214
215
    /// Expressions that extracts the vector-value of a Vector<DOFVector> at QPs
    template<template<class> class Vector, typename T, typename Name>
216
217
    struct ValueOf<Vector<DOFVector<T>*>, Name,
		  typename boost::enable_if<typename traits::is_vector<Vector<T> >::type>::type >
218
219
220
221
222
223
224
225
226
      : public LazyOperatorTermBase
    {
      typedef Vector<T> value_type;
      typedef Name  id;

      Vector<DOFVector<T>*> vecDV;
      mutable mtl::dense_vector<value_type> vec;
      mutable Vector<mtl::dense_vector<T> > coeff;

227
      ValueOf(Vector<DOFVector<T>*>& vector) : vecDV(vector)
228
229
230
231
232
233
234
235
236
237
      {
	resize(coeff, num_rows(vecDV));
      }

      template<typename List>
      inline void insertFeSpaces(List& feSpaces) const
      {
	for (size_t i = 0; i < num_rows(vecDV); i++)
	  feSpaces.insert(at(vecDV, i)->getFeSpace());
      }
238

239
240
241
242
243
244
245
      inline int getDegree() const
      {
	return at(vecDV, 0)->getFeSpace()->getBasisFcts()->getDegree();
      }

      template<typename OT>
      inline void initElement(OT* ot, const ElInfo* elInfo,
246
			      SubAssembler* subAssembler, Quadrature *quad,
247
248
249
250
251
252
253
			      const BasisFunction *basisFct = NULL)
      {
	Vector<mtl::dense_vector<T> > helper; resize(helper, num_rows(vecDV));
	for (size_t i = 0; i < num_rows(vecDV); i++) {
	  if (ot && subAssembler)
	    ot->getVectorAtQPs(at(vecDV, i), elInfo, subAssembler, quad, at(helper, i));
	  else if (quad)
254
	    at(vecDV, i)->getVecAtQPs(elInfo, quad, NULL, at(helper, i));
255
256
	  else if (basisFct) {
	    const BasisFunction *localBasisFct = at(vecDV, i)->getFeSpace()->getBasisFcts();
257

258
259
260
	    // get coefficients of DOFVector
	    at(coeff, i).change_dim(localBasisFct->getNumber());
	    at(vecDV, i)->getLocalVector(elInfo->getElement(), at(coeff, i));
261

262
263
	    // eval basisfunctions of DOFVector at coords of given basisFct
	    size_t nBasisFct = basisFct->getNumber();
264
	    at(helper, i).change_dim(nBasisFct);
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
	    for (size_t j = 0; j < nBasisFct; j++)
	      at(helper, i)[j] = localBasisFct->evalUh(*basisFct->getCoords(j), at(coeff, i));
	  }
	}
	vec.change_dim(num_rows(at(helper, 0)));
	for (size_t iq = 0; iq < num_rows(at(helper, 0)); iq++) {
	  value_type tmp; resize(tmp, num_rows(vecDV));
	  for (size_t i = 0; i < num_rows(vecDV); i++)
	    at(tmp, i) = at(helper, i)[iq];
	  vec[iq] = tmp;
	}
      }


      template<typename OT>
      inline void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
281
			      SubAssembler* subAssembler, Quadrature *quad,
282
			      const BasisFunction *basisFct = NULL)
283
      {
284
285
286
287
	initElement(ot, smallElInfo, subAssembler, quad, basisFct);
      }

      inline value_type operator()(const int& iq) const { return vec[iq]; }
288

289
290
      std::string str() const { return std::string("value_(") + at(vecDV, 0)->getName() + ")"; }
    };
291
292


293
294
295
    /// Expression that extracts the component of a vector-values DOFVector at QPs
    template<typename Vector>
    struct ComponentOf : public LazyOperatorTermBase {};
296

297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
    template<template<class> class Vector, typename T>
    struct ComponentOf<DOFVector<Vector<T> > > : public LazyOperatorTermBase
    {
      typedef T value_type;

      DOFVector<Vector<T> >* vecDV;
      mutable mtl::dense_vector<Vector<T> > vec;
      mutable mtl::dense_vector<Vector<T> > coeff;
      int I;

      ComponentOf(DOFVector<Vector<T> >& vector, int I_) : vecDV(&vector), I(I_) {}
      ComponentOf(DOFVector<Vector<T> >* vector, int I_) : vecDV(vector), I(I_) {}

      template<typename List>
      inline void insertFeSpaces(List& feSpaces) const
      {
	feSpaces.insert(vecDV->getFeSpace());
      }
315

316
317
318
319
320
321
322
      inline int getDegree() const
      {
	return vecDV->getFeSpace()->getBasisFcts()->getDegree();
      }

      template<typename OT>
      inline void initElement(OT* ot, const ElInfo* elInfo,
323
			      SubAssembler* subAssembler, Quadrature *quad,
324
325
326
327
328
329
330
331
			      const BasisFunction *basisFct = NULL)
      {
	if (ot && subAssembler)
	  ot->getVectorAtQPs(vecDV, elInfo, subAssembler, quad, vec);
	else if (quad)
	  vecDV->getVecAtQPs(elInfo, quad, NULL, vec);
	else if (basisFct) {
	  const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
332

333
334
335
	  // get coefficients of DOFVector
	  coeff.change_dim(localBasisFct->getNumber());
	  vecDV->getLocalVector(elInfo->getElement(), coeff);
336

337
338
	  // eval basisfunctions of DOFVector at coords of given basisFct
	  size_t nBasisFct = basisFct->getNumber();
339
	  vec.change_dim(nBasisFct);
340
341
342
343
344
345
346
347
	  for (size_t i = 0; i < nBasisFct; i++)
	    vec[i] = localBasisFct->evalUh(*basisFct->getCoords(i), coeff);
	}
      }


      template<typename OT>
      inline void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
348
			      SubAssembler* subAssembler, Quadrature *quad,
349
350
351
352
353
354
355
356
			      const BasisFunction *basisFct = NULL)
      {
	if (ot && subAssembler)
	  ot->getVectorAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
	else if (quad)
	  vecDV->getVecAtQPs(smallElInfo, largeElInfo, quad, NULL, vec);
	else if (basisFct) {
	  const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
357

358
359
360
	  // get coefficients of DOFVector
	  coeff.change_dim(localBasisFct->getNumber());
	  vecDV->getLocalVector(smallElInfo->getElement(), coeff);
361

362
363
	  // eval basisfunctions of DOFVector at coords of given basisFct
	  size_t nBasisFct = basisFct->getNumber();
364
	  vec.change_dim(nBasisFct);
365
366
367
368
369
370
	  for (size_t i = 0; i < nBasisFct; i++)
	    vec[i] = localBasisFct->evalUh(*basisFct->getCoords(i), coeff);
	}
      }

      inline value_type operator()(const int& iq) const { return vec[iq][I]; }
371

372
373
374
375
376
377
378
379
380
381
382
      std::string str() const { return std::string("comp<") + boost::lexical_cast<std::string>(I) + ">(" + vecDV->getName() + ")"; }
    };


  } // end namespace expressions

  // value of a DOFVector<T>
  // ___________________________________________________________________________

  // with Name
  template<typename Name, typename T>
383
  expressions::ValueOf<DOFVector<T>, Name > valueOf(DOFVector<T>& vector)
384
385
386
  { return expressions::ValueOf<DOFVector<T>, Name >(vector); }

  template<typename Name, typename T>
387
  expressions::ValueOf<DOFVector<T>, Name > valueOf(DOFVector<T>* vector)
388
389
390
  { return expressions::ValueOf<DOFVector<T>, Name >(vector); }

  template<typename Name, template<class> class Matrix, typename T>
391
  typename boost::enable_if<typename traits::is_matrix<Matrix<T> >::type,
392
    expressions::ValueOf<Matrix<DOFVector<T>*>, Name > >::type
393
  valueOf(Matrix<DOFVector<T>*> &mat)
394
395
396
  { return expressions::ValueOf<Matrix<DOFVector<T>*>, Name >(mat); }

  template<typename Name, template<class> class Vector, typename T>
397
  typename boost::enable_if<typename traits::is_vector<Vector<T> >::type,
398
    expressions::ValueOf<Vector<DOFVector<T>*>, Name > >::type
399
  valueOf(Vector<DOFVector<T>*> &vector)
400
401
402
403
404
  { return expressions::ValueOf<Vector<DOFVector<T>*>, Name >(vector); }


  // without Name
  template<typename T>
405
  expressions::ValueOf<DOFVector<T>, _unknown > valueOf(DOFVector<T>& vector)
406
407
408
  { return expressions::ValueOf<DOFVector<T>, _unknown >(vector); }

  template<typename T>
409
  expressions::ValueOf<DOFVector<T>, _unknown > valueOf(DOFVector<T>* vector)
410
411
412
  { return expressions::ValueOf<DOFVector<T>, _unknown >(vector); }

  template<template<class> class Matrix, typename T>
413
  typename boost::enable_if<typename traits::is_matrix<Matrix<T> >::type,
414
    expressions::ValueOf<Matrix<DOFVector<T>*>, _unknown > >::type
415
  valueOf(Matrix<DOFVector<T>*> &mat)
416
417
418
  { return expressions::ValueOf<Matrix<DOFVector<T>*>, _unknown >(mat); }

  template<template<class> class Vector, typename T>
419
  typename boost::enable_if<typename traits::is_vector<Vector<T> >::type,
420
    expressions::ValueOf<Vector<DOFVector<T>*>, _unknown > >::type
421
  valueOf(Vector<DOFVector<T>*> &vector)
422
423
  { return expressions::ValueOf<Vector<DOFVector<T>*>, _unknown >(vector); }

424

425
426
427
428
  // component of a DOFVector<Vector>
  // ___________________________________________________________________________

  template<template<class> class Vector, typename T>
429
  expressions::ComponentOf<DOFVector<Vector<T> > > componentOf(DOFVector<Vector<T> >& vector, int I)
430
431
432
  { return expressions::ComponentOf<DOFVector<Vector<T> > >(vector, I); }

  template<template<class> class Vector, typename T>
433
  expressions::ComponentOf<DOFVector<Vector<T> > > componentOf(DOFVector<Vector<T> >* vector, int I)
434
435
436
437
438
  { return expressions::ComponentOf<DOFVector<Vector<T> > >(vector, I); }

} // end namespace AMDiS

#endif // AMDIS_VALUE_OF_HPP