SecondOrderAssembler.cc 8.78 KB
Newer Older
1
#include <vector>
2
#include <boost/numeric/mtl/mtl.hpp>
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#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)
  {}

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

    // check if a assembler is needed at all
33
    if (op->secondOrder[myRank].size() == 0)
34
35
36
37
38
      return NULL;

    SecondOrderAssembler *newAssembler;

    std::vector<SubAssembler*> *subAssemblers =
39
40
41
      optimized ?
      &optimizedSubAssemblers :
      &standardSubAssemblers;
42

43
    std::vector<OperatorTerm*> opTerms = op->zeroOrder[myRank];
44
45
46
47

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

    // check if a new assembler is needed
48
    for (int i = 0; i < static_cast<int>(subAssemblers->size()); i++) {
49
50
51
52
      std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
    
      sort(assTerms.begin(), assTerms.end());

53
      if (opTerms == assTerms && (*subAssemblers)[i]->getQuadrature() == quad)	
54
55
56
57
58
	return dynamic_cast<SecondOrderAssembler*>((*subAssemblers)[i]);
    }

    // check if all terms are pw_const
    bool pwConst = true;
59
    for (int i = 0; i < static_cast<int>(op->secondOrder[myRank].size()); i++) {
60
61
62
63
64
65
66
      if (!op->secondOrder[myRank][i]->isPWConst()) {
	pwConst = false;
	break;
      }
    }  

    // create new assembler
67
    if (!optimized) {
68
      newAssembler = new Stand2(op, assembler, quad);
69
    } else {
70
      if (pwConst) {
71
	newAssembler = new Pre2(op, assembler, quad);
72
      } else {
73
	newAssembler = new Quad2(op, assembler, quad);
74
      }
75
    }
76
77
78
79
80
81
82
83
84
85
86

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

  Pre2::Pre2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
  {
    q11 = Q11PsiPhi::provideQ11PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
				      owner->getColFESpace()->getBasisFcts(), 
				      quadrature);
87
88
    tmpLALt.resize(omp_get_overall_max_threads());
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
89
90
      tmpLALt[i] = new DimMat<double>*;
      *(tmpLALt[i]) = new DimMat<double>(dim, NO_INIT);
91
92
93
94
95
96
    }
  }

  Pre2::~Pre2()
  {
    for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
97
98
      delete *(tmpLALt[i]);
      delete tmpLALt[i];
99
100
101
    }
  }

102
  void Pre2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
103
104
105
106
107
108
  {
    const int **nEntries;
    const int *k, *l;
    const double *values;

    int myRank = omp_get_thread_num();
Thomas Witkowski's avatar
Thomas Witkowski committed
109
    DimMat<double> **LALt = tmpLALt[0];
110
111
112
    DimMat<double> &tmpMat = *LALt[0];
    tmpMat.set(0.0);

113
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
114
115
116
117
118
119
120
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, 1, LALt);
    }

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

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

123
124
125
126
127
128
129
130
      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]];
	}
131
	mat[i][i] += val;
132
133
134
135
136
137
138
139
140

	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]];
	  }
141
142
	  mat[i][j] += val;
	  mat[j][i] += val;
143
144
145
146
147
148
149
150
151
152
153
154
	}
      }
    } 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]];
	  }
155
	  mat[i][j] += val;
156
157
158
159
160
161
162
163
	}
      }
    }
  }

  Quad2::Quad2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
  {
164
    tmpLALt.resize(omp_get_overall_max_threads());
165
166
167
168
169
170
171
172
  }

  Quad2::~Quad2()
  {
    if (!firstCall) {
      int nPoints = quadrature->getNumPoints();
      for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
	for (int j = 0; j < nPoints; j++) {
173
	  delete tmpLALt[i][j];
174
	}
175
	delete [] tmpLALt[i];
176
177
178
179
      }
    }
  }

180
  void Quad2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
181
182
183
184
185
  {
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();

    if (firstCall) {
186
      tmpLALt[myRank] = new DimMat<double>*[nPoints];
187
188
      for (int j = 0; j < nPoints; j++)
	tmpLALt[myRank][j] = new DimMat<double>(dim, NO_INIT);      
189
190
191
192
193
194
195
196
197
198
    
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
      firstCall = false;
    }

    DimMat<double> **LALt = tmpLALt[myRank];
    
199
    for (int i = 0; i < nPoints; i++)
200
      LALt[i]->set(0.0);
201
202

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

205
206
207
    VectorOfFixVecs< DimVec<double> > *grdPsi, *grdPhi;

    if (symmetric) {
208
209
      // === Symmetric assembling. ===

210
211
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");

212
213
      for (int iq = 0; iq < nPoints; iq++) {
	(*LALt[iq]) *= elInfo->getDet();
Thomas Witkowski's avatar
Thomas Witkowski committed
214
	
215
216
217
218
	grdPsi = psiFast->getGradient(iq);
	grdPhi = phiFast->getGradient(iq);

	for (int i = 0; i < nRow; i++) {
219
	  mat[i][i] += quadrature->getWeight(iq) * 
220
221
222
223
224
	    xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[i]);

	  for (int j = i + 1; j < nCol; j++) {
	    double val = quadrature->getWeight(iq) * 
	      xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
225
226
	    mat[i][j] += val;
	    mat[j][i] += val;
227
228
229
	  }
	}
      }
230
231
    } else {      
      // === Non symmetric assembling. ===
232
233
234
235
236
237
238
239
      for (int iq = 0; iq < nPoints; iq++) {
	(*LALt[iq]) *= elInfo->getDet();

	grdPsi = psiFast->getGradient(iq);
	grdPhi = phiFast->getGradient(iq);

	for (int i = 0; i < nRow; i++) {
	  for (int j = 0; j < nCol; j++) {
240
	    mat[i][j] += quadrature->getWeight(iq) * 
241
242
243
244
245
246
247
248
249
250
251
	      xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
	  }
	}
      }
    }
  }

  Stand2::Stand2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, false)
  {}

252
  void Stand2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
253
254
255
256
257
258
259
260
261
  {
    DimVec<double> grdPsi(dim, NO_INIT);
    VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);

    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();

    int nPoints = quadrature->getNumPoints();

262
    DimMat<double> **LALt = new DimMat<double>*[nPoints];
263
    for (int iq = 0; iq < nPoints; iq++) {
264
      LALt[iq] = new DimMat<double>(dim, NO_INIT);
265
266
267
268
269
270
271
272
273
274
      LALt[iq]->set(0.0);
    }

    int myRank = omp_get_thread_num();

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

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

277
278
279
      for (int iq = 0; iq < nPoints; iq++) {
	(*LALt[iq]) *= elInfo->getDet();

280
	for (int i = 0; i < nCol; i++)
281
282
283
284
	  (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);

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

286
	  mat[i][i] += quadrature->getWeight(iq) * 
287
288
289
	    (grdPsi * ((*LALt[iq]) * grdPhi[i]));

	  for (int j = i + 1; j < nCol; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
290
291
	    double val = quadrature->getWeight(iq) * 
	      (grdPsi * ((*LALt[iq]) * grdPhi[j]));
292
293
	    mat[i][j] += val;
	    mat[j][i] += val;
294
295
296
297
298
299
300
	  }
	}
      }
    } else {      /*  non symmetric assembling   */
      for (int iq = 0; iq < nPoints; iq++) {
	(*LALt[iq]) *= elInfo->getDet();

301
	for (int i = 0; i < nCol; i++)
302
303
304
305
306
	  (*(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++) {
307
	    mat[i][j] += quadrature->getWeight(iq) *
308
309
310
311
312
313
314
	      (grdPsi * ((*LALt[iq]) * grdPhi[j]));
	  }
	}
      }
    }

    for (int iq = 0; iq < nPoints; iq++) 
315
      delete LALt[iq];
316

317
    delete [] LALt;
318
319
320
  }

}