SubElementAssembler.cc 5.46 KB
Newer Older
1
2
3
4
5
6
7
#include <vector>
#include "ElementMatrix.h"
#include "ElementVector.h"
#include "SubElementAssembler.h"
#include "ScalableQuadrature.h"
#include "SubPolytope.h"

8
9
10
11
12
13
14
15
16
17
18
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.
     */
19

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

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

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

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

  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);
    }
71
72
  }

73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
  void SubElementAssembler::getSubElementVector(SubElInfo *subElInfo, 
						const ElInfo *elInfo, 
						ElementVector *userVec)
  {
    double corrFactor;

    /** 
     * 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.
     */
    corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
    for (int i = 0; i < nRow; i++) {
      (*userVec)[i] *= corrFactor;
    }
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
  void SubElementAssembler::getSubElementMatrix(SubElInfo *subElInfo, 
						const ElInfo *elInfo, 
						ElementMatrix *userMat)
  {
    /** 
     * 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++) {
	(*userMat)[i][j] *= corrFactor;
      }
121
122
123
    }
  }

124
125
126
127
128
129
130
131
132
133
  void SubElementAssembler::getSubPolytopeVector(SubPolytope *subPolytope,
						 SubElementAssembler *subElementAssembler,
						 const ElInfo *elInfo,
						 ElementVector *subPolVec)
  {
    /**
     * Note: There is no reset of subPolVec.
     */
    std::vector<SubElInfo *>::iterator it;
    ElementVector *subElVec = NEW ElementVector(nRow);
134
135

    /**
136
     * Assemble for each subelement of subpolytope.
137
     */
138
139
140
141
142
143
144
145
146
147
148
149
    for (it = subPolytope->getSubElementsBegin(); 
	 it != subPolytope->getSubElementsEnd();
	 it++) {
      subElVec->set(0.0);
      subElementAssembler->getSubElementVector(*it, elInfo, subElVec);

      /**
       * Add results for subelement to total result for subpolytope.
       */
      for (int i = 0; i < nRow; i++) {
	(*subPolVec)[i] += (*subElVec)[i];
      }
150
151
    }

152
153
    DELETE subElVec;
  }
154

155
156
157
158
159
160
161
162
163
164
  void SubElementAssembler::getSubPolytopeMatrix(SubPolytope *subPolytope,
						 SubElementAssembler *subElementAssembler,
						 const ElInfo *elInfo,
						 ElementMatrix *subPolMat)
  {
    /**
     * Note: There is no reset of subPolMat.
     */
    std::vector<SubElInfo *>::iterator it;
    ElementMatrix *subElMat = NEW ElementMatrix(nRow, nCol);
165
166

    /**
167
     * Assemble for each subelement of subpolytope.
168
     */
169
170
171
172
173
174
175
176
177
178
179
180
181
    for (it = subPolytope->getSubElementsBegin(); 
	 it != subPolytope->getSubElementsEnd();
	 it++) {
      subElMat->set(0.0);
      subElementAssembler->getSubElementMatrix(*it, elInfo, subElMat);

      /**
       * Add results for subelement to total result for subpolytope.
       */
      for (int i = 0; i < nRow; i++) {
	for (int j = 0; j < nCol; j++) {
	  (*subPolMat)[i][j] += (*subElMat)[i][j];
	}
182
183
      }
    }
184
185

    DELETE subElMat;
186
187
188
  }

}