SubElementAssembler.cc 5.21 KB
Newer Older
1
#include <vector>
2
#include <boost/numeric/mtl/mtl.hpp>
3
4
5
6
#include "SubElementAssembler.h"
#include "ScalableQuadrature.h"
#include "SubPolytope.h"

7
8
9
10
11
12
13
14
15
16
17
namespace AMDiS {

  SubElementAssembler::SubElementAssembler(Operator *op, 
					   const FiniteElemSpace *rowFESpace_,
					   const FiniteElemSpace *colFESpace_)
    : StandardAssembler(op, NULL, NULL, NULL, NULL, rowFESpace_, colFESpace_)
  {
    /** 
     * Create a scalable quadrature for subassembler and replace the original 
     * quadrature of the subassembler with the scalable quadrature.
     */
18

19
20
21
22
23
24
25
26
    if (zeroOrderAssembler) {
      checkQuadratures();
      zeroOrderScalableQuadrature = 
	NEW ScalableQuadrature(zeroOrderAssembler->getQuadrature());
      zeroOrderAssembler->setQuadrature(zeroOrderScalableQuadrature);
    } else {
      zeroOrderScalableQuadrature = NULL;
    }
27

28
29
30
31
32
33
34
35
    if (firstOrderAssemblerGrdPsi) {
      checkQuadratures();
      firstOrderGrdPsiScalableQuadrature = 
	NEW ScalableQuadrature(firstOrderAssemblerGrdPsi->getQuadrature());
      firstOrderAssemblerGrdPsi->setQuadrature(firstOrderGrdPsiScalableQuadrature);
    } else {
      firstOrderGrdPsiScalableQuadrature = NULL;
    }
36

37
38
39
40
41
42
43
44
    if (firstOrderAssemblerGrdPhi) {
      checkQuadratures();
      firstOrderGrdPhiScalableQuadrature = 
	NEW ScalableQuadrature(firstOrderAssemblerGrdPhi->getQuadrature());
      firstOrderAssemblerGrdPhi->setQuadrature(firstOrderGrdPhiScalableQuadrature);
    } else {
      firstOrderGrdPhiScalableQuadrature = NULL;
    }
45

46
47
48
49
50
51
52
53
    if (secondOrderAssembler) {
      checkQuadratures();
      secondOrderScalableQuadrature = 
	NEW ScalableQuadrature(secondOrderAssembler->getQuadrature());
      secondOrderAssembler->setQuadrature(secondOrderScalableQuadrature);
    } else {
      secondOrderScalableQuadrature = NULL;
    }
54
  }
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

  void SubElementAssembler::scaleQuadratures(const SubElInfo& subElInfo)
  {
    if (zeroOrderAssembler) {
      zeroOrderScalableQuadrature->scaleQuadrature(subElInfo);
    }
    if (firstOrderAssemblerGrdPsi) {
      firstOrderGrdPsiScalableQuadrature->scaleQuadrature(subElInfo);
    }
    if (firstOrderAssemblerGrdPhi) {
      firstOrderGrdPhiScalableQuadrature->scaleQuadrature(subElInfo);
    }
    if (secondOrderAssembler) {
      secondOrderScalableQuadrature->scaleQuadrature(subElInfo);
    }
70
71
  }

72
73
  void SubElementAssembler::getSubElementVector(SubElInfo *subElInfo, 
						const ElInfo *elInfo, 
74
						ElementVector& userVec)
75
76
77
78
79
80
81
82
83
84
85
86
87
  {
    /** 
     * Manipulate the quadratures of the SubAssemblers for subelement.
     */
    scaleQuadratures(*subElInfo);

    calculateElementVector(elInfo, userVec);

    /** 
     * The integration has been performed with a quadrature living on element. The 
     * determinant of element has been used instead of the determinant of subelement. Thus
     * the result must be corrected with respect to subelement.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
88
    double corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
89
90
    for (int i = 0; i < nRow; i++)
      userVec[i] *= corrFactor;
91
92
  }

93
94
  void SubElementAssembler::getSubElementMatrix(SubElInfo *subElInfo, 
						const ElInfo *elInfo, 
95
						ElementMatrix& userMat)
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
  {
    /** 
     * Manipulate the quadratures of the SubAssemblers for subelement.
     */
    scaleQuadratures(*subElInfo);

    /** 
     * Integrate using the manipulated quadrature.
     */
    calculateElementMatrix(elInfo, userMat);

    /** 
     * The integration has been performed with a quadrature living on element. The 
     * determinant of element has been used instead of the determinant of subelement. 
     * Thus the result must be corrected with respect to subelement.
     */
    double corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
115
	userMat[i][j] *= corrFactor;
116
      }
117
118
119
    }
  }

120
121
122
  void SubElementAssembler::getSubPolytopeVector(SubPolytope *subPolytope,
						 SubElementAssembler *subElementAssembler,
						 const ElInfo *elInfo,
123
						 ElementVector& subPolVec)
124
  {
125
    /// Note: There is no reset of subPolVec.
126
    std::vector<SubElInfo *>::iterator it;
127
    ElementVector subElVec(nRow);
128

129
    /// Assemble for each subelement of subpolytope.
130
131
132
    for (it = subPolytope->getSubElementsBegin(); 
	 it != subPolytope->getSubElementsEnd();
	 it++) {
133
      set_to_zero(subElVec);
134
135
      subElementAssembler->getSubElementVector(*it, elInfo, subElVec);

136
137
      /// Add results for subelement to total result for subpolytope.
      subPolVec += subElVec;
138
    }
139
  }
140

141
142
143
  void SubElementAssembler::getSubPolytopeMatrix(SubPolytope *subPolytope,
						 SubElementAssembler *subElementAssembler,
						 const ElInfo *elInfo,
144
						 ElementMatrix& subPolMat)
145
146
147
148
149
  {
    /**
     * Note: There is no reset of subPolMat.
     */
    std::vector<SubElInfo *>::iterator it;
150
    ElementMatrix subElMat(nRow, nCol);
151
152

    /**
153
     * Assemble for each subelement of subpolytope.
154
     */
155
156
157
    for (it = subPolytope->getSubElementsBegin(); 
	 it != subPolytope->getSubElementsEnd();
	 it++) {
158
      set_to_zero(subElMat);
159
160
161
162
163
      subElementAssembler->getSubElementMatrix(*it, elInfo, subElMat);

      /**
       * Add results for subelement to total result for subpolytope.
       */
164
165
166
      for (int i = 0; i < nRow; i++)
	for (int j = 0; j < nCol; j++)
	  subPolMat[i][j] += subElMat[i][j];
167
168
169
170
    }
  }

}