ZeroOrderAssembler.cc 9.52 KB
Newer Older
1
#include <vector>
2
#include <boost/numeric/mtl/mtl.hpp>
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
#include "Assembler.h"
#include "ZeroOrderAssembler.h"
#include "Operator.h"
#include "QPsiPhi.h"
#include "FiniteElemSpace.h"
#include "Quadrature.h"
#include "DOFVector.h"
#include "OpenMP.h"

namespace AMDiS {

  std::vector<SubAssembler*> ZeroOrderAssembler::optimizedSubAssemblers;
  std::vector<SubAssembler*> ZeroOrderAssembler::standardSubAssemblers;

  ZeroOrderAssembler::ZeroOrderAssembler(Operator *op,
					 Assembler *assembler,
					 Quadrature *quad,
					 bool optimized)
    : SubAssembler(op, assembler, quad, 0, optimized)
  {}

  ZeroOrderAssembler* ZeroOrderAssembler::getSubAssembler(Operator* op,
							  Assembler *assembler,
							  Quadrature *quad,
							  bool optimized)
  {
29
30
    int myRank = omp_get_thread_num();

31
    // check if a assembler is needed at all
32
    if (op->zeroOrder[myRank].size() == 0) {
33
      return NULL;   
34
    }
35
36
37
38

    ZeroOrderAssembler *newAssembler;

    std::vector<SubAssembler*> *subAssemblers =
39
      optimized ? &optimizedSubAssemblers : &standardSubAssemblers;
40

41
    std::vector<OperatorTerm*> opTerms = op->zeroOrder[myRank];
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

    sort(opTerms.begin(), opTerms.end());

    // check if a new assembler is needed
    if (quad) {
      for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
	std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());

	sort(assTerms.begin(), assTerms.end());

	if ((opTerms == assTerms) && 
	    ((*subAssemblers)[i]->getQuadrature() == quad)) {
	  return dynamic_cast<ZeroOrderAssembler*>((*subAssemblers)[i]);
	}
      }
    }
 
    // check if all terms are pw_const
    bool pwConst = true;
61
    for (int i = 0; i < static_cast<int>(op->zeroOrder[myRank].size()); i++) {
62
63
64
65
66
67
68
      if (!op->zeroOrder[myRank][i]->isPWConst()) {
	pwConst = false;
	break;
      }
    }  

    // create new assembler
69
70
71
72
73
74
75
76
77
    if (!optimized) {
	newAssembler = NEW StandardZOA(op, assembler, quad);
    } else {
      if (pwConst) {
	  newAssembler = NEW PrecalcZOA(op, assembler, quad);
      } else {
	  newAssembler = NEW FastQuadZOA(op, assembler, quad);
      }
    }
78
79
80
81
82
83
84

    subAssemblers->push_back(newAssembler);
    return newAssembler;
  }

  StandardZOA::StandardZOA(Operator *op, Assembler *assembler, Quadrature *quad)
    : ZeroOrderAssembler(op, assembler, quad, false)
Thomas Witkowski's avatar
Thomas Witkowski committed
85
  {}
86

87
  void StandardZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
88
89
90
91
92
  {
    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();

    int nPoints = quadrature->getNumPoints();
93
94
    std::vector<double> c(nPoints, 0.0);
    std::vector<double> phival(nCol);
95
96
97
98
99
100
101
102

    int myRank = omp_get_thread_num();
    std::vector<OperatorTerm*>::iterator termIt;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
    }
      
    if (symmetric) {
103
104
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
            
105
106
107
108
109
110
111
112
113
114
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	// calculate phi at QPs only once!
	for (int i = 0; i < nCol; i++) {
	  phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
	}

	for (int i = 0; i < nRow; i++) {
	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
115
	  mat[i][i] += quadrature->getWeight(iq) * c[iq] * psival * phival[i];
116
117
	  for (int j = i + 1; j < nCol; j++) {
	    double val = quadrature->getWeight(iq) * c[iq] * psival * phival[j];
118
119
	    mat[i][j] += val;
	    mat[j][i] += val;
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
	  }
	}
      }
    } else {      //  non symmetric assembling 
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	// calculate phi at QPs only once!
	for (int i = 0; i < nCol; i++) {
	  phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
	}

	for (int i = 0; i < nRow; i++) {
	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
	  for (int j = 0; j < nCol; j++) {
135
	    mat[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j];
136
137
138
139
140
141
	  }
	}
      }
    }
  }

142
  void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
143
144
  {
    int nPoints = quadrature->getNumPoints();
145
    std::vector<double> c(nPoints, 0.0);
146
147
148
149
150
151
152
153
154
155
156
157
158

    int myRank = omp_get_thread_num();
    std::vector<OperatorTerm*>::iterator termIt;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
    }

    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] *= elInfo->getDet();

      for (int i = 0; i < nRow; i++) {
	double psi = (*(owner->getRowFESpace()->getBasisFcts()->getPhi(i)))
	  (quadrature->getLambda(iq));
159
	vec[i] += quadrature->getWeight(iq) * c[iq] * psi;
160
161
162
163
164
165
      }
    }
  }

  FastQuadZOA::FastQuadZOA(Operator *op, Assembler *assembler, Quadrature *quad)
    : ZeroOrderAssembler(op, assembler, quad, true)
166
  {}
167

168
  void FastQuadZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
169
170
171
172
173
  {
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();

    if (firstCall) {
174
175
176
177
178
179
180
181
182
183
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
	const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
	psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
	basFcts = owner->getColFESpace()->getBasisFcts();
	phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
	firstCall = false;
      }
184
185
    }

186
    std::vector<double> c(nPoints, 0.0);
187
188
189
190
191
192
    std::vector<OperatorTerm*>::iterator termIt;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
    }

    if (symmetric) {
193
194
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");

195
196
197
198
199
200
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	const double *psi = psiFast->getPhi(iq);
	const double *phi = phiFast->getPhi(iq);
	for (int i = 0; i < nRow; i++) {
201
	  mat[i][i] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[i];
202
203
	  for (int j = i + 1; j < nCol; j++) {
	    double val = quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
204
205
	    mat[i][j] += val;
	    mat[j][i] += val;
206
207
208
209
210
211
212
213
214
215
216
	  }
	}
      }
    } else {      /*  non symmetric assembling   */
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	const double *psi = psiFast->getPhi(iq);
	const double *phi = phiFast->getPhi(iq);
	for (int i = 0; i < nRow; i++) {
	  for (int j = 0; j < nCol; j++) {
217
	    mat[i][j] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
218
219
220
221
222
	  }
	}
      }
    }
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
223

224
  void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
225
  {
226
227
228
    int myRank = omp_get_thread_num();
    int nPoints = quadrature->getNumPoints();

229
    if (firstCall) {
230
231
232
233
234
235
236
237
238
239
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
	const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
	psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
	basFcts = owner->getColFESpace()->getBasisFcts();
	phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
	firstCall = false;
      }
240
241
    }

242
    std::vector<double> c(nPoints, 0.0);
243
244
245
246
247
248
249
250
251
252
    std::vector<OperatorTerm*>::iterator termIt;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
    }

    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] *= elInfo->getDet();

      const double *psi = psiFast->getPhi(iq);
      for (int i = 0; i < nRow; i++) {
253
	vec[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
254
255
256
257
258
259
260
261
262
      }
    }
  }

  PrecalcZOA::PrecalcZOA(Operator *op, Assembler *assembler, Quadrature *quad) 
    : ZeroOrderAssembler(op, assembler, quad, true)
  {
  }

263
  void PrecalcZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
264
265
  {
    if (firstCall) {
266
267
268
269
270
271
272
273
274
275
276
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
	q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					  owner->getColFESpace()->getBasisFcts(), 
					  quadrature);
	q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
				 quadrature);
	firstCall = false;
      }
277
278
    }

279
    std::vector<double> c(1, 0.0);
280
281
282
283
    int myRank = omp_get_thread_num();
    int size = static_cast<int>(terms[myRank].size());

    for (int i = 0; i < size; i++) {
284
      (static_cast<ZeroOrderTerm*>((terms[myRank][i])))->getC(elInfo, 1, c);
285
286
    }

287
    c[0] *= elInfo->getDet();
288
    
289
290
    if (symmetric) {
      for (int i = 0; i < nRow; i++) {
291
	mat[i][i] += c[0] * q00->getValue(i,i);
292
	for (int j = i + 1; j < nCol; j++) {
293
	  double val = c[0] * q00->getValue(i, j);
294
295
	  mat[i][j] += val;
	  mat[j][i] += val;
296
297
298
	}
      }
    } else {
299
      for (int i = 0; i < nRow; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
300
	for (int j = 0; j < nCol; j++) {
301
	  mat[i][j] += c[0] * q00->getValue(i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
302
	}
303
      }
304
305
306
    }
  }

307
  void PrecalcZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
308
309
  {
    if (firstCall) {
310
311
312
313
314
315
316
317
318
319
320
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
	q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					  owner->getColFESpace()->getBasisFcts(), 
					  quadrature);
	q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
				 quadrature);	
	firstCall = false;
      }
321
322
323
324
325
    }

    std::vector<OperatorTerm*>::iterator termIt;

    int myRank = omp_get_thread_num();
326
    std::vector<double> c(1, 0.0);
327
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
328
      (static_cast<ZeroOrderTerm*>( *termIt))->getC(elInfo, 1, c);
329

330
    c[0] *= elInfo->getDet();
331
332

    for (int i = 0; i < nRow; i++)
333
      vec[i] += c[0] * q0->getValue(i);
334
335
336
  }

}