DUNE-FEM (unstable)

raviartthomas4cube2dlocalinterpolation.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS4_CUBE2D_LOCALINTERPOLATION_HH
6#define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS4_CUBE2D_LOCALINTERPOLATION_HH
7
8#include <vector>
9
11
12namespace Dune
13{
14
23 template<class LB>
25 {
26
27 public:
28
34 RT4Cube2DLocalInterpolation (unsigned int s = 0)
35 {
36 sign0 = sign1 = sign2 = sign3 = 1.0;
37 if (s & 1)
38 {
39 sign0 *= -1.0;
40 }
41 if (s & 2)
42 {
43 sign1 *= -1.0;
44 }
45 if (s & 4)
46 {
47 sign2 *= -1.0;
48 }
49 if (s & 8)
50 {
51 sign3 *= -1.0;
52 }
53
54 n0[0] = -1.0;
55 n0[1] = 0.0;
56 n1[0] = 1.0;
57 n1[1] = 0.0;
58 n2[0] = 0.0;
59 n2[1] = -1.0;
60 n3[0] = 0.0;
61 n3[1] = 1.0;
62 }
63
72 template<typename F, typename C>
73 void interpolate (const F& f, std::vector<C>& out) const
74 {
75 // f gives v*outer normal at a point on the edge!
76 typedef typename LB::Traits::RangeFieldType Scalar;
77 typedef typename LB::Traits::DomainFieldType Vector;
78
79 out.resize(60);
80 fill(out.begin(), out.end(), 0.0);
81
82 const int qOrder = 12;
84
85 for (typename QuadratureRule<Scalar,1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
86 {
87 Scalar qPos = it->position();
88 typename LB::Traits::DomainType localPos;
89
90 localPos[0] = 0.0;
91 localPos[1] = qPos;
92 auto y = f(localPos);
93 out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0;
94 out[1] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight();
95 out[2] += (y[0]*n0[0] + y[1]*n0[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign0;
96 out[3] += (y[0]*n0[0] + y[1]*n0[1])*(20.0*qPos*qPos*qPos - 30.0*qPos*qPos + 12.0*qPos - 1.0)*it->weight();
97 out[4] += (y[0]*n0[0] + y[1]*n0[1])*(1.0-20.0*qPos+90.0*pow(qPos,2)-140.0*pow(qPos,3)+70.0*pow(qPos,4))*it->weight()*sign0;
98
99 localPos[0] = 1.0;
100 localPos[1] = qPos;
101 y = f(localPos);
102 out[5] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1;
103 out[6] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight();
104 out[7] += (y[0]*n1[0] + y[1]*n1[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign1;
105 out[8] += (y[0]*n1[0] + y[1]*n1[1])*(-20.0*qPos*qPos*qPos + 30.0*qPos*qPos - 12.0*qPos + 1.0)*it->weight();
106 out[9] += (y[0]*n1[0] + y[1]*n1[1])*(1.0-20.0*qPos+90.0*pow(qPos,2)-140.0*pow(qPos,3)+70.0*pow(qPos,4))*it->weight()*sign1;
107
108 localPos[0] = qPos;
109 localPos[1] = 0.0;
110 y = f(localPos);
111 out[10] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2;
112 out[11] += (y[0]*n2[0] + y[1]*n2[1])*(1.0 - 2.0*qPos)*it->weight();
113 out[12] += (y[0]*n2[0] + y[1]*n2[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign2;
114 out[13] += (y[0]*n2[0] + y[1]*n2[1])*(-20.0*qPos*qPos*qPos + 30.0*qPos*qPos - 12.0*qPos + 1.0)*it->weight();
115 out[14] += (y[0]*n2[0] + y[1]*n2[1])*(1.0-20.0*qPos+90.0*pow(qPos,2)-140.0*pow(qPos,3)+70.0*pow(qPos,4))*it->weight()*sign2;
116
117 localPos[0] = qPos;
118 localPos[1] = 1.0;
119 y = f(localPos);
120 out[15] += (y[0]*n3[0] + y[1]*n3[1])*it->weight()*sign3;
121 out[16] += (y[0]*n3[0] + y[1]*n3[1])*(2.0*qPos - 1.0)*it->weight();
122 out[17] += (y[0]*n3[0] + y[1]*n3[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign3;
123 out[18] += (y[0]*n3[0] + y[1]*n3[1])*(20.0*qPos*qPos*qPos - 30.0*qPos*qPos + 12.0*qPos - 1.0)*it->weight();
124 out[19] += (y[0]*n3[0] + y[1]*n3[1])*(1.0-20.0*qPos+90.0*pow(qPos,2)-140.0*pow(qPos,3)+70.0*pow(qPos,4))*it->weight()*sign3;
125 }
126
128
129 for (typename QuadratureRule<Vector,2>::const_iterator it = rule2.begin();
130 it != rule2.end(); ++it)
131 {
132 FieldVector<double,2> qPos = it->position();
133
134 auto y = f(qPos);
135 std::vector<std::vector<double> > l(2,std::vector<double> (5));
136 l[0][0]=1.0;
137 l[1][0]=1.0;
138 l[0][1]=2.0*qPos[0]-1.0;
139 l[1][1]=2.0*qPos[1]-1.0;
140 l[0][2]=6.0*qPos[0]*qPos[0]-6.0*qPos[0]+1.0;
141 l[1][2]=6.0*qPos[1]*qPos[1]-6.0*qPos[1]+1.0;
142 l[0][3]=20.0*qPos[0]*qPos[0]*qPos[0] - 30.0*qPos[0]*qPos[0] + 12.0*qPos[0] - 1.0;
143 l[1][3]=20.0*qPos[1]*qPos[1]*qPos[1] - 30.0*qPos[1]*qPos[1] + 12.0*qPos[1] - 1.0;
144 l[0][4]=1.0-20.0*qPos[0]+90.0*pow(qPos[0],2)-140.0*pow(qPos[0],3)+70.0*pow(qPos[0],4);
145 l[1][4]=1.0-20.0*qPos[1]+90.0*pow(qPos[1],2)-140.0*pow(qPos[1],3)+70.0*pow(qPos[1],4);
146
147 for (int i=0;i<4;i++)
148 for (int j=0;j<5;j++)
149 out[20+i*5+j] +=y[0]*l[0][i]*l[1][j]*it->weight();
150
151 for (int i=0;i<5;i++)
152 for (int j=0;j<4;j++)
153 out[40+i*4+j] +=y[1]*l[0][i]*l[1][j]*it->weight();
154 }
155 }
156
157 private:
158 typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3;
159 typename LB::Traits::DomainType n0, n1, n2, n3;
160 };
161}
162
163#endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS3_CUBE2D_LOCALINTERPOLATION_HH
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:214
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:326
Second order Raviart-Thomas shape functions on the reference triangle.
Definition: raviartthomas4cube2dlocalinterpolation.hh:25
RT4Cube2DLocalInterpolation(unsigned int s=0)
Make set number s, where 0 <= s < 16.
Definition: raviartthomas4cube2dlocalinterpolation.hh:34
void interpolate(const F &f, std::vector< C > &out) const
Interpolate a given function with shape functions.
Definition: raviartthomas4cube2dlocalinterpolation.hh:73
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: interface.hh:235
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)