FirstOrderAssembler.cc 14.1 KB
Newer Older
1 2 3 4 5 6 7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9 10 11 12 13 14 15 16 17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20 21


22
#include <vector>
23
#include <boost/numeric/mtl/mtl.hpp>
24 25 26 27 28 29 30 31
#include "Assembler.h"
#include "FirstOrderAssembler.h"
#include "Operator.h"
#include "QPsiPhi.h"
#include "FiniteElemSpace.h"
#include "Quadrature.h"
#include "DOFVector.h"

32 33
using namespace std;

34 35
namespace AMDiS {

36
  ThreadPrivate<vector<SubAssembler*> >
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
37 38 39
  FirstOrderAssembler::optimizedSubAssemblersGrdPhi;
  ThreadPrivate<vector<SubAssembler*> >
  FirstOrderAssembler::optimizedSubAssemblersGrdPsi;
40

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
41 42 43 44
  ThreadPrivate<vector<SubAssembler*> >
  FirstOrderAssembler::standardSubAssemblersGrdPhi;
  ThreadPrivate<vector<SubAssembler*> >
  FirstOrderAssembler::standardSubAssemblersGrdPsi;
45 46 47 48 49 50 51 52


  FirstOrderAssembler::FirstOrderAssembler(Operator *op,
					   Assembler *assembler,
					   Quadrature *quad,
					   bool optimized,
					   FirstOrderType type)
    : SubAssembler(op, assembler, quad, 1, optimized, type)
53
  {
54
    FUNCNAME_DBG("FirstOrderAssmebler::FirstOrderAssembler()");
55

56 57 58 59
    TEST_EXIT_DBG(dim > 0)("Should not happen!\n");

    Lb.resize(1);
    Lb[0].change_dim(dim + 1);
60
  }
61

62 63 64 65 66 67

  FirstOrderAssembler* FirstOrderAssembler::getSubAssembler(Operator* op,
							    Assembler *assembler,
							    Quadrature *quad,
							    FirstOrderType type,
							    bool optimized)
68
  {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
69
    vector<SubAssembler*> &subAssemblers =
70
      optimized ?
71 72
      (type == GRD_PSI ?
       optimizedSubAssemblersGrdPsi.get() :
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
73
       optimizedSubAssemblersGrdPhi.get()) :
74
      (type == GRD_PSI ?
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
75 76
       standardSubAssemblersGrdPsi.get() :
       standardSubAssemblersGrdPhi.get());
77 78

    vector<OperatorTerm*> opTerms =
79
      (type == GRD_PSI) ? op->firstOrderGrdPsi : op->firstOrderGrdPhi;
80 81

    // check if a assembler is needed at all
82
    if (opTerms.size() == 0)
83
      return NULL;
84 85 86 87 88 89

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

    FirstOrderAssembler *newAssembler;

    // check if a new assembler is needed
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
90 91
    for (unsigned int i = 0; i < subAssemblers.size(); i++) {
      vector<OperatorTerm*> assTerms = *(subAssemblers[i]->getTerms());
92 93

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
95 96
      if (opTerms == assTerms && subAssemblers[i]->getQuadrature() == quad)
	return dynamic_cast<FirstOrderAssembler*>(subAssemblers[i]);
97 98 99 100
    }

    // check if all terms are pw_const
    bool pwConst = true;
101
    for (unsigned int i = 0; i < opTerms.size(); i++) {
102 103 104 105
      if (!(opTerms[i])->isPWConst()) {
	pwConst = false;
	break;
      }
106
    }
107 108

    // create new assembler
109
    if (!optimized) {
110
      newAssembler =
111
	(type == GRD_PSI) ?
Thomas Witkowski's avatar
Thomas Witkowski committed
112
	dynamic_cast<FirstOrderAssembler*>(new Stand10(op, assembler, quad)) :
113
	dynamic_cast<FirstOrderAssembler*>(new Stand01(op, assembler, quad));
114 115
    } else {
       if (pwConst) {
116
	 newAssembler =
117 118 119
	   (type == GRD_PSI) ?
	   dynamic_cast<FirstOrderAssembler*>(new Pre10(op, assembler, quad)) :
	   dynamic_cast<FirstOrderAssembler*>(new Pre01(op, assembler, quad));
120
       } else {
121
	 newAssembler =
122 123 124
	   (type == GRD_PSI) ?
	   dynamic_cast<FirstOrderAssembler*>(new Quad10(op, assembler, quad)) :
	   dynamic_cast<FirstOrderAssembler*>(new Quad01(op, assembler, quad));
125 126
       }
    }
127

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
128
    subAssemblers.push_back(newAssembler);
129
    return newAssembler;
130
  }
131

132

133
  Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad)
134
    : FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
Thomas Witkowski's avatar
Thomas Witkowski committed
135
  {
136 137
    name = "standard first order assembler";

138 139
    psi = rowFeSpace->getBasisFcts();
    phi = colFeSpace->getBasisFcts();
Thomas Witkowski's avatar
Thomas Witkowski committed
140
  }
141

142

143
  void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
144
  {
145
    mtl::dense_vector<double> grdPsi(dim + 1, 0.0);
146
    int nPoints = quadrature->getNumPoints();
147
    Lb.resize(nPoints);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
148
    vector<double> phival(nCol);
149

150 151 152 153 154 155
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].change_dim(dim + 1);
      Lb[iq] = 0.0;
    }
    for (unsigned int i = 0; i < terms.size(); i++)
      (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
156

157 158 159
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

160
      for (int i = 0; i < nCol; i++)
161
	phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
162 163 164

      for (int i = 0; i < nRow; i++) {
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
165
	for (int j = 0; j < nCol; j++)
166
	  mat[i][j] += quadrature->getWeight(iq) * phival[j] * dot(Lb[iq], grdPsi);
167 168 169 170
      }
    }
  }

171

172
  void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
173
  {
174
    mtl::dense_vector<double> grdPsi(dim + 1, 0.0);
175 176 177
    int nPoints = quadrature->getNumPoints();
    Lb.resize(nPoints);

178 179 180 181 182 183 184
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].change_dim(dim + 1);
      Lb[iq] = 0.0;
    }

    for (unsigned int i = 0; i < terms.size(); i++)
      (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
185

186 187 188 189 190
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

      for (int i = 0; i < nRow; i++) {
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
191
	vec[i] += quadrature->getWeight(iq) * dot(Lb[iq], grdPsi);
192 193 194 195
      }
    }
  }

196

197
  Quad10::Quad10(Operator *op, Assembler *assembler, Quadrature *quad)
198
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
199 200 201
  {
    name = "fast quadrature first order assembler";
  }
202 203


204
  void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
205 206
  {
    if (firstCall) {
Praetorius, Simon's avatar
Praetorius, Simon committed
207 208 209 210 211 212 213 214 215 216
      #ifdef _OPENMP
      #pragma omp critical
      #endif
      {
        const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
        psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
        basFcts = colFeSpace->getBasisFcts();
        phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
        firstCall = false;
      }
217 218 219
    }

    int nPoints = quadrature->getNumPoints();
220
    Lb.resize(nPoints);
221 222 223 224
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].change_dim(dim + 1);
      Lb[iq] = 0.0;
    }
225

226 227 228
    for (unsigned int i = 0; i < terms.size(); i++)
      (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);

229
    const mtl::dense2D<double>& phi = phiFast->getPhi();
230 231 232 233

    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

234
      const vector<mtl::dense_vector<double> >& grdPsi = psiFast->getGradient(iq);
235
      double factor = quadrature->getWeight(iq);
236

237
      for (int i = 0; i < nRow; i++) {
238
	for (int j = 0; j < nCol; j++)
239
	  mat[i][j] += factor * phi[iq][j] * dot(Lb[iq], grdPsi[i]);
240
      }
241 242 243
    }
  }

244

245
  void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
246 247
  {
    if (firstCall) {
Praetorius, Simon's avatar
Praetorius, Simon committed
248 249 250 251 252 253 254 255 256 257
      #ifdef _OPENMP
      #pragma omp critical
      #endif
      {
        const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
        psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
        basFcts = colFeSpace->getBasisFcts();
        phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
        firstCall = false;
      }
258 259 260 261
    }

    int nPoints = quadrature->getNumPoints();
    Lb.resize(nPoints);
262 263 264 265
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].change_dim(dim + 1);
      Lb[iq] = 0.0;
    }
266

267 268
    for (unsigned int i = 0; i < terms.size(); i++)
      (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
269

270 271
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();
272
      const vector<mtl::dense_vector<double> >& grdPsi = psiFast->getGradient(iq);
273

Thomas Witkowski's avatar
Thomas Witkowski committed
274
      for (int i = 0; i < nRow; i++)
275
	vec[i] += quadrature->getWeight(iq) * dot(Lb[iq], grdPsi[i]);
276 277
    }
  }
278

279

280
  Pre10::Pre10(Operator *op, Assembler *assembler, Quadrature *quad)
281 282
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
  {
283
    name = "precalculated first order assembler";
284 285 286
  }


287
  void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
288 289 290 291 292
  {
    const int *k;
    const double *values;

    if (firstCall) {
Praetorius, Simon's avatar
Praetorius, Simon committed
293 294 295 296 297 298 299 300 301 302
      #ifdef _OPENMP
      #pragma omp critical
      #endif
      {
        q10 = Q10PsiPhi::provideQ10PsiPhi(rowFeSpace->getBasisFcts(),
                                          colFeSpace->getBasisFcts(),
                                          quadrature);
        q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature);
        firstCall = false;
      }
303 304 305
    }

    const int **nEntries = q10->getNumberEntries();
306
    // Do not need do resize Lb, because it's size is always at least one.
307
    Lb[0] = 0.0;
308

309 310
    for (unsigned int i = 0; i < terms.size(); i++)
      (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
311 312 313 314 315 316 317 318

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

    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
	k = q10->getKVec(i, j);
	values = q10->getValVec(i, j);
	double val = 0.0;
319
	for (int m = 0; m < nEntries[i][j]; m++)
320
	  val += values[m] * Lb[0][k[m]];
321
	mat[i][j] += val;
322 323 324 325
      }
    }
  }

326

327
  Stand01::Stand01(Operator *op, Assembler *assembler, Quadrature *quad)
328
    : FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
329
  {
330
    FUNCNAME_DBG("Stand01::Stand01()");
331 332 333 334 335 336

    TEST_EXIT_DBG(dim > 0)("Should not happen!\n");

    grdPhi.resize(nCol);
    for (int i = 0; i < nCol; i++)
      grdPhi[i].change_dim(dim + 1);
Thomas Witkowski's avatar
Thomas Witkowski committed
337

338 339
    psi = rowFeSpace->getBasisFcts();
    phi = colFeSpace->getBasisFcts();
340
  }
341

342

343
  void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
344 345
  {
    int nPoints = quadrature->getNumPoints();
346
    Lb.resize(nPoints);
347 348 349 350
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].change_dim(dim + 1);
      Lb[iq] = 0.0;
    }
351

352 353
    for (int i = 0; i < static_cast<int>(terms.size()); i++)
      (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
354

355 356 357
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

358
      for (int i = 0; i < nCol; i++)
359 360 361 362 363
	(*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);

      for (int i = 0; i < nRow; i++) {
	double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
	for (int j = 0; j < nCol; j++)
364
	  mat[i][j] += quadrature->getWeight(iq) * psival * dot(Lb[iq], grdPhi[j]);
365
      }
366
    }
367 368
  }

369

370
  Quad01::Quad01(Operator *op, Assembler *assembler, Quadrature *quad)
371
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
372 373 374 375
  {
    name = "fast quadrature first order assembler";
  }

376

377
  void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
378 379
  {
    if (firstCall) {
Praetorius, Simon's avatar
Praetorius, Simon committed
380 381 382 383 384 385 386 387 388 389
      #ifdef _OPENMP
      #pragma omp critical
      #endif
      {
        const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
        psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
        basFcts = colFeSpace->getBasisFcts();
        phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
        firstCall = false;
      }
390 391 392
    }

    int nPoints = quadrature->getNumPoints();
393
    Lb.resize(nPoints);
394 395 396 397
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].change_dim(dim + 1);
      Lb[iq] = 0.0;
    }
398

399 400
    for (int i = 0; i < static_cast<int>(terms.size()); i++)
      (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
401

402
    const mtl::dense2D<double>& psi = psiFast->getPhi();
403

404
    mtl::dense2D<double> facMat(nPoints, nCol);
405
    for (int iq = 0; iq < nPoints; iq++) {
406 407 408
      double weight = quadrature->getWeight(iq);
      mtl::dense_vector<double>& Lb_iq = Lb[iq];
      const vector<mtl::dense_vector<double> >& grdPsi = phiFast->getGradient(iq);
409

410 411 412 413
      Lb_iq *= elInfo->getDet();

      for (int i = 0; i < nCol; i++)
	facMat[iq][i] = weight * dot(Lb_iq, grdPsi[i]);
414
    }
415 416

    mat += trans(psi) * facMat;
417 418
  }

419

420
  Pre01::Pre01(Operator *op, Assembler *assembler, Quadrature *quad)
421
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
422 423 424 425
  {
    name = "precalculated first order assembler";
  }

426

427
  void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
428 429 430 431 432
  {
    const int *l;
    const double *values;

    if (firstCall) {
Praetorius, Simon's avatar
Praetorius, Simon committed
433 434 435 436 437 438 439 440 441 442
      #ifdef _OPENMP
      #pragma omp critical
      #endif
      {
        q01 = Q01PsiPhi::provideQ01PsiPhi(rowFeSpace->getBasisFcts(),
                                          colFeSpace->getBasisFcts(),
                                          quadrature);
        q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature);
        firstCall = false;
      }
443 444 445
    }

    const int **nEntries = q01->getNumberEntries();
446
    // Do not need to resize Lb, because it's size is always at least one!
447
    Lb[0] = 0.0;
448

449 450
    for (int i = 0; i < static_cast<int>( terms.size()); i++)
      (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
451 452 453 454 455 456 457 458 459 460 461

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

    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
	l = q01->getLVec(i, j);
	values = q01->getValVec(i, j);
	double val = 0.0;
	for (int m = 0; m < nEntries[i][j]; m++) {
	  val += values[m] * Lb[0][l[m]];
	}
462
	mat[i][j] += val;
463 464 465 466
      }
    }
  }

467

468
  void Pre10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
469 470 471 472 473
  {
    const int *k;
    const double *values;

    if (firstCall) {
Praetorius, Simon's avatar
Praetorius, Simon committed
474 475 476 477 478 479 480 481 482 483
      #ifdef _OPENMP
      #pragma omp critical
      #endif
      {
        q10 = Q10PsiPhi::provideQ10PsiPhi(rowFeSpace->getBasisFcts(),
                                          colFeSpace->getBasisFcts(),
                                          quadrature);
        q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature);
        firstCall = false;
      }
484 485 486
    }

    const int *nEntries = q1->getNumberEntries();
487
    // Do not need to resize Lb, because it's size is always at least one!
488
    Lb[0] = 0.0;
489

490 491
    for (int i = 0; i < static_cast<int>(terms.size()); i++)
      (static_cast<FirstOrderTerm*>(terms[i]))->getLb(elInfo, Lb);
492 493 494 495 496 497 498 499 500 501

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

    for (int i = 0; i < nRow; i++) {
      k = q1->getKVec(i);
      values = q1->getValVec(i);
      double val = 0.0;
      for (int m = 0; m < nEntries[i]; m++) {
	val += values[m] * Lb[0][k[m]];
      }
502
      vec[i] += val;
503 504 505 506
    }
  }

}