Assembler.cc 11.9 KB
Newer Older
1
2
3
#include <vector>
#include <algorithm>
#include <boost/numeric/mtl/mtl.hpp>
4
5
6
7
8
#include "Assembler.h"
#include "Operator.h"
#include "Element.h"
#include "QPsiPhi.h"
#include "DOFVector.h"
9
#include "OpenMP.h"
10
11
12

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
13
14
15
  Assembler::Assembler(Operator *op,
		       const FiniteElemSpace *row,
		       const FiniteElemSpace *col) 
16
    : operat(op),
Thomas Witkowski's avatar
Thomas Witkowski committed
17
18
      rowFESpace(row),
      colFESpace(col ? col : row),
19
20
21
22
23
      nRow(rowFESpace->getBasisFcts()->getNumber()),
      nCol(colFESpace->getBasisFcts()->getNumber()),
      remember(true),
      rememberElMat(false),
      rememberElVec(false),
24
25
      elementMatrix(nRow, nCol),
      elementVector(nRow),
26
27
28
29
      lastMatEl(NULL),
      lastVecEl(NULL),
      lastTraverseId(-1)
  
30
  {}
Thomas Witkowski's avatar
Thomas Witkowski committed
31
32

  Assembler::~Assembler()
33
  {}
34
35

  void Assembler::calculateElementMatrix(const ElInfo *elInfo, 
36
					 ElementMatrix& userMat,
37
38
39
40
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

41
    if (remember && (factor != 1.0 || operat->uhOld))
42
      rememberElMat = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
43

44
    Element *el = elInfo->getElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
45

46
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
47
48
49
      initElement(elInfo);

    if (el != lastMatEl || !operat->isOptimized()) {
50
51
52
      if (rememberElMat)
	set_to_zero(elementMatrix);

53
54
55
      lastMatEl = el;
    } else {
      if (rememberElMat) {
56
	userMat += factor * elementMatrix;
57
58
59
	return;
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
60
 
61
    ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
62
63
64
65
66
67
68
69
70
71

    if (secondOrderAssembler)
      secondOrderAssembler->calculateElementMatrix(elInfo, mat);
    if (firstOrderAssemblerGrdPsi)
      firstOrderAssemblerGrdPsi->calculateElementMatrix(elInfo, mat);
    if (firstOrderAssemblerGrdPhi)
      firstOrderAssemblerGrdPhi->calculateElementMatrix(elInfo, mat);
    if (zeroOrderAssembler)
      zeroOrderAssembler->calculateElementMatrix(elInfo, mat);

Thomas Witkowski's avatar
Thomas Witkowski committed
72
    if (rememberElMat && &userMat != &elementMatrix)
73
      userMat += factor * elementMatrix;
74
75
  }

76
77
  void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
					 const ElInfo *colElInfo,
78
79
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
80
					 ElementMatrix& userMat,
81
82
83
84
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

85
    if (remember && ((factor != 1.0) || (operat->uhOld)))
86
87
      rememberElMat = true;
  
Thomas Witkowski's avatar
Thomas Witkowski committed
88
    Element *el = smallElInfo->getElement();   
89
    lastVecEl = lastMatEl = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
90
   
91
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
Thomas Witkowski's avatar
Thomas Witkowski committed
92
      initElement(smallElInfo, largeElInfo);
93
94

    if (el != lastMatEl || !operat->isOptimized()) {
95
96
97
      if (rememberElMat)
	set_to_zero(elementMatrix);

98
99
100
      lastMatEl = el;
    } else {
      if (rememberElMat) {
101
	userMat += factor * elementMatrix;
102
103
104
105
	return;
      }
    }
  
106
    ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
107

108
    if (secondOrderAssembler) {
109
      secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
110
111
112
113
114
115
116
117
118
119

      ElementMatrix m(nRow, nRow);
      smallElInfo->getSubElemCoordsMat_so(m, rowFESpace->getBasisFcts()->getDegree());
      ElementMatrix tmpMat(nRow, nRow);
      tmpMat = m * mat;
      mat = tmpMat;
    }

    if (firstOrderAssemblerGrdPsi) {
      ERROR_EXIT("Not yet implemented for not zero order assembler!\n");
120
      firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
121
122
123
124
    }

    if (firstOrderAssemblerGrdPhi) {
      ERROR_EXIT("Not yet implemented for not zero order assembler!\n");
125
      firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
126
    }
127

128
129
    if (zeroOrderAssembler) {
      zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
130

131
132
133
134
135
136
      ElementMatrix m(nRow, nRow);
      smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
      ElementMatrix tmpMat(nRow, nRow);
      tmpMat = m * mat;
      mat = tmpMat;
    }
137

138
139
    if (rememberElMat)
      userMat += factor * elementMatrix;
140
141
  }

142
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
143
					 ElementVector& userVec,
144
145
146
147
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

148
    if (remember && factor != 1.0)
149
150
151
152
      rememberElVec = true;

    Element *el = elInfo->getElement();

153
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
154
      initElement(elInfo);
155
    
Thomas Witkowski's avatar
Thomas Witkowski committed
156
    if (el != lastVecEl || !operat->isOptimized()) {
157
158
159
      if (rememberElVec)
	set_to_zero(elementVector);
	
160
161
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
162
      if (rememberElVec) {
163
	userVec += factor * elementVector;
164
165
166
	return;
      }
    }
167
168

    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
169
    if (operat->uhOld && remember) {
170
      matVecAssemble(elInfo, vec);
171
      if (rememberElVec)
172
	userVec += factor * elementVector;
173

174
175
      return;
    } 
176
177

    if (firstOrderAssemblerGrdPsi)
178
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
179
    if (zeroOrderAssembler)
180
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
181
      
182
183
    if (rememberElVec)
      userVec += factor * elementVector;
184
185
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
186
187
188
189
  void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
					 const ElInfo *auxElInfo,
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
190
					 ElementVector& userVec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
191
192
193
194
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

195
    if (remember && factor != 1.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
196
197
198
199
      rememberElVec = true;

    Element *el = mainElInfo->getElement();

200
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
Thomas Witkowski's avatar
Thomas Witkowski committed
201
202
203
      initElement(auxElInfo);
    
    if (el != lastVecEl || !operat->isOptimized()) {
204
205
206
      if (rememberElVec)
	set_to_zero(elementVector);

Thomas Witkowski's avatar
Thomas Witkowski committed
207
208
209
      lastVecEl = el;
    } else {
      if (rememberElVec) {
210
	userVec += factor * elementVector;
Thomas Witkowski's avatar
Thomas Witkowski committed
211
212
213
	return;
      }
    }
214
    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
215
216
217
218
219
220
221
222
223

    if (operat->uhOld && remember) {

      if (smallElInfo->getLevel() == largeElInfo->getLevel()) {
	matVecAssemble(auxElInfo, vec);
      } else {
	matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec);
      }

224
      if (rememberElVec)
225
226
	userVec += factor * elementVector;

Thomas Witkowski's avatar
Thomas Witkowski committed
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
      return;
    } 

    ERROR_EXIT("Not yet implemented!\n");
//     if (firstOrderAssemblerGrdPsi) {
//       firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
//     }
//     if (zeroOrderAssembler) {
//       zeroOrderAssembler->calculateElementVector(elInfo, vec);
//     }
//     if (rememberElVec) {
//       axpy(factor, *elementVector, *userVec);
//     }      

  }

243
  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec)
244
245
246
247
  {
    FUNCNAME("Assembler::matVecAssemble()");

    Element *el = elInfo->getElement(); 
248
    double *uhOldLoc = new double[nRow];
249

250
251
    operat->uhOld->getLocalVector(el, uhOldLoc);
    
252
253
    if (el != lastMatEl)
      calculateElementMatrix(elInfo, elementMatrix);
254
    
255
    for (int i = 0; i < nRow; i++) {
256
      double val = 0.0;
257
      for (int j = 0; j < nRow; j++) {
258
	val += elementMatrix[i][j] * uhOldLoc[j];
259
      }
260
      vec[i] += val;
261
    }   
262

Thomas Witkowski's avatar
Thomas Witkowski committed
263
264
265
266
267
    delete [] uhOldLoc;
  }

  void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
				 const ElInfo *smallElInfo, const ElInfo *largeElInfo,
268
				 ElementVector& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
  {
    FUNCNAME("Assembler::matVecAssemble()");

    TEST_EXIT(rowFESpace->getBasisFcts() == colFESpace->getBasisFcts())
      ("Works only for equal basis functions for different components!\n");

    TEST_EXIT(operat->uhOld->getFESpace()->getMesh() == auxElInfo->getMesh())
      ("Da stimmt was nicht!\n");

    Element *mainEl = mainElInfo->getElement(); 
    Element *auxEl = auxElInfo->getElement();

    const BasisFunction *basFcts = rowFESpace->getBasisFcts();
    int nBasFcts = basFcts->getNumber();
    double *uhOldLoc = new double[nBasFcts];
    double *uhOldLoc2 = new double[nBasFcts];

    operat->uhOld->getLocalVector(auxEl, uhOldLoc);

288
289
    mtl::dense2D<double> m(nRow, nRow);
    smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
290
291
292
293

    for (int i = 0; i < nBasFcts; i++) {
      uhOldLoc2[i] = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
294
	uhOldLoc2[i] += m[j][i] * uhOldLoc[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
295
296
297
298
299
      }      
    }

    
    if (mainEl != lastMatEl) {
300
301
      calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, 
			     elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
304
305
306
    }
    
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
307
	val += elementMatrix[i][j] * uhOldLoc2[j];
Thomas Witkowski's avatar
Thomas Witkowski committed
308
      }
309
      vec[i] += val;
Thomas Witkowski's avatar
Thomas Witkowski committed
310
311
312
313
    }
    
    delete [] uhOldLoc;
    delete [] uhOldLoc2;
314
315
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
316
317
318
  void Assembler::initElement(const ElInfo *smallElInfo, 
			      const ElInfo *largeElInfo,
			      Quadrature *quad)
319
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
320
    if (secondOrderAssembler) 
Thomas Witkowski's avatar
Thomas Witkowski committed
321
      secondOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
322
    if (firstOrderAssemblerGrdPsi)
Thomas Witkowski's avatar
Thomas Witkowski committed
323
      firstOrderAssemblerGrdPsi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
324
    if (firstOrderAssemblerGrdPhi)
Thomas Witkowski's avatar
Thomas Witkowski committed
325
      firstOrderAssemblerGrdPhi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
326
    if (zeroOrderAssembler)
Thomas Witkowski's avatar
Thomas Witkowski committed
327
      zeroOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
328
329
330
  }

  void Assembler::checkQuadratures()
Thomas Witkowski's avatar
Thomas Witkowski committed
331
332
  { 
    if (secondOrderAssembler) {
333
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
      if (!secondOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
336
337
338
339
340
	int degree = operat->getQuadratureDegree(2);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	secondOrderAssembler->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
341
    if (firstOrderAssemblerGrdPsi) {
342
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
343
344
      if (!firstOrderAssemblerGrdPsi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
345
346
347
348
349
	int degree = operat->getQuadratureDegree(1, GRD_PSI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPsi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
350
    if (firstOrderAssemblerGrdPhi) {
351
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
352
353
      if (!firstOrderAssemblerGrdPhi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
354
355
356
357
358
	int degree = operat->getQuadratureDegree(1, GRD_PHI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPhi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
359
    if (zeroOrderAssembler) {
360
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
      if (!zeroOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
363
364
365
366
367
368
369
	int degree = operat->getQuadratureDegree(0);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	zeroOrderAssembler->setQuadrature(quadrature);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
370
371
372
373
374
  void Assembler::finishAssembling()
  {
    lastVecEl = NULL;
    lastMatEl = NULL;
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423

  OptimizedAssembler::OptimizedAssembler(Operator  *op,
					 Quadrature *quad2,
					 Quadrature *quad1GrdPsi,
					 Quadrature *quad1GrdPhi,
					 Quadrature *quad0,
					 const FiniteElemSpace *row,
					 const FiniteElemSpace *col) 
    : Assembler(op, row, col)
  {
    bool opt = (row == col);

    // create sub assemblers
    secondOrderAssembler = 
      SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
    firstOrderAssemblerGrdPsi = 
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
    firstOrderAssemblerGrdPhi = 
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
    zeroOrderAssembler = 
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, opt);

    checkQuadratures();
  }

  StandardAssembler::StandardAssembler(Operator *op,
				       Quadrature *quad2,
				       Quadrature *quad1GrdPsi,
				       Quadrature *quad1GrdPhi,
				       Quadrature *quad0,
				       const FiniteElemSpace *row,
				       const FiniteElemSpace *col) 
    : Assembler(op, row, col)
  {
    remember = false;

    // create sub assemblers
    secondOrderAssembler = 
      SecondOrderAssembler::getSubAssembler(op, this, quad2, false);
    firstOrderAssemblerGrdPsi = 
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, false);
    firstOrderAssemblerGrdPhi = 
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, false);
    zeroOrderAssembler = 
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);

    checkQuadratures();
  }

424
}