Operator.inc.hpp 19.8 KB
Newer Older
1
2
3
4
5
#pragma once

namespace AMDiS
{
  template <class MeshView>
6
7
8
    template <class RowFeSpace, class ColFeSpace>
  void Operator<MeshView>::init(RowFeSpace const& rowFeSpace,
				ColFeSpace const& colFeSpace)
9
  {
10
    AMDIS_FUNCNAME("Operator::init()");
11

12
13
14
    using IdType = typename Operator<MeshView>::IdType;
    lastMatrixId = std::numeric_limits<IdType>::max();
    lastVectorId = std::numeric_limits<IdType>::max();
15

16
17
//     auto const& rowFE = rowView.tree().finiteElement();
//     auto const& colFE = colView.tree().finiteElement();
18

19
20
//     psiDegree = rowFE.localBasis().order();
//     phiDegree = colFE.localBasis().order();
21

22
23
    psiDegree = getPolynomialDegree<RowFeSpace>;
    phiDegree = getPolynomialDegree<ColFeSpace>;
24

25
26
    // TODO: calc quadrature degree here.
  }
27

28

29
  template <class MeshView>
30
31
    template <class Element>
  int Operator<MeshView>::getQuadratureDegree(Element const& element, int order, FirstOrderType firstOrderType)
32
33
34
  {
    std::list<OperatorTermType*>* terms = NULL;

35
36
37
    auto const& t = element.type();
    auto geometry = element.geometry();

38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
    switch(order)
    {
    case 0:
      terms = &zeroOrder;
      break;
    case 1:
      if (firstOrderType == GRD_PHI)
        terms = &firstOrderGrdPhi;
      else
        terms = &firstOrderGrdPsi;
      break;
    case 2:
      terms = &secondOrder;
      break;
    }
    int maxTermDegree = 0;

    for (OperatorTermType* term : *terms)
56
57
58
59
60
61
62
63
      maxTermDegree = std::max(maxTermDegree, term->getDegree(t));

    int degree = psiDegree + phiDegree + maxTermDegree;
    if (t.isSimplex())
      degree -= order;

    if (!geometry.affine())
      degree += 1; // oder += (order+1)
64

65
    return degree;
66
  }
67
68


69
70
  template <class MeshView>
    template <class LocalView>
71
  typename Operator<MeshView>::IdType
72
73
74
75
  Operator<MeshView>::getElementId(LocalView const& localView)
  {
    auto const& element = localView.element();
    auto const& grid = localView.globalBasis().gridView().grid();
76

77
78
79
    auto const& idSet = grid.localIdSet();
    return idSet.id(element);
  }
80

81

82
  template <class MeshView>
83
    template <class RowView, class ColView>
84
  bool Operator<MeshView>::getElementMatrix(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
85
86
                                            ColView const& colView,
                                            ElementMatrix& elementMatrix,
87
                                            double* factor)
88
  {
89
    AMDIS_FUNCNAME("Operator::getElementMatrix()");
90

91
    double fac = factor ? *factor : 1.0;
92

93
94
    if (fac == 0.0)
      return false;
95

96
97
    const auto nRows = rowView.tree().finiteElement().size();
    const auto nCols = colView.tree().finiteElement().size();
98

99
100
    auto currentId = getElementId(rowView);
    bool useCachedMatrix = (lastMatrixId == currentId
101
                            && num_rows(cachedElementMatrix) == nRows
102
                            && num_cols(cachedElementMatrix) == nCols);
103

104
105
    bool assign = true;
    if (!useCachedMatrix) {
106

107
108
      cachedElementMatrix.change_dim(nRows, nCols);
      set_to_zero(cachedElementMatrix);
109

110
111
112
113
114
115
116
117
      if (!zeroOrder.empty())
        assembleZeroOrderTerms(rowView, colView, cachedElementMatrix);
      if (!firstOrderGrdPhi.empty())
        assembleFirstOrderTermsGrdPhi(rowView, colView, cachedElementMatrix);
      if (!firstOrderGrdPsi.empty())
        assembleFirstOrderTermsGrdPsi(rowView, colView, cachedElementMatrix);
      if (!secondOrder.empty())
        assembleSecondOrderTerms(rowView, colView, cachedElementMatrix);
118

119
120
121
      assign = !zeroOrder.empty() || !firstOrderGrdPhi.empty() ||
               !firstOrderGrdPsi.empty() || !secondOrder.empty();
    }
122

123
124
    if (assign)
      elementMatrix += fac * cachedElementMatrix;
125

126
    lastMatrixId = currentId;
127

128
    return true;
129
  }
130
131


132
  template <class MeshView>
133
    template <class RowView>
134
  bool Operator<MeshView>::getElementVector(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
135
					                                  ElementVector& elementVector,
136
                                            double* factor)
137
  {
138
    AMDIS_FUNCNAME("Operator::getElementVector()");
139

140
    double fac = factor ? *factor : 1.0;
141

142
143
    if (fac == 0.0)
      return false;
144
145

    test_exit( firstOrderGrdPhi.empty() && secondOrder.empty(),
146
               "Derivatives on ansatz-functions not allowed on the vector-side!");
147
148
149

    const auto nRows = rowView.tree().finiteElement().size();

150
151
    auto currentId = getElementId(rowView);
    bool useCachedVector = (lastVectorId == currentId && size(cachedElementVector) == nRows);
152

153
    bool assign = true;
154
    if (!useCachedVector) {
155
156
      cachedElementVector.change_dim(nRows);
      set_to_zero(cachedElementVector);
157

158
159
160
161
      if (!zeroOrder.empty())
        assembleZeroOrderTerms(rowView, cachedElementVector);
      if (!firstOrderGrdPsi.empty())
        assembleFirstOrderTermsGrdPsi(rowView, cachedElementVector);
162

163
164
      assign = !zeroOrder.empty() || !firstOrderGrdPsi.empty();
    }
165

166
167
    if (assign)
      elementVector += fac * cachedElementVector;
168

169
    lastVectorId = currentId;
170

171
    return true;
172
  }
173
174


175
  template <class MeshView>
176
    template <class RowView, class ColView>
177
  void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
178
179
                                                  ColView const& colView,
                                                  ElementMatrix& elementMatrix)
180
  {
181
    AMDIS_FUNCNAME("Operator::assembleZeroOrderTerms(elementMatrix)");
182

183
184
185
186
    using Element = typename RowView::Element;
    auto element = rowView.element();
    const int dim = Element::dimension;
    auto geometry = element.geometry();
187

188
189
    auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
    auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
190

191
    int order = getQuadratureDegree(element, 0);
192
    auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
193

194
    for (auto* operatorTerm : zeroOrder)
195
      operatorTerm->init(element, quad);
196

197
    for (std::size_t iq = 0; iq < quad.size(); ++iq) {
198
      // Position of the current quadrature point in the reference element
Praetorius, Simon's avatar
Praetorius, Simon committed
199
      Dune::FieldVector<double,dim> const& quadPos = quad[iq].position();
200
201

      // The multiplicative factor in the integral transformation formula
202
      const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
203

204
205
      std::vector<Dune::FieldVector<double,1> > rowShapeValues, colShapeValues;
      rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
206

207
208
209
210
211
      if (psiDegree == phiDegree)
        colShapeValues = rowShapeValues;
      else
        colLocalBasis.evaluateFunction(quadPos, colShapeValues);

212
213
      for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
        for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
214
215
216
          const int local_i = rowView.tree().localIndex(i);
          const int local_j = colView.tree().localIndex(j);
          for (auto* operatorTerm : zeroOrder)
217
            elementMatrix[local_i][local_j]
218
              += operatorTerm->evalZot(iq, rowShapeValues[i], colShapeValues[j]) * factor;
219
220
        }
      }
221
222
223
    }
  }

224
225


226
  template <class MeshView>
227
    template <class RowView>
228
  void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
229
						                                      ElementVector& elementvector)
230
  {
231
    AMDIS_FUNCNAME("Operator::assembleZeroOrderTerms(elementvector)");
232

233
234
235
236
    using Element = typename RowView::Element;
    auto element = rowView.element();
    const int dim = Element::dimension;
    auto geometry = element.geometry();
237

238
    auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
239

240
    int order = getQuadratureDegree(element, 0);
241
    auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
242

243
    for (auto* operatorTerm : zeroOrder)
Praetorius, Simon's avatar
Praetorius, Simon committed
244
	    operatorTerm->init(element, quad);
245

246
    for (std::size_t iq = 0; iq < quad.size(); ++iq) {
247
248
249
250
      // Position of the current quadrature point in the reference element
      const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();

      // The multiplicative factor in the integral transformation formula
251
      const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
252

253
254
255
      std::vector<Dune::FieldVector<double,1> > rowShapeValues;
      rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);

256
      for (std::size_t i = 0; i < size(elementvector); ++i) {
257
258
        const int local_i = rowView.tree().localIndex(i);
        for (auto* operatorTerm : zeroOrder)
259
          elementvector[local_i] += operatorTerm->evalZot(iq, rowShapeValues[i]) * factor;
260
      }
261
262
    }
  }
263

264
265

  template <class MeshView>
266
    template <class RowView, class ColView>
267
  void Operator<MeshView>::assembleFirstOrderTermsGrdPhi(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
268
269
                                                         ColView const& colView,
                                                         ElementMatrix& elementMatrix)
270
  {
271
    AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPhi(elementMatrix)");
272

273
274
    auto element = rowView.element();
    auto geometry = element.geometry();
275
276
    const int dim = RowView::GridView::dimension;
    const int dow = RowView::GridView::dimensionworld;
277

278
279
    auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
    auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
280

281
    int order = getQuadratureDegree(element, 1, GRD_PHI);
282
    auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
283

284
    for (auto* operatorTerm : firstOrderGrdPhi)
285
      operatorTerm->init(element, quad);
286

287
    for (std::size_t iq = 0; iq < quad.size(); ++iq) {
288
289
290
291
292
293
294
      // Position of the current quadrature point in the reference element
      const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();

      // The transposed inverse Jacobian of the map from the reference element to the element
      const auto jacobian = geometry.jacobianInverseTransposed(quadPos);

      // The multiplicative factor in the integral transformation formula
295
      const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
296
297
298

      std::vector<Dune::FieldVector<double,1> > rowShapeValues;
      rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
299

300
301
302
303
304
      // The gradients of the shape functions on the reference element
      std::vector<Dune::FieldMatrix<double,1,dim> > colReferenceGradients;
      colLocalBasis.evaluateJacobian(quadPos, colReferenceGradients);

      // Compute the shape function gradients on the real element
305
      std::vector<Dune::FieldVector<double,dow> > colGradients(colReferenceGradients.size());
306

307
      for (std::size_t i = 0; i < colGradients.size(); ++i)
308
309
        jacobian.mv(colReferenceGradients[i][0], colGradients[i]);

310
311
      for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
        for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
312
313
314
          const int local_i = rowView.tree().localIndex(i);
          const int local_j = colView.tree().localIndex(j);
          for (auto* operatorTerm : firstOrderGrdPhi)
315
            elementMatrix[local_i][local_j]
316
              += operatorTerm->evalFot1(iq, rowShapeValues[i], colGradients[j]) * factor;
317
318
        }
      }
319
320
    }
  }
321

322
323

  template <class MeshView>
324
    template <class RowView, class ColView>
325
  void Operator<MeshView>::assembleFirstOrderTermsGrdPsi(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
326
327
                                                         ColView const& colView,
                                                         ElementMatrix& elementMatrix)
328
  {
329
    AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPsi(elementMatrix)");
330

331
332
    auto element = rowView.element();
    auto geometry = element.geometry();
333
334
    const int dim = RowView::GridView::dimension;
    const int dow = RowView::GridView::dimensionworld;
335

336
337
    auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
    auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
338

339
    int order = getQuadratureDegree(element, 1, GRD_PSI);
340
    auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
341

342
    for (auto* operatorTerm : firstOrderGrdPsi)
343
      operatorTerm->init(element, quad);
344

345
    for (std::size_t iq = 0; iq < quad.size(); ++iq) {
346
347
348
349
350
351
352
      // Position of the current quadrature point in the reference element
      const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();

      // The transposed inverse Jacobian of the map from the reference element to the element
      const auto jacobian = geometry.jacobianInverseTransposed(quadPos);

      // The multiplicative factor in the integral transformation formula
353
      const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
354

355
356
357
358
359
      // The gradients of the shape functions on the reference element
      std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
      rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);

      // Compute the shape function gradients on the real element
360
      std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());
361

362
      for (std::size_t i = 0; i < rowGradients.size(); ++i)
363
364
365
366
367
        jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);

      std::vector<Dune::FieldVector<double,1> > colShapeValues;
      colLocalBasis.evaluateFunction(quadPos, colShapeValues);

368
369
      for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
        for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
370
371
372
          const int local_i = rowView.tree().localIndex(i);
          const int local_j = colView.tree().localIndex(j);
          for (auto* operatorTerm : firstOrderGrdPsi)
373
            elementMatrix[local_i][local_j]
374
              += operatorTerm->evalFot2(iq, rowGradients[i], colShapeValues[j]) * factor;
375
376
        }
      }
377
378
379
    }
  }

380
381


382
  template <class MeshView>
383
    template <class RowView>
384
  void Operator<MeshView>::assembleFirstOrderTermsGrdPsi(RowView const& rowView,
385
                                                         ElementVector& elementvector)
386
387
  {
    AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPsi(elementvector)");
388

389
    auto element = rowView.element();
390
391
    const int dim = RowView::GridView::dimension;
    const int dow = RowView::GridView::dimensionworld;
392
    auto geometry = element.geometry();
393

394
    auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
395

396
    int order = getQuadratureDegree(element, 1, GRD_PSI);
397
    auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
398

399
400
    for (auto* operatorTerm : firstOrderGrdPsi)
        operatorTerm->init(element, quad);
401

402
    for (std::size_t iq = 0; iq < quad.size(); ++iq) {
403
404
405
406
407
408
409
      // Position of the current quadrature point in the reference element
      const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();

      // The transposed inverse Jacobian of the map from the reference element to the element
      const auto jacobian = geometry.jacobianInverseTransposed(quadPos);

      // The multiplicative factor in the integral transformation formula
410
      const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
411

412
413
414
415
416
      // The gradients of the shape functions on the reference element
      std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
      rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);

      // Compute the shape function gradients on the real element
417
      std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());
418

419
      for (std::size_t i = 0; i < rowGradients.size(); ++i)
420
421
        jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);

422
      for (std::size_t i = 0; i < size(elementvector); ++i) {
423
424
        const int local_i = rowView.tree().localIndex(i);
        for (auto* operatorTerm : firstOrderGrdPsi)
425
          elementvector[local_i] += operatorTerm->evalFot2(iq, rowGradients[i]) * factor;
426
427
428
      }
    }
  }
429

430

431
  template <class MeshView>
432
    template <class RowView, class ColView>
433
  void Operator<MeshView>::assembleSecondOrderTerms(RowView const& rowView,
Praetorius, Simon's avatar
Praetorius, Simon committed
434
435
                                                    ColView const& colView,
                                                    ElementMatrix& elementMatrix)
436
  {
437
    AMDIS_FUNCNAME("Operator::assembleSecondOrderTerms(elementMatrix)");
438

439
    auto element = rowView.element();
440
441
    const int dim = RowView::GridView::dimension;
    const int dow = RowView::GridView::dimensionworld;
442
    auto geometry = element.geometry();
443

444
445
    auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
//     auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
446

447
    int order = getQuadratureDegree(element, 2);
448
    auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
449

450
    for (auto* operatorTerm : secondOrder)
451
      operatorTerm->init(element, quad);
452

453
454
    // TODO: currently only the implementation for equal fespaces
    assert( psiDegree == phiDegree );
455

456
    for (std::size_t iq = 0; iq < quad.size(); ++iq) {
457
458
459
460
461
462
463
      // Position of the current quadrature point in the reference element
      const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();

      // The transposed inverse Jacobian of the map from the reference element to the element
      const auto jacobian = geometry.jacobianInverseTransposed(quadPos);

      // The multiplicative factor in the integral transformation formula
464
      const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
465
466
467
468
469
470

      // The gradients of the shape functions on the reference element
      std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients;
      rowLocalBasis.evaluateJacobian(quadPos, referenceGradients);

      // Compute the shape function gradients on the real element
471
      std::vector<Dune::FieldVector<double,dow> > gradients(referenceGradients.size());
472

473
      for (std::size_t i = 0; i < gradients.size(); ++i)
474
475
        jacobian.mv(referenceGradients[i][0], gradients[i]);

476
477
      for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
        for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
478
479
480
          const int local_i = rowView.tree().localIndex(i);
          const int local_j = colView.tree().localIndex(j);
          for (auto* operatorTerm : secondOrder)
481
            elementMatrix[local_i][local_j]
482
              += operatorTerm->evalSot(iq, gradients[i], gradients[j]) * factor;
483
484
        }
      }
485
486
    }
  }
487
488
489



490
  template <class MeshView>
491
492
    template <class Term>
  Operator<MeshView>& Operator<MeshView>::addZOTImpl(Term const& term)
493
  {
494
495
496
    zeroOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
    return *this;
  }
497
498


499
500
  template <class MeshView>
    template <class Term>
501
  Operator<MeshView>& Operator<MeshView>::addFOTImpl(Term const& term,
502
503
504
505
506
507
508
509
510
                                                     FirstOrderType firstOrderType)
  {
    using OpTerm = GenericOperatorTerm<MeshView, Term>;
    if (firstOrderType == GRD_PHI)
      firstOrderGrdPhi.push_back(new OpTerm(term));
    else
      firstOrderGrdPsi.push_back(new OpTerm(term));
    return *this;
  }
511
512


513
  template <class MeshView>
514
    template <class Term>
515
  Operator<MeshView>& Operator<MeshView>::addFOTImpl(Term const& term,
516
                                                     std::size_t i,
517
518
                                                     FirstOrderType firstOrderType)
  {
519
    using OpTerm = GenericOperatorTerm<MeshView, Term, VectorComponent>;
520

521
    if (firstOrderType == GRD_PHI)
522
      firstOrderGrdPhi.push_back(new OpTerm(term, {i}));
523
    else
524
      firstOrderGrdPsi.push_back(new OpTerm(term, {i}));
525
526
    return *this;
  }
527
528


529
530
531
532
533
534
535
  template <class MeshView>
    template <class Term>
  Operator<MeshView>& Operator<MeshView>::addSOTImpl(Term const& term)
  {
    secondOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
    return *this;
  }
536
537


538
  template <class MeshView>
539
    template <class Term>
540
  Operator<MeshView>& Operator<MeshView>::addSOTImpl(Term const& term,
541
                                                     std::size_t i, std::size_t j)
542
  {
543
544
    using OpTerm = GenericOperatorTerm<MeshView, Term, MatrixComponent>;
    secondOrder.push_back(new OpTerm(term, {i,j}));
545
    return *this;
546
  }
547

548
} // end namespace AMDiS