SubAssembler.hh 2.63 KB
Newer Older
1
2
#include "Quadrature.h"
#include "FiniteElemSpace.h"
3
#include "ElInfo.h"
4
5
6
7

namespace AMDiS {

  template<typename T>
8
9
10
11
  void SubAssembler::getVectorAtQPs(DOFVectorBase<T>* vec, 
				    const ElInfo* elInfo,
				    Quadrature *quad,
				    mtl::dense_vector<T>& vecAtQPs)
12
13
14
15
16
17
18
19
20
  {
    FUNCNAME("SubAssembler::getVectorAtQPs()");

    TEST_EXIT_DBG(vec)("No DOF vector!\n");
    TEST_EXIT_DBG(elInfo->getMesh() == vec->getFeSpace()->getMesh())
      ("Vector and FE space do not fit together!\n");

    Quadrature *localQuad = quad ? quad : quadrature;

21
22
23
24
    if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid) {
      vecAtQPs = valuesAtQPs[vec]->values;
      return;
    }
25
26
27

    if (!valuesAtQPs[vec])
      valuesAtQPs[vec] = new ValuesAtQPs;
28
    valuesAtQPs[vec]->values.change_dim(localQuad->getNumPoints());
Thomas Witkowski's avatar
Thomas Witkowski committed
29
    vecAtQPs.change_dim(localQuad->getNumPoints());
30

31
    mtl::dense_vector<T>& values = valuesAtQPs[vec]->values;
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
    bool sameFeSpaces = 
      vec->getFeSpace() == rowFeSpace || vec->getFeSpace() == colFeSpace;

    if (opt && !quad && sameFeSpaces) {
      const BasisFunction *psi = rowFeSpace->getBasisFcts();
      const BasisFunction *phi = colFeSpace->getBasisFcts();
      if (vec->getFeSpace()->getBasisFcts() == psi)
	psiFast = updateFastQuadrature(psiFast, psi, INIT_PHI);
      else if(vec->getFeSpace()->getBasisFcts() == phi)
	phiFast = updateFastQuadrature(phiFast, phi, INIT_PHI);      
    }

    // calculate new values
    const BasisFunction *basFcts = vec->getFeSpace()->getBasisFcts();

    if (opt && !quad && sameFeSpaces) {
      if (psiFast->getBasisFunctions() == basFcts) {
	vec->getVecAtQPs(elInfo, NULL, psiFast, values);
      } else if (phiFast->getBasisFunctions() == basFcts) {
	vec->getVecAtQPs(elInfo, NULL, phiFast, values);
      } else {
	vec->getVecAtQPs(elInfo, localQuad, NULL, values);
      }
    } else {
      vec->getVecAtQPs(elInfo, localQuad, NULL, values);
    }

    valuesAtQPs[vec]->valid = true;
60
    vecAtQPs = values;
61
62
63
64
  }


  template<typename T>
65
66
67
68
69
  void SubAssembler::getVectorAtQPs(DOFVectorBase<T>* vec, 
				    const ElInfo* smallElInfo,
				    const ElInfo* largeElInfo,
				    Quadrature *quad,
				    mtl::dense_vector<T>& vecAtQPs)
70
71
72
73
74
75
76
77
78
79
  {
    FUNCNAME("SubAssembler::getVectorAtQPs()");

    TEST_EXIT_DBG(vec)("No DOF vector!\n");

    Quadrature *localQuad = quad ? quad : quadrature;

    if (!valuesAtQPs[vec])
      valuesAtQPs[vec] = new ValuesAtQPs;

80
81
    valuesAtQPs[vec]->values.change_dim(localQuad->getNumPoints());
    mtl::dense_vector<T>& values = valuesAtQPs[vec]->values;
82
83
    valuesAtQPs[vec]->valid = true;
    vec->getVecAtQPs(smallElInfo, largeElInfo, localQuad, NULL, values);
84
    vecAtQPs = values;
85
86
87
  }

}