Assembler.cc 12 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
78
  void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
					 const ElInfo *colElInfo,
79
80
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
81
					 ElementMatrix& userMat,
82
83
84
85
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

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

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

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

109
    if (secondOrderAssembler) {
110
      secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
111
112
113

      ElementMatrix m(nRow, nRow);
      smallElInfo->getSubElemCoordsMat_so(m, rowFESpace->getBasisFcts()->getDegree());
114

115
116
      ElementMatrix tmpMat(nRow, nRow);
      tmpMat = m * mat;
117
      mat = tmpMat;      
118
119
120
    }

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

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

130
131
    if (zeroOrderAssembler) {
      zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
132

133
134
      ElementMatrix m(nRow, nRow);
      smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
135

136
      ElementMatrix tmpMat(nRow, nRow);
137
138
139
140
141
142
      
      if (smallElInfo == colElInfo)
	tmpMat = m * mat;
      else 
	tmpMat = mat * trans(m);

143
144
      mat = tmpMat;
    }
145

146
147
    if (rememberElMat)
      userMat += factor * elementMatrix;
148
149
  }

150

151
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
152
					 ElementVector& userVec,
153
154
155
156
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

157
    if (remember && factor != 1.0)
158
159
160
161
      rememberElVec = true;

    Element *el = elInfo->getElement();

162
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
163
      initElement(elInfo);
164
    
Thomas Witkowski's avatar
Thomas Witkowski committed
165
    if (el != lastVecEl || !operat->isOptimized()) {
166
167
168
      if (rememberElVec)
	set_to_zero(elementVector);
	
169
170
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
171
      if (rememberElVec) {
172
	userVec += factor * elementVector;
173
174
175
	return;
      }
    }
176
177

    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
178
    if (operat->uhOld && remember) {
179
      matVecAssemble(elInfo, vec);
180
      if (rememberElVec)
181
	userVec += factor * elementVector;
182

183
184
      return;
    } 
185
186

    if (firstOrderAssemblerGrdPsi)
187
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
188
    if (zeroOrderAssembler)
189
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
190
      
191
192
    if (rememberElVec)
      userVec += factor * elementVector;
193
194
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
195
196
197
198
  void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
					 const ElInfo *auxElInfo,
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
199
					 ElementVector& userVec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
200
201
202
203
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

204
    if (remember && factor != 1.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
205
206
207
208
      rememberElVec = true;

    Element *el = mainElInfo->getElement();

209
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
Thomas Witkowski's avatar
Thomas Witkowski committed
210
211
212
      initElement(auxElInfo);
    
    if (el != lastVecEl || !operat->isOptimized()) {
213
214
215
      if (rememberElVec)
	set_to_zero(elementVector);

Thomas Witkowski's avatar
Thomas Witkowski committed
216
217
218
      lastVecEl = el;
    } else {
      if (rememberElVec) {
219
	userVec += factor * elementVector;
Thomas Witkowski's avatar
Thomas Witkowski committed
220
221
222
	return;
      }
    }
223
    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
224
225
226
227
228
229
230
231
232

    if (operat->uhOld && remember) {

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

233
      if (rememberElVec)
234
235
	userVec += factor * elementVector;

Thomas Witkowski's avatar
Thomas Witkowski committed
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
      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);
//     }      

  }

252
  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec)
253
254
255
256
  {
    FUNCNAME("Assembler::matVecAssemble()");

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

259
260
    operat->uhOld->getLocalVector(el, uhOldLoc);
    
261
262
    if (el != lastMatEl)
      calculateElementMatrix(elInfo, elementMatrix);
263
    
264
    for (int i = 0; i < nRow; i++) {
265
      double val = 0.0;
266
      for (int j = 0; j < nRow; j++) {
267
	val += elementMatrix[i][j] * uhOldLoc[j];
268
      }
269
      vec[i] += val;
270
    }   
271

Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
274
275
276
    delete [] uhOldLoc;
  }

  void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
				 const ElInfo *smallElInfo, const ElInfo *largeElInfo,
277
				 ElementVector& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
  {
    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);

297
298
    mtl::dense2D<double> m(nRow, nRow);
    smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
299
300
301
302

    for (int i = 0; i < nBasFcts; i++) {
      uhOldLoc2[i] = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
303
	uhOldLoc2[i] += m[j][i] * uhOldLoc[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
304
305
306
307
308
      }      
    }

    
    if (mainEl != lastMatEl) {
309
310
      calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, 
			     elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
311
312
313
314
315
    }
    
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
316
	val += elementMatrix[i][j] * uhOldLoc2[j];
Thomas Witkowski's avatar
Thomas Witkowski committed
317
      }
318
      vec[i] += val;
Thomas Witkowski's avatar
Thomas Witkowski committed
319
320
321
322
    }
    
    delete [] uhOldLoc;
    delete [] uhOldLoc2;
323
324
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
327
  void Assembler::initElement(const ElInfo *smallElInfo, 
			      const ElInfo *largeElInfo,
			      Quadrature *quad)
328
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
329
    if (secondOrderAssembler) 
Thomas Witkowski's avatar
Thomas Witkowski committed
330
      secondOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
331
    if (firstOrderAssemblerGrdPsi)
Thomas Witkowski's avatar
Thomas Witkowski committed
332
      firstOrderAssemblerGrdPsi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
333
    if (firstOrderAssemblerGrdPhi)
Thomas Witkowski's avatar
Thomas Witkowski committed
334
      firstOrderAssemblerGrdPhi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
335
    if (zeroOrderAssembler)
Thomas Witkowski's avatar
Thomas Witkowski committed
336
      zeroOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
337
338
339
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
379
380
381
382
383
  void Assembler::finishAssembling()
  {
    lastVecEl = NULL;
    lastMatEl = NULL;
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
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
424
425
426
427
428
429
430
431
432

  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();
  }

433
}