SecondOrderAssembler.cc 9.17 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include <vector>
14
#include <boost/numeric/mtl/mtl.hpp>
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#include "Assembler.h"
#include "SecondOrderAssembler.h"
#include "Operator.h"
#include "QPsiPhi.h"
#include "FiniteElemSpace.h"
#include "Quadrature.h"
#include "DOFVector.h"
#include "OpenMP.h"

namespace AMDiS {

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

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

36
37
38
39
40

  SecondOrderAssembler* SecondOrderAssembler::getSubAssembler(Operator* op,
							      Assembler *assembler,
							      Quadrature *quad,
							      bool optimized)
41
42
43
44
  {
    int myRank = omp_get_thread_num();

    // check if a assembler is needed at all
45
    if (op->secondOrder[myRank].size() == 0)
46
47
48
49
50
      return NULL;

    SecondOrderAssembler *newAssembler;

    std::vector<SubAssembler*> *subAssemblers =
51
52
53
      optimized ?
      &optimizedSubAssemblers :
      &standardSubAssemblers;
54

55
    std::vector<OperatorTerm*> opTerms = op->zeroOrder[myRank];
56
57
58
59

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

    // check if a new assembler is needed
60
    for (unsigned int i = 0; i < subAssemblers->size(); i++) {
61
62
63
64
      std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
    
      sort(assTerms.begin(), assTerms.end());

65
      if (opTerms == assTerms && (*subAssemblers)[i]->getQuadrature() == quad)	
66
67
68
69
70
	return dynamic_cast<SecondOrderAssembler*>((*subAssemblers)[i]);
    }

    // check if all terms are pw_const
    bool pwConst = true;
71
    for (unsigned int i = 0; i < op->secondOrder[myRank].size(); i++) {
72
73
74
75
76
77
78
      if (!op->secondOrder[myRank][i]->isPWConst()) {
	pwConst = false;
	break;
      }
    }  

    // create new assembler
79
    if (!optimized) {
80
      newAssembler = new Stand2(op, assembler, quad);
81
    } else {
82
      if (pwConst) {
83
	newAssembler = new Pre2(op, assembler, quad);
84
      } else {
85
	newAssembler = new Quad2(op, assembler, quad);
86
      }
87
    }
88
89
90
91
92

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

93

94
95
96
  Pre2::Pre2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
  {
97
98
    name = "precalculated second order assembler";

99
100
    q11 = Q11PsiPhi::provideQ11PsiPhi(rowFeSpace->getBasisFcts(), 
				      colFeSpace->getBasisFcts(), 
101
				      quadrature);
102
103
    tmpLALt.resize(omp_get_overall_max_threads());
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
104
105
      tmpLALt[i] = new DimMat<double>*;
      *(tmpLALt[i]) = new DimMat<double>(dim, NO_INIT);
106
107
108
    }
  }

109

110
111
  Pre2::~Pre2()
  {
112
    for (unsigned int i = 0; i < tmpLALt.size(); i++) {
113
114
      delete *(tmpLALt[i]);
      delete tmpLALt[i];
115
116
117
    }
  }

118

119
  void Pre2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
120
121
122
123
124
125
  {
    const int **nEntries;
    const int *k, *l;
    const double *values;

    int myRank = omp_get_thread_num();
Thomas Witkowski's avatar
Thomas Witkowski committed
126
    DimMat<double> **LALt = tmpLALt[0];
127
128
129
    DimMat<double> &tmpMat = *LALt[0];
    tmpMat.set(0.0);

130
131
    for (unsigned int i = 0; i < terms[myRank].size(); i++)
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, 1, LALt);    
132
133
134
135
136

    tmpMat *= elInfo->getDet();
    nEntries = q11->getNumberEntries();

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

139
140
141
142
143
144
145
146
      for (int i = 0; i < nRow; i++) {
	k = q11->getKVec(i, i);
	l = q11->getLVec(i, i);
	values = q11->getValVec(i, i);
	double val = 0.0;
	for (int m = 0; m < nEntries[i][i]; m++) {
	  val += values[m] * tmpMat[k[m]][l[m]];
	}
147
	mat[i][i] += val;
148
149
150
151
152
153
154
155
156

	for (int j = i + 1; j < nCol; j++) {
	  k = q11->getKVec(i, j);
	  l = q11->getLVec(i, j);
	  values = q11->getValVec(i, j);
	  val = 0.0;
	  for (int m = 0; m < nEntries[i][j]; m++) {
	    val += values[m] * tmpMat[k[m]][l[m]];
	  }
157
158
	  mat[i][j] += val;
	  mat[j][i] += val;
159
160
161
162
163
164
165
166
167
168
169
170
	}
      }
    } else {  /*  A not symmetric or psi != phi        */
      for (int i = 0; i < nRow; i++) {
	for (int j = 0; j < nCol; j++) {
	  k = q11->getKVec(i, j);
	  l = q11->getLVec(i, j);
	  values = q11->getValVec(i, j);
	  double val = 0.0;
	  for (int m = 0; m < nEntries[i][j]; m++) {
	    val += values[m] * tmpMat[k[m]][l[m]];
	  }
171
	  mat[i][j] += val;
172
173
174
175
176
	}
      }
    }
  }

177

178
179
180
  Quad2::Quad2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
  {
181
182
    name = "fast quadrature second order assembler";

183
    tmpVec.resize(omp_get_overall_max_threads());
184
    tmpLALt.resize(omp_get_overall_max_threads());
185
186
  }

187

188
189
190
191
  Quad2::~Quad2()
  {
    if (!firstCall) {
      int nPoints = quadrature->getNumPoints();
192
      for (unsigned int i = 0; i < tmpLALt.size(); i++) {
193
	for (int j = 0; j < nPoints; j++)
194
	  delete tmpLALt[i][j];
195
	
196
	delete [] tmpLALt[i];
197
198
199
200
      }
    }
  }

201

202
  void Quad2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
203
204
205
206
207
  {
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();

    if (firstCall) {
208
      tmpVec[myRank] = new DimVec<double>(dim, NO_INIT);
209
      tmpLALt[myRank] = new DimMat<double>*[nPoints];
210
211
      for (int j = 0; j < nPoints; j++)
	tmpLALt[myRank][j] = new DimMat<double>(dim, NO_INIT);      
212
213

#ifdef _OPENMP
214
#pragma omp critical
215
216
#endif
      {   
217
218
219
220
        psiFast = 
	  updateFastQuadrature(psiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI);
        phiFast = 
	  updateFastQuadrature(phiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI);
221
      }
222
223
224
225
      firstCall = false;
    }

    DimMat<double> **LALt = tmpLALt[myRank];
226
    DimVec<double> &dimVec = *(tmpVec[myRank]);
227
    
228
    for (int i = 0; i < nPoints; i++)
229
      LALt[i]->set(0.0);
230

231
    for (unsigned int i = 0; i < terms[myRank].size(); i++)
232
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
233

234
235
    VectorOfFixVecs< DimVec<double> > *grdPsi, *grdPhi;

236
    if (symmetric) {
237
      // === Symmetric assembling. ===
238
239
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");

240
241
      for (int iq = 0; iq < nPoints; iq++) {
	(*LALt[iq]) *= elInfo->getDet();
242
243
	grdPsi = psiFast->getGradient(iq);
	grdPhi = phiFast->getGradient(iq);
244

245
	for (int i = 0; i < nCol; i++) {	  
246
	  amv(quadrature->getWeight(iq), (*LALt[iq]), (*grdPhi)[i], dimVec);
247

248
	  mat[i][i] += (*grdPsi)[i] * dimVec;
249
	  for (int j = i + 1; j < nRow; j++) {
250
	    double tmp = (*grdPhi)[j] * dimVec;
251
252
	    mat[i][j] += tmp;
	    mat[j][i] += tmp;
253
254
255
	  }
	}
      }
256
257
    } else {      
      // === Non symmetric assembling. ===
258
259
260
261
262
      for (int iq = 0; iq < nPoints; iq++) {
	(*LALt[iq]) *= elInfo->getDet();
	grdPsi = psiFast->getGradient(iq);
	grdPhi = phiFast->getGradient(iq);

263
264
	for (int i = 0; i < nRow; i++)
	  for (int j = 0; j < nCol; j++)
265
	    mat[i][j] += quadrature->getWeight(iq) * 
266
267
268
269
270
	      xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
      }
    }
  }

271

272
273
  Stand2::Stand2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, false)
274
275
276
277
  {
    name = "standard second order assembler";
  }

278

279
  void Stand2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
280
281
282
283
  {
    DimVec<double> grdPsi(dim, NO_INIT);
    VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);

284
285
    const BasisFunction *psi = rowFeSpace->getBasisFcts();
    const BasisFunction *phi = colFeSpace->getBasisFcts();
286
287
288

    int nPoints = quadrature->getNumPoints();

289
    DimMat<double> **LALt = new DimMat<double>*[nPoints];
290
    for (int iq = 0; iq < nPoints; iq++) {
291
      LALt[iq] = new DimMat<double>(dim, NO_INIT);
292
293
294
295
296
      LALt[iq]->set(0.0);
    }

    int myRank = omp_get_thread_num();

297
298
    for (unsigned int i = 0; i < terms[myRank].size(); i++)
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
299
300

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

303
304
305
      for (int iq = 0; iq < nPoints; iq++) {
	(*LALt[iq]) *= elInfo->getDet();

306
	for (int i = 0; i < nCol; i++)
307
308
309
310
	  (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);

	for (int i = 0; i < nRow; i++) {
	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
311

312
	  mat[i][i] += quadrature->getWeight(iq) * 
313
314
315
	    (grdPsi * ((*LALt[iq]) * grdPhi[i]));

	  for (int j = i + 1; j < nCol; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
316
317
	    double val = quadrature->getWeight(iq) * 
	      (grdPsi * ((*LALt[iq]) * grdPhi[j]));
318
319
	    mat[i][j] += val;
	    mat[j][i] += val;
320
321
322
323
324
325
326
	  }
	}
      }
    } else {      /*  non symmetric assembling   */
      for (int iq = 0; iq < nPoints; iq++) {
	(*LALt[iq]) *= elInfo->getDet();

327
	for (int i = 0; i < nCol; i++)
328
329
330
331
332
	  (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);

	for (int i = 0; i < nRow; i++) {
	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
	  for (int j = 0; j < nCol; j++) {
333
	    mat[i][j] += quadrature->getWeight(iq) *
334
335
336
337
338
339
340
	      (grdPsi * ((*LALt[iq]) * grdPhi[j]));
	  }
	}
      }
    }

    for (int iq = 0; iq < nPoints; iq++) 
341
      delete LALt[iq];
342
    
343
    delete [] LALt;
344
345
346
  }

}