Assembler.cc 46.1 KB
Newer Older
1
2
3
4
5
6
7
#include "Assembler.h"
#include "Operator.h"
#include "Element.h"
#include "QPsiPhi.h"
#include "ElementMatrix.h"
#include "ElementVector.h"
#include "DOFVector.h"
8
#include "OpenMP.h"
9
10
11
12
13
#include <vector>
#include <algorithm>

namespace AMDiS {

14
15
16
17
  std::vector<SubAssembler*> ZeroOrderAssembler::optimizedSubAssemblers;
  std::vector<SubAssembler*> FirstOrderAssembler::optimizedSubAssemblersGrdPhi;
  std::vector<SubAssembler*> FirstOrderAssembler::optimizedSubAssemblersGrdPsi;
  std::vector<SubAssembler*> SecondOrderAssembler::optimizedSubAssemblers;
18
  
19
20
21
22
  std::vector<SubAssembler*> ZeroOrderAssembler::standardSubAssemblers;
  std::vector<SubAssembler*> FirstOrderAssembler::standardSubAssemblersGrdPhi;
  std::vector<SubAssembler*> FirstOrderAssembler::standardSubAssemblersGrdPsi;
  std::vector<SubAssembler*> SecondOrderAssembler::standardSubAssemblers;
23
24
25
26
27
28
29
30
31
32

  SubAssembler::SubAssembler(Operator *op,
			     Assembler *assembler,
			     Quadrature *quadrat,
			     int order, 
			     bool optimized,
			     FirstOrderType type) 
    : nRow(0),
      nCol(0),
      coordsAtQPs(NULL),
Thomas Witkowski's avatar
Thomas Witkowski committed
33
      coordsNumAllocated(0),
34
35
36
37
38
39
40
41
42
43
44
45
46
47
      quadrature(quadrat),
      psiFast(NULL),
      phiFast(NULL),
      owner(assembler),
      symmetric(true),
      opt(optimized),
      firstCall(true)
  {
    const BasisFunction *psi = assembler->rowFESpace->getBasisFcts();
    const BasisFunction *phi = assembler->colFESpace->getBasisFcts();

    nRow = psi->getNumber();
    nCol = phi->getNumber();

48
49
50
    int maxThreads = omp_get_max_threads();
    terms.resize(maxThreads);

51
    switch (order) {
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
    case 0:
      terms = op->zeroOrder;
      break;
    case 1:
      if(type == GRD_PHI)
	terms = op->firstOrderGrdPhi;
      else 
	terms = op->firstOrderGrdPsi;
      break;
    case 2:
      terms = op->secondOrder;
      break;
    }

    // check if all terms are symmetric
    symmetric = true;
68
69
    for (int i = 0; i < static_cast<int>(terms[0].size()); i++) {
      if (!(terms[0][i])->isSymmetric()) {
70
71
72
73
74
75
76
77
78
79
80
81
	symmetric = false;
	break;
      }
    }  

    dim = assembler->rowFESpace->getMesh()->getDim();
  }

  FastQuadrature *SubAssembler::updateFastQuadrature(FastQuadrature *quadFast,
						     const BasisFunction *psi,
						     Flag updateFlag)
  {
82
83
84
85
    if (!quadFast) {
      quadFast = FastQuadrature::provideFastQuadrature(psi,
						       *quadrature,
						       updateFlag);
86
    } else {
87
      if (!quadFast->initialized(updateFlag))
88
89
90
91
92
93
94
95
96
97
98
99
100
	quadFast->init(updateFlag);
    }

    return quadFast;
  }

  void SubAssembler::initElement(const ElInfo* elInfo,
				 Quadrature *quad)
  {
    // set corrdsAtQPs invalid
    coordsValid = false;

    // set values at QPs invalid
101
    std::map<const DOFVectorBase<double>*, ValuesAtQPs*>::iterator it1;
102
    for (it1 = valuesAtQPs.begin(); it1 != valuesAtQPs.end(); ++it1) {
103
104
105
106
      ((*it1).second)->valid = false;
    }

    // set gradients at QPs invalid
107
    std::map<const DOFVectorBase<double>*, GradientsAtQPs*>::iterator it2;
108
    for (it2 = gradientsAtQPs.begin(); it2 != gradientsAtQPs.end(); ++it2) {
109
110
111
      ((*it2).second)->valid = false;
    }

112
    int myRank = omp_get_thread_num();
113
    // calls initElement of each term
114
    std::vector<OperatorTerm*>::iterator it;
115
    for (it = terms[myRank].begin(); it != terms[myRank].end(); ++it) {
116
117
118
119
120
121
122
123
124
      (*it)->initElement(elInfo, this, quad);
    }
  }

  WorldVector<double>* SubAssembler::getCoordsAtQPs(const ElInfo* elInfo,
						    Quadrature *quad)
  {
    Quadrature *localQuad = quad ? quad : quadrature;
  
125
    const int nPoints = localQuad->getNumPoints();
126
127

    // already calculated for this element ?
128
    if (coordsValid) {
129
130
131
      return coordsAtQPs;
    }
   
Thomas Witkowski's avatar
Thomas Witkowski committed
132
    if (coordsAtQPs)  {
133
      if (coordsNumAllocated != nPoints) {
Thomas Witkowski's avatar
Thomas Witkowski committed
134
	DELETE [] coordsAtQPs;
135
136
        coordsAtQPs = NEW WorldVector<double>[nPoints];
	coordsNumAllocated = nPoints;
Thomas Witkowski's avatar
Thomas Witkowski committed
137
138
      }
    } else {
139
140
      coordsAtQPs = NEW WorldVector<double>[nPoints];
      coordsNumAllocated = nPoints;
Thomas Witkowski's avatar
Thomas Witkowski committed
141
    }
142
143

    // set new values
Thomas Witkowski's avatar
Thomas Witkowski committed
144
    WorldVector<double>* k = &(coordsAtQPs[0]);
145
    for (int l = 0; k < &(coordsAtQPs[nPoints]); ++k, ++l) {
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
      elInfo->coordToWorld(localQuad->getLambda(l), k);
    }

    // mark values as valid
    coordsValid = true;

    return coordsAtQPs;
  }

  double* SubAssembler::getVectorAtQPs(DOFVectorBase<double>* dv, 
				       const ElInfo* elInfo,
				       Quadrature *quad)
  {
    FUNCNAME("SubAssembler::getVectorAtQPs()");

    const DOFVectorBase<double>* vec = dv ? dv : owner->operat->getUhOld();

163
    TEST_EXIT_DBG(vec)("no dof vector!\n");
164

165
    if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid) 
166
167
168
169
      return valuesAtQPs[vec]->values.getValArray();

    Quadrature *localQuad = quad ? quad : quadrature;

170
    if (!valuesAtQPs[vec]) {
171
172
173
174
175
176
177
178
179
180
      valuesAtQPs[vec] = new ValuesAtQPs;
    }
    valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());

    double *values = valuesAtQPs[vec]->values.getValArray();

    bool sameFESpaces = 
      (vec->getFESpace() == owner->rowFESpace) || 
      (vec->getFESpace() == owner->colFESpace);

181
    if (opt && !quad && sameFESpaces) {
182
183
      const BasisFunction *psi = owner->rowFESpace->getBasisFcts();
      const BasisFunction *phi = owner->colFESpace->getBasisFcts();
184
      if (vec->getFESpace()->getBasisFcts() == psi) {
185
186
187
188
189
190
191
192
193
	psiFast = updateFastQuadrature(psiFast, psi, INIT_PHI);
      } else if(vec->getFESpace()->getBasisFcts() == phi) {
	phiFast = updateFastQuadrature(phiFast, phi, INIT_PHI);
      }
    }

    // calculate new values
    const BasisFunction *basFcts = vec->getFESpace()->getBasisFcts();

194
195
    if (opt && !quad && sameFESpaces) {
      if (psiFast->getBasisFunctions() == basFcts) {
196
197
198
199
200
201
202
203
204
	vec->getVecAtQPs(elInfo, NULL, psiFast, values);
      } else if(phiFast->getBasisFunctions() == basFcts) {
	vec->getVecAtQPs(elInfo, NULL, phiFast, values);
      } else {
	vec->getVecAtQPs(elInfo, localQuad, NULL, values);
      }
    } else {
      vec->getVecAtQPs(elInfo, localQuad, NULL, values);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
205
  
206
207
208
209
210
211
212
213
214
215
216
217
218
    valuesAtQPs[vec]->valid = true;

    return values;
  }

  WorldVector<double>* SubAssembler::getGradientsAtQPs(DOFVectorBase<double>* dv, 
						       const ElInfo* elInfo,
						       Quadrature *quad)
  {
    FUNCNAME("SubAssembler::getGradientsAtQPs()");

    const DOFVectorBase<double>* vec = dv ? dv : owner->operat->getUhOld();

219
    TEST_EXIT_DBG(vec)("no dof vector!\n");
220

221
    if (gradientsAtQPs[vec] && gradientsAtQPs[vec]->valid) 
222
223
224
225
      return gradientsAtQPs[vec]->values.getValArray();

    Quadrature *localQuad = quad ? quad : quadrature;

226
    if (!gradientsAtQPs[vec]) {
227
228
229
230
231
232
233
234
235
236
237
238
239
      gradientsAtQPs[vec] = new GradientsAtQPs;
    } 
    gradientsAtQPs[vec]->values.resize(localQuad->getNumPoints());

    WorldVector<double> *values = gradientsAtQPs[vec]->values.getValArray();

    const BasisFunction *psi = owner->rowFESpace->getBasisFcts();
    const BasisFunction *phi = owner->colFESpace->getBasisFcts();

    bool sameFESpaces = 
      (vec->getFESpace() == owner->rowFESpace) || 
      (vec->getFESpace() == owner->colFESpace);

240
241
    if (opt && !quad && sameFESpaces) {
      if (vec->getFESpace()->getBasisFcts() == psi) {
242
243
244
245
246
247
248
249
250
	psiFast = updateFastQuadrature(psiFast, psi, INIT_GRD_PHI);
      } else if(vec->getFESpace()->getBasisFcts() == phi) {
	phiFast = updateFastQuadrature(phiFast, phi, INIT_GRD_PHI);
      }
    }
  
    // calculate new values
    const BasisFunction *basFcts = vec->getFESpace()->getBasisFcts();

251
    if (opt && !quad && sameFESpaces) {
252
      if (psiFast->getBasisFunctions() == basFcts) {
253
254
255
256
257
258
259
	vec->getGrdAtQPs(elInfo, NULL, psiFast, values);
      } else {
	vec->getGrdAtQPs(elInfo, NULL, phiFast, values);
      }
    } else {
      vec->getGrdAtQPs(elInfo, localQuad, NULL, values);
    }
260
   
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
    gradientsAtQPs[vec]->valid = true;

    return values;
  }

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

  FirstOrderAssembler::FirstOrderAssembler(Operator *op,
					   Assembler *assembler,
					   Quadrature *quad,
					   bool optimized,
					   FirstOrderType type)
    : SubAssembler(op, assembler, quad, 1, optimized, type)
  {}

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

  ZeroOrderAssembler* 
  ZeroOrderAssembler::getSubAssembler(Operator* op,
				      Assembler *assembler,
				      Quadrature *quad,
				      bool optimized)
  {
    // check if a assembler is needed at all
295
    if (op->zeroOrder.size() == 0) {
296
297
298
299
300
      return NULL;
    }

    ZeroOrderAssembler *newAssembler;

301
    std::vector<SubAssembler*> *subAssemblers =
302
303
304
305
	optimized ?
	&optimizedSubAssemblers :
    &standardSubAssemblers;

306
    int myRank = omp_get_thread_num();
307
    std::vector<OperatorTerm*> opTerms  = op->zeroOrder[myRank];
308
309
310
311

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

    // check if a new assembler is needed
312
313
    if (quad) {
      for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
314
	std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
315
316
317

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

318
319
320
321
322
	if ((opTerms == assTerms) && 
	    ((*subAssemblers)[i]->getQuadrature() == quad)) {
	
	  return dynamic_cast<ZeroOrderAssembler*>((*subAssemblers)[i]);
	}
323
324
      }
    }
325
 
326
327
    // check if all terms are pw_const
    bool pwConst = true;
328
329
    for (int i = 0; i < static_cast<int>( op->zeroOrder[myRank].size()); i++) {
      if (!op->zeroOrder[myRank][i]->isPWConst()) {
330
331
332
333
334
335
	pwConst = false;
	break;
      }
    }  

    // create new assembler
336
    if (!optimized) {
337
338
      newAssembler = NEW Stand0(op, assembler, quad);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
339
      if (pwConst) {
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
	newAssembler = NEW Pre0(op, assembler, quad);
      } else {
	newAssembler = NEW Quad0(op, assembler, quad);
      }
    }

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

  FirstOrderAssembler* 
  FirstOrderAssembler::getSubAssembler(Operator* op,
				       Assembler *assembler,
				       Quadrature *quad,
				       FirstOrderType type,
				       bool optimized)
  {
357
    std::vector<SubAssembler*> *subAssemblers =
358
359
360
361
362
363
364
365
	optimized ?
	(type == GRD_PSI ? 
	 &optimizedSubAssemblersGrdPsi : 
	 &optimizedSubAssemblersGrdPhi) :
    (type == GRD_PSI ? 
     &standardSubAssemblersGrdPsi :
     &standardSubAssemblersGrdPhi);

366
    int myRank = omp_get_thread_num();
367
    std::vector<OperatorTerm*> opTerms 
368
369
370
	= (type == GRD_PSI) ? 
	    op->firstOrderGrdPsi[myRank] : 
            op->firstOrderGrdPhi[myRank];
371
372

    // check if a assembler is needed at all
373
    if (opTerms.size() == 0) {
374
375
376
377
378
379
380
381
      return NULL;
    }

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

    FirstOrderAssembler *newAssembler;

    // check if a new assembler is needed
382
    for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
383

384
      std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
385
    
386
      sort(assTerms.begin(), assTerms.end());   
387

388
389
      if ((opTerms == assTerms) && 
	  ((*subAssemblers)[i]->getQuadrature() == quad)) {
390

391
392
	return dynamic_cast<FirstOrderAssembler*>((*subAssemblers)[i]);
      }
393
394
395
396
    }

    // check if all terms are pw_const
    bool pwConst = true;
397
398
    for (int i = 0; i < static_cast<int>( opTerms.size()); i++) {
      if (!(opTerms[i])->isPWConst()) {
399
400
401
402
403
404
	pwConst = false;
	break;
      }
    }  

    // create new assembler
405
    if (!optimized) {
406
407
408
409
410
      newAssembler = 
	(type == GRD_PSI) ?
	dynamic_cast<FirstOrderAssembler*>(NEW Stand10(op, assembler, quad)) :
	dynamic_cast<FirstOrderAssembler*>(NEW Stand01(op, assembler, quad));    
    } else {
411
      if (pwConst) {
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
	newAssembler = 
	  (type == GRD_PSI) ?
	  dynamic_cast<FirstOrderAssembler*>(NEW Pre10(op, assembler, quad)) :
	  dynamic_cast<FirstOrderAssembler*>(NEW Pre01(op, assembler, quad));
      } else {
	newAssembler = 
	  (type == GRD_PSI) ?
	  dynamic_cast<FirstOrderAssembler*>( NEW Quad10(op, assembler, quad)) :
	  dynamic_cast<FirstOrderAssembler*>( NEW Quad01(op, assembler, quad));
      }
    }

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

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

436
    // check if a assembler is needed at all
437
    if (op->secondOrder[myRank].size() == 0) {
438
439
440
441
442
      return NULL;
    }

    SecondOrderAssembler *newAssembler;

443
    std::vector<SubAssembler*> *subAssemblers =
444
445
446
447
	optimized ?
	&optimizedSubAssemblers :
    &standardSubAssemblers;

448
    std::vector<OperatorTerm*> opTerms  = op->zeroOrder[myRank];
449
450
451
452

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

    // check if a new assembler is needed
453
    for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
454
      std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
455
456
457
    
      sort(assTerms.begin(), assTerms.end());

458
459
460
461
462
      if ((opTerms == assTerms) && 
	  ((*subAssemblers)[i]->getQuadrature() == quad)) {
	
	return dynamic_cast<SecondOrderAssembler*>((*subAssemblers)[i]);
      }
463
464
465
466
    }

    // check if all terms are pw_const
    bool pwConst = true;
467
468
    for (int i = 0; i < static_cast<int>( op->secondOrder[myRank].size()); i++) {
      if (!op->secondOrder[myRank][i]->isPWConst()) {
469
470
471
472
473
474
	pwConst = false;
	break;
      }
    }  

    // create new assembler
475
    if (!optimized) {
476
477
      newAssembler = NEW Stand2(op, assembler, quad);
    } else {
478
      if (pwConst) {
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
	newAssembler = NEW Pre2(op, assembler, quad);
      } else {
	newAssembler = NEW Quad2(op, assembler, quad);
      }
    }

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

  Stand0::Stand0(Operator *op, Assembler *assembler, Quadrature *quad)
    : ZeroOrderAssembler(op, assembler, quad, false)
  {
  }

  void Stand0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();

    double *phival = GET_MEMORY(double, nCol);
500
501
    int nPoints = quadrature->getNumPoints();
    double *c = GET_MEMORY(double, nPoints);
Thomas Witkowski's avatar
Thomas Witkowski committed
502

503
    for (int iq = 0; iq < nPoints; iq++) {
504
505
506
      c[iq] = 0.0;
    }

507
    int myRank = omp_get_thread_num();
508
    std::vector<OperatorTerm*>::iterator termIt;
509
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
510
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
511
512
513
    }
      
    if (symmetric) {
514
      for (int iq = 0; iq < nPoints; iq++) {
515
516
517
	c[iq] *= elInfo->getDet();

	// calculate phi at QPs only once!
518
	for (int i = 0; i < nCol; i++) {
519
520
521
	  phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
	}

522
	for (int i = 0; i < nRow; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
523
	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
Thomas Witkowski's avatar
Thomas Witkowski committed
524
	  (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psival * phival[i];
525
	  for (int j = i + 1; j < nCol; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
526
	    double val = quadrature->getWeight(iq) * c[iq] * psival * phival[j];
527
528
529
530
531
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
532
    } else {      //  non symmetric assembling 
533
      for (int iq = 0; iq < nPoints; iq++) {
534
535
536
	c[iq] *= elInfo->getDet();

	// calculate phi at QPs only once!
537
	for (int i = 0; i < nCol; i++) {
538
539
540
	  phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
	}

541
	for (int i = 0; i < nRow; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
542
	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
543
	  for (int j = 0; j < nCol; j++) {
544
	    (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j];
545
546
547
548
	  }
	}
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
549

550
    FREE_MEMORY(phival, double, nCol);
551
    FREE_MEMORY(c, double, nPoints);
552
553
554
555
  }

  void Stand0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
556
    int nPoints = quadrature->getNumPoints();
557

558
    double *c = GET_MEMORY(double, nPoints);
559

560
    for (int iq = 0; iq < nPoints; iq++) {
561
562
563
      c[iq] = 0.0;
    }

564
    int myRank = omp_get_thread_num();
565
    std::vector<OperatorTerm*>::iterator termIt;
566
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
567
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
568
569
    }

570
    for (int iq = 0; iq < nPoints; iq++) {
571
572
      c[iq] *= elInfo->getDet();

573
574
      for (int i = 0; i < nRow; i++) {
	double psi = (*(owner->getRowFESpace()->getBasisFcts()->getPhi(i)))
575
	  (quadrature->getLambda(iq));
576
	(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi;
577
578
      }
    }
579
    
580
    FREE_MEMORY(c, double, nPoints);
581
582
583
584
585
  }

  Quad0::Quad0(Operator *op, Assembler *assembler, Quadrature *quad)
    : ZeroOrderAssembler(op, assembler, quad, true)
  {
586
587
588
589
590
591
592
593
    cPtrs.resize(omp_get_max_threads());
  }

  Quad0::~Quad0()
  {
    for (int i = 0; i < omp_get_max_threads(); i++) {
      FREE_MEMORY(cPtrs[i], double, quadrature->getNumPoints());
    }
594
595
596
597
  }

  void Quad0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
598
599
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();
600

601
    if (firstCall) {
602
      cPtrs[myRank] = GET_MEMORY(double, nPoints);
603
604
605
606
607
608
609
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;
    }

610
611
    double *c = cPtrs[myRank];
    for (int iq = 0; iq < nPoints; iq++) {
612
613
614
      c[iq] = 0.0;
    }

615
    std::vector<OperatorTerm*>::iterator termIt;
616
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
617
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
618
619
620
    }

    if (symmetric) {
621
      for (int iq = 0; iq < nPoints; iq++) {
622
623
	c[iq] *= elInfo->getDet();

624
625
	const double *psi = psiFast->getPhi(iq);
	const double *phi = phiFast->getPhi(iq);
626
627
628
629
	for (int i = 0; i < nRow; i++) {
	  (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[i];
	  for (int j = i + 1; j < nCol; j++) {
	    double val = quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
630
631
632
633
634
635
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
    } else {      /*  non symmetric assembling   */
636
      for (int iq = 0; iq < nPoints; iq++) {
637
638
	c[iq] *= elInfo->getDet();

639
640
	const double *psi = psiFast->getPhi(iq);
	const double *phi = phiFast->getPhi(iq);
641
642
643
	for (int i = 0; i < nRow; i++) {
	  for (int j = 0; j < nCol; j++) {
	    (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
644
645
646
647
648
649
650
651
	  }
	}
      }
    }
  }

  void Quad0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
652
    if (firstCall) {
653
654
655
656
657
658
659
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;
    }

660
661
    int nPoints = quadrature->getNumPoints();
    double *c = GET_MEMORY(double, nPoints);
662

663
    for (int iq = 0; iq < nPoints; iq++) {
664
665
666
      c[iq] = 0.0;
    }

667
    int myRank = omp_get_thread_num();
668
    std::vector<OperatorTerm*>::iterator termIt;
669
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
670
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
671
672
    }

673
    for (int iq = 0; iq < nPoints; iq++) {
674
675
      c[iq] *= elInfo->getDet();

676
677
678
      const double *psi = psiFast->getPhi(iq);
      for (int i = 0; i < nRow; i++) {
	(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
679
680
      }
    }
681
    FREE_MEMORY(c, double, nPoints);
682
683
684
685
686
687
688
689
690
  }

  Pre0::Pre0(Operator *op, Assembler *assembler, Quadrature *quad) 
    : ZeroOrderAssembler(op, assembler, quad, true)
  {
  }

  void Pre0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
691
    if (firstCall) {
692
693
694
695
696
697
698
699
      q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

700
701
    //    c[0] = 0.0;
    double c = 0.0;
702
    int myRank = omp_get_thread_num();
703
    int size = static_cast<int>(terms[myRank].size());
704

705
706
    for (int i = 0; i < size; i++) {
      (static_cast<ZeroOrderTerm*>((terms[myRank][i])))->getC(elInfo, 1, &c);
707
708
    }

709
    c *= elInfo->getDet();
710
711

    if (symmetric) {
712
      for (int i = 0; i < nRow; i++) {
713
	(*mat)[i][i] += c * q00->getValue(i,i);
714
	for (int j = i + 1; j < nCol; j++) {
715
	  double val = c * q00->getValue(i, j);
716
717
718
719
720
	  (*mat)[i][j] += val;
	  (*mat)[j][i] += val;
	}
      }
    } else {
721
722
      for (int i = 0; i < nRow; i++)
	for (int j = 0; j < nCol; j++)
723
	  (*mat)[i][j] += c * q00->getValue(i, j);
724
725
726
727
728
    }
  }

  void Pre0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
729
    if (firstCall) {
730
731
732
733
734
735
736
737
      q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

738
    std::vector<OperatorTerm*>::iterator termIt;
739

740
    int myRank = omp_get_thread_num();
741
    double c = 0.0;
742
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
743
      (static_cast<ZeroOrderTerm*>( *termIt))->getC(elInfo, 1, &c);
744
745
    }

746
    c *= elInfo->getDet();
747

748
    for (int i = 0; i < nRow; i++)
749
      (*vec)[i] += c * q0->getValue(i);
750
751
752
753
754
755
756
757
758
  }

  Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
  {}


  void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
759
    DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
760
761
762
763
764
    double *phival = GET_MEMORY(double, nCol);

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

765
    int nPoints = quadrature->getNumPoints();
766

767
768
    VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
    for (int iq = 0; iq < nPoints; iq++) {
769
770
      Lb[iq].set(0.0);
    }
771
772
773

    int myRank = omp_get_thread_num();
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
774
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
775
776
    }
  
777
    for (int iq = 0; iq < nPoints; iq++) {
778
779
      Lb[iq] *= elInfo->getDet();

780
      for (int i = 0; i < nCol; i++) {
781
782
783
	phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
      }

784
      for (int i = 0; i < nRow; i++) {
785
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
786
787
	for (int j = 0; j < nCol; j++) {
	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j];
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
	}
      }
    }
    FREE_MEMORY(phival, double, nCol);
  }


  Quad10::Quad10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
  {
  }


  void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    VectorOfFixVecs<DimVec<double> > *grdPsi;
    const double *phi;

806
    if (firstCall) {
807
808
809
810
811
812
813
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;
    }

814
    int nPoints = quadrature->getNumPoints();
815

816
817
    VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
    for (int iq = 0; iq < nPoints; iq++) {
818
819
      Lb[iq].set(0.0);
    }
820

821
822
    int myRank = omp_get_thread_num();
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
823
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
824
825
    }
  
826
    for (int iq = 0; iq < nPoints; iq++) {
827
828
829
830
831
      Lb[iq] *= elInfo->getDet();

      phi    = phiFast->getPhi(iq);
      grdPsi = psiFast->getGradient(iq);

832
833
834
      for (int i = 0; i < nRow; i++) {
	for (int j = 0; j < nCol; j++)
	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
835
836
837
838
839
840
841
842
843
844
      }
    }
  }


  Pre10::Pre10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
  {
  }

845

846
847
  void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
848
    VectorOfFixVecs<DimVec<double> > Lb(dim, 1, NO_INIT);
849
850
851
    const int *k;
    const double *values;

852
    if (firstCall) {
853
854
855
856
857
858
859
860
861
      q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    const int **nEntries = q10->getNumberEntries();
862
    int myRank = omp_get_thread_num();
863
    Lb[0].set(0.0);
864
865
866

    for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, 1, Lb);
867
868
869
870
    }

    Lb[0] *= elInfo->getDet();

871
872
873
    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
	k = q10->getKVec(i, j);
874
	values = q10->getValVec(i, j);
875
	double val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
876
877
878
	for (int m = 0; m < nEntries[i][j]; m++) {
	  val += values[m] * Lb[0][k[m]];
	}
879
880
881
882
883
884
885
886
887
888
	(*mat)[i][j] += val;
      }
    }
  }


  Stand01::Stand01(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
  {}

889

890
891
892
893
894
895
896
  void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);

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

897
898
    int nPoints = quadrature->getNumPoints();
    VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
899
    int myRank = omp_get_thread_num();
900

901
    for (int iq = 0; iq < nPoints; iq++) {
902
903
      Lb[iq].set(0.0);
    }
904
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
905
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
906
907
    }
  
908
    for (int iq = 0; iq < nPoints; iq++) {
909
910
      Lb[iq] *= elInfo->getDet();

911
      for (int i = 0; i < nCol; i++) {
912
	(*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
913
914
      }

915
      for (int i = 0; i < nRow; i++) {
916
	double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
917
918
	for (int j = 0; j < nCol; j++)
	  (*mat)[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * psival) * grdPhi[j]);
919
      }
920
    } 
921
922
923
924
  }

  void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
925
    DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
926
    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
927
928
    int nPoints = quadrature->getNumPoints();
    VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
929
    int myRank = omp_get_thread_num();
930

931
    for (int iq = 0; iq < nPoints; iq++) {
932
933
      Lb[iq].set(0.0);
    }
934
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
935
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
936
937
    }
  
938
    for (int iq = 0; iq < nPoints; iq++) {
939
940
      Lb[iq] *= elInfo->getDet();

941
      for (int i = 0; i < nRow; i++) {
942
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
943
944
945
946
947
948
949
950
951
952
953
954
955
956
	(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
      }
    }
  }

  Quad01::Quad01(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
  {
  }

  void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    VectorOfFixVecs<DimVec<double> > *grdPhi;

957
    if (firstCall) {
958
959
960
961
962
963
964
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
      firstCall = false;
    }

965
966
    int nPoints = quadrature->getNumPoints();
    VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
967
    int myRank = omp_get_thread_num();
968

969
    for (int iq = 0; iq < nPoints; iq++) {
970
971
      Lb[iq].set(0.0);
    }
972
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
973
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
974
975
    }
  
976
    for (int iq = 0; iq < nPoints; iq++) {
977
978
      Lb[iq] *= elInfo->getDet();

979
      const double *psi = psiFast->getPhi(iq);
980
981
      grdPhi = phiFast->getGradient(iq);

982
983
984
      for (int i = 0; i < nRow; i++) {
	for (int j = 0; j < nCol; j++)
	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPhi)[j]) * psi[i];
985
986
987
988
989
990
991
992
      }
    }
  }

  void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    VectorOfFixVecs<DimVec<double> > *grdPsi;

993
    if (firstCall) {
994
995
996
997
998
999
1000
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;
    }

1001
1002
    int nPoints = quadrature->getNumPoints();
    VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
1003
    int myRank = omp_get_thread_num();
1004

1005
    for (int iq = 0; iq < nPoints; iq++) {
1006
1007
      Lb[iq].<