Assembler.cc 11.7 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()");

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

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

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

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

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

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

74
75
    if (rememberElMat)
      userMat += factor * elementMatrix;
76
77
  }

78
79
  void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
					 const ElInfo *colElInfo,
80
81
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
82
					 ElementMatrix& userMat,
83
84
85
86
87
88
89
90
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

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

    if (el != lastMatEl || !operat->isOptimized()) {
99
100
101
      if (rememberElMat)
	set_to_zero(elementMatrix);

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

112
113
    TEST_EXIT(zeroOrderAssembler)("Not yet implemented for not zero order assembler!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
114
    if (secondOrderAssembler)
115
      secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
Thomas Witkowski's avatar
Thomas Witkowski committed
116
    if (firstOrderAssemblerGrdPsi)
117
      firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
Thomas Witkowski's avatar
Thomas Witkowski committed
118
    if (firstOrderAssemblerGrdPhi)
119
120
121
122
123
124
125
126
127
128
      firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
    if (zeroOrderAssembler)
      zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);

    const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
    ElementMatrix tmpMat(nRow, nRow);
    tmpMat = m * mat;
    mat = tmpMat;


129
130
    if (rememberElMat)
      userMat += factor * elementMatrix;
131
132
  }

133
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
134
					 ElementVector& userVec,
135
136
137
138
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

Thomas Witkowski's avatar
Thomas Witkowski committed
139
    if (remember && factor != 1.0) {
140
141
142
143
144
      rememberElVec = true;
    }

    Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
145
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
146
147
      initElement(elInfo);
    }
148
    
Thomas Witkowski's avatar
Thomas Witkowski committed
149
    if (el != lastVecEl || !operat->isOptimized()) {
150
151
152
      if (rememberElVec)
	set_to_zero(elementVector);
	
153
154
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
155
      if (rememberElVec) {
156
	userVec += factor * elementVector;
157
158
159
	return;
      }
    }
160
161

    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
162
    if (operat->uhOld && remember) {
163
      matVecAssemble(elInfo, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
164
      if (rememberElVec) {
165
	userVec += factor * elementVector;
166
167
168
      }
      return;
    } 
169
170

    if (firstOrderAssemblerGrdPsi)
171
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
172
    if (zeroOrderAssembler)
173
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
174
175
    if (rememberElVec)
      userVec += factor * elementVector;
176
177
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
180
181
  void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
					 const ElInfo *auxElInfo,
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
182
					 ElementVector& userVec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

    if (remember && factor != 1.0) {
      rememberElVec = true;
    }

    Element *el = mainElInfo->getElement();

    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
      initElement(auxElInfo);
    }
    
    if (el != lastVecEl || !operat->isOptimized()) {
198
199
200
      if (rememberElVec)
	set_to_zero(elementVector);

Thomas Witkowski's avatar
Thomas Witkowski committed
201
202
203
      lastVecEl = el;
    } else {
      if (rememberElVec) {
204
	userVec += factor * elementVector;
Thomas Witkowski's avatar
Thomas Witkowski committed
205
206
207
	return;
      }
    }
208
    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
209
210
211
212
213
214
215
216
217
218

    if (operat->uhOld && remember) {

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

      if (rememberElVec) {
219
	userVec += factor * elementVector;
Thomas Witkowski's avatar
Thomas Witkowski committed
220
      }
221

Thomas Witkowski's avatar
Thomas Witkowski committed
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
      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);
//     }      

  }

238
  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec)
239
240
241
242
243
  {
    FUNCNAME("Assembler::matVecAssemble()");

    Element *el = elInfo->getElement(); 
    const BasisFunction *basFcts = rowFESpace->getBasisFcts();
244
245
    int nBasFcts = basFcts->getNumber();
    double *uhOldLoc = new double[nBasFcts];
246

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

Thomas Witkowski's avatar
Thomas Witkowski committed
261
262
263
264
265
    delete [] uhOldLoc;
  }

  void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
				 const ElInfo *smallElInfo, const ElInfo *largeElInfo,
266
				 ElementVector& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
  {
    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);

286
    const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
Thomas Witkowski's avatar
Thomas Witkowski committed
287
288
289
290

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

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

    delete [] uhOldLoc;
    delete [] uhOldLoc2;
312
313
  }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
368
369
370
371
372
  void Assembler::finishAssembling()
  {
    lastVecEl = NULL;
    lastMatEl = NULL;
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
373
374
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

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

422
}