Assembler.cc 11.6 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

Thomas Witkowski's avatar
Thomas Witkowski committed
112
113
114
115
116
117
118
119
120
121
122
    if (secondOrderAssembler)
      secondOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo, 
						   smallElInfo, largeElInfo, mat);   
    if (firstOrderAssemblerGrdPsi)
      firstOrderAssemblerGrdPsi->calculateElementMatrix(rowElInfo, colElInfo, 
							smallElInfo, largeElInfo, mat);    
    if (firstOrderAssemblerGrdPhi)
      firstOrderAssemblerGrdPhi->calculateElementMatrix(rowElInfo, colElInfo, 
							smallElInfo, largeElInfo, mat);    
    if (zeroOrderAssembler) 
      zeroOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo, 
123
						 smallElInfo, largeElInfo, mat);
Thomas Witkowski's avatar
Thomas Witkowski committed
124
  
125
126
    if (rememberElMat)
      userMat += factor * elementMatrix;
127
128
  }

129
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
130
					 ElementVector& userVec,
131
132
133
134
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

Thomas Witkowski's avatar
Thomas Witkowski committed
135
    if (remember && factor != 1.0) {
136
137
138
139
140
      rememberElVec = true;
    }

    Element *el = elInfo->getElement();

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

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

    if (firstOrderAssemblerGrdPsi)
167
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
168
    if (zeroOrderAssembler)
169
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
170
171
    if (rememberElVec)
      userVec += factor * elementVector;
172
173
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
174
175
176
177
  void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
					 const ElInfo *auxElInfo,
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
178
					 ElementVector& userVec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
					 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()) {
194
195
196
      if (rememberElVec)
	set_to_zero(elementVector);

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

    if (operat->uhOld && remember) {

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

      if (rememberElVec) {
215
	userVec += factor * elementVector;
Thomas Witkowski's avatar
Thomas Witkowski committed
216
      }
217

Thomas Witkowski's avatar
Thomas Witkowski committed
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
      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);
//     }      

  }

234
  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec)
235
236
237
238
239
  {
    FUNCNAME("Assembler::matVecAssemble()");

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
257
258
259
260
261
    delete [] uhOldLoc;
  }

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

282
    const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
Thomas Witkowski's avatar
Thomas Witkowski committed
283
284
285
286

    for (int i = 0; i < nBasFcts; i++) {
      uhOldLoc2[i] = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
287
	uhOldLoc2[i] += m[j][i] * uhOldLoc[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
288
289
290
291
292
      }      
    }

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

    delete [] uhOldLoc;
    delete [] uhOldLoc2;
308
309
  }

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

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

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

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

418
}