Operator.inc.hpp 7.31 KB
Newer Older
1
2
3
4
5
6
7
#pragma once

#include <vector>

namespace AMDiS
{
  template <class MeshView>
8
9
10
    template <class RowFeSpace, class ColFeSpace>
  void Operator<MeshView>::init(RowFeSpace const& rowFeSpace,
				ColFeSpace const& colFeSpace)
11
  {
12
13
//     const auto& rowFE = rowView.tree().finiteElement();
//     const auto& colFE = colView.tree().finiteElement();
14
    
15
16
17
18
19
//     psiDegree = rowFE.localBasis().order();
//     phiDegree = colFE.localBasis().order();
      
      psiDegree = GetDegree<RowFeSpace>::value;
      phiDegree = GetDegree<ColFeSpace>::value;
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
    
    // TODO: calc quadrature degree here.
  }
  

  template <class MeshView>
    template <class RowView, class ColView, class ElementMatrix>
  void Operator<MeshView>::getElementMatrix(RowView const& rowView,
					    ColView const& colView,
					    ElementMatrix& elementMatrix)
  {    
    assembleZeroOrderTerms(rowView, colView, elementMatrix);
    assembleSecondOrderTerms(rowView, colView, elementMatrix);
  }
  
35
36
37
38
39
40
41
42
43
44
  
  template <class MeshView>
    template <class RowView, class ElementVector>
  void Operator<MeshView>::getElementVector(RowView const& rowView,
					    ElementVector& elementVector)
  {    
    assembleZeroOrderTerms(rowView, elementVector);
  }
  
      
45
46
47
48
49
50
51
52
53
54
55
  template <class MeshView>
    template <class RowView, class ColView, class ElementMatrix>
  void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
						  ColView const& colView,
						  ElementMatrix& elementMatrix)
  {
    using Element = typename RowView::Element;
    auto element = rowView.element();
    const int dim = Element::dimension;
    auto geometry = element.geometry();
    
56
57
    const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
    const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
58
59
    
    int order = getQuadratureDegree(0);
60
    const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
61
62
63
64
65
66
67
    
    std::vector<double> functionValues(quad.size(), 0.0);
    for (auto* operatorTerm : zeroOrder)
	operatorTerm->evalAtQPs(element, quad, functionValues);
    
    for (size_t iq = 0; iq < quad.size(); ++iq) {
	// Position of the current quadrature point in the reference element
68
	const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
69
70
71
72

	// The multiplicative factor in the integral transformation formula
	const double integrationElement = geometry.integrationElement(quadPos);
	
73
	std::vector<Dune::FieldVector<double,1> > rowShapeValues, colShapeValues;
74
	rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
75
76
77
78
	
	if (psiDegree == phiDegree)
	  colShapeValues = rowShapeValues;
	else
79
	  colLocalBasis.evaluateFunction(quadPos, colShapeValues);
80
81
82
83
84
85
86
87
88
89
90
91

	for (size_t i = 0; i < elementMatrix.N(); ++i) {
	    for (size_t j = 0; j < elementMatrix.M(); ++j) {
		const int local_i = rowView.tree().localIndex(i);
		const int local_j = colView.tree().localIndex(j);
		elementMatrix[local_i][local_j] += functionValues[iq] * (rowShapeValues[i] * colShapeValues[j]) * quad[iq].weight() * integrationElement;
	    }
	}
    }
  }
  

92
93
94
95
96
97
98
99
100
101
102
      
  template <class MeshView>
    template <class RowView, class ElementVector>
  void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
						  ElementVector& elementvector)
  {
    using Element = typename RowView::Element;
    auto element = rowView.element();
    const int dim = Element::dimension;
    auto geometry = element.geometry();
    
103
    const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
    
    int order = getQuadratureDegree(0);
    const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
    
    std::vector<double> functionValues(quad.size(), 0.0);
    for (auto* operatorTerm : zeroOrder)
	operatorTerm->evalAtQPs(element, quad, functionValues);
    
    for (size_t iq = 0; iq < quad.size(); ++iq) {
	// 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
	const double integrationElement = geometry.integrationElement(quadPos);
	
	std::vector<Dune::FieldVector<double,1> > rowShapeValues;
120
	rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
121
122
123
124
125
126
127
128
129

	for (size_t i = 0; i < elementvector.size(); ++i) {
	    const int local_i = rowView.tree().localIndex(i);
	    elementvector[local_i] += functionValues[iq] * rowShapeValues[i] * quad[iq].weight() * integrationElement;
	}
    }
  }
  

130
131
132
133
134
135
136
137
138
139
140
  template <class MeshView>
    template <class RowView, class ColView, class ElementMatrix>
  void Operator<MeshView>::assembleSecondOrderTerms(RowView const& rowView,
						    ColView const& colView,
						    ElementMatrix& elementMatrix)
  {
    using Element = typename RowView::Element;
    auto element = rowView.element();
    const int dim = Element::dimension;
    auto geometry = element.geometry();
    
141
142
    const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
    const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
143
    
144
    int order = getQuadratureDegree(2);
145
    const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
146
147
148
149
150
151
152
153
154
155
    
    std::vector<double> functionValues(quad.size(), 0.0);
    for (auto* operatorTerm : secondOrder)
	operatorTerm->evalAtQPs(element, quad, functionValues);
    
    // TODO: currently only the implementation for equal fespaces
    assert( psiDegree == phiDegree );
    
    for (size_t iq = 0; iq < quad.size(); ++iq) {
	// Position of the current quadrature point in the reference element
156
	const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
157
158
159
160
161
162
163
164

	// 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
	const double integrationElement = geometry.integrationElement(quadPos);

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

	// Compute the shape function gradients on the real element
169
	std::vector<Dune::FieldVector<double,dim> > gradients(referenceGradients.size());
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187

	for (size_t i = 0; i < gradients.size(); ++i)
	    jacobian.mv(referenceGradients[i][0], gradients[i]);

	for (size_t i = 0; i < elementMatrix.N(); ++i) {
	    for (size_t j = 0; j < elementMatrix.M(); ++j) {
		const int local_i = rowView.tree().localIndex(i);
		const int local_j = colView.tree().localIndex(j);
		elementMatrix[local_i][local_j] += functionValues[iq] * (gradients[i] * gradients[j]) * quad[iq].weight() * integrationElement;
	    }
	}
    }
  }
  

  template <class MeshView>
  int Operator<MeshView>::getQuadratureDegree(int order/*, FirstOrderType firstOrderType*/)
  {
188
    std::list<OperatorTermType*>* terms = NULL;
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207

    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)
208
	maxTermDegree = std::max(maxTermDegree, term->getDegree());
209
210
211
212
213

    return psiDegree + phiDegree - order + maxTermDegree;
  }
  
} // end namespace AMDiS