SubAssembler.hh 2.41 KB
Newer Older
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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#include "Quadrature.h"
#include "FiniteElemSpace.h"

namespace AMDiS {

  template<typename T>
  T* SubAssembler::getVectorAtQPs(DOFVectorBase<T>* vec, 
				  const ElInfo* elInfo,
				  Quadrature *quad)
  {
    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;

    if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid) 
      return &(valuesAtQPs[vec]->values[0]);

    if (!valuesAtQPs[vec])
      valuesAtQPs[vec] = new ValuesAtQPs;
    valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());

    T *values = &(valuesAtQPs[vec]->values[0]);
    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;

    return values;
  }


  template<typename T>
  T* SubAssembler::getVectorAtQPs(DOFVectorBase<T>* vec, 
				  const ElInfo* smallElInfo,
				  const ElInfo* largeElInfo,
				  Quadrature *quad)
  {
    FUNCNAME("SubAssembler::getVectorAtQPs()");

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

    Quadrature *localQuad = quad ? quad : quadrature;

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

    valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());
    T *values = &(valuesAtQPs[vec]->values[0]);
    valuesAtQPs[vec]->valid = true;
    vec->getVecAtQPs(smallElInfo, largeElInfo, localQuad, NULL, values);

    return values;
  }

}