/[escript]/trunk/speckley/src/RectangleIntegrals.cpp
ViewVC logotype

Contents of /trunk/speckley/src/RectangleIntegrals.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5478 - (show annotations)
Wed Feb 18 01:57:49 2015 UTC (4 years, 7 months ago) by sshaw
File size: 9692 byte(s)
introducing shiny new reduced function spaces to speckley
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2015 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16 #define ESNEEDPYTHON
17 #include "esysUtils/first.h"
18
19 #include <esysUtils/index.h>
20 #include <speckley/Rectangle.h>
21
22 namespace speckley {
23 void Rectangle::integral_order2(std::vector<double>& integrals, const escript::Data& arg) const {
24 const double weights[] = {0.333333333333, 1.33333333333, 0.333333333333};
25 const int numComp = arg.getDataPointSize();
26 const double volume_product = 0.25*m_dx[0]*m_dx[1];
27 for (int ei = 0; ei < m_NE[1]; ++ei) {
28 for (int ej = 0; ej < m_NE[0]; ++ej) {
29 const double *e = arg.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
30 double result = 0;
31 for (int comp = 0; comp < numComp; ++comp) {
32 for (int i = 0; i < 3; ++i) {
33 for (int j = 0; j < 3; ++j) {
34 result += weights[i] * weights[j] * e[INDEX3(comp,i,j,numComp,3)];
35 }
36 }
37 integrals[comp] += result;
38 }
39 }
40 }
41 for (int comp = 0; comp < numComp; ++comp) {
42 integrals[comp] *= volume_product;
43 }
44 }
45
46 void Rectangle::integral_order3(std::vector<double>& integrals, const escript::Data& arg) const {
47 const double weights[] = {0.166666666667, 0.833333333333, 0.833333333333, 0.166666666667};
48 const int numComp = arg.getDataPointSize();
49 const double volume_product = 0.25*m_dx[0]*m_dx[1];
50 for (int ei = 0; ei < m_NE[1]; ++ei) {
51 for (int ej = 0; ej < m_NE[0]; ++ej) {
52 const double *e = arg.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
53 double result = 0;
54 for (int comp = 0; comp < numComp; ++comp) {
55 for (int i = 0; i < 4; ++i) {
56 for (int j = 0; j < 4; ++j) {
57 result += weights[i] * weights[j] * e[INDEX3(comp,i,j,numComp,4)];
58 }
59 }
60 integrals[comp] += result;
61 }
62 }
63 }
64 for (int comp = 0; comp < numComp; ++comp) {
65 integrals[comp] *= volume_product;
66 }
67 }
68
69 void Rectangle::integral_order4(std::vector<double>& integrals, const escript::Data& arg) const {
70 const double weights[] = {0.1, 0.544444444444, 0.711111111111, 0.544444444444, 0.1};
71 const int numComp = arg.getDataPointSize();
72 const double volume_product = 0.25*m_dx[0]*m_dx[1];
73 for (int ei = 0; ei < m_NE[1]; ++ei) {
74 for (int ej = 0; ej < m_NE[0]; ++ej) {
75 const double *e = arg.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
76 double result = 0;
77 for (int comp = 0; comp < numComp; ++comp) {
78 for (int i = 0; i < 5; ++i) {
79 for (int j = 0; j < 5; ++j) {
80 result += weights[i] * weights[j] * e[INDEX3(comp,i,j,numComp,5)];
81 }
82 }
83 integrals[comp] += result;
84 }
85 }
86 }
87 for (int comp = 0; comp < numComp; ++comp) {
88 integrals[comp] *= volume_product;
89 }
90 }
91
92 void Rectangle::integral_order5(std::vector<double>& integrals, const escript::Data& arg) const {
93 const double weights[] = {0.0666666666667, 0.378474956298, 0.554858377035, 0.554858377035, 0.378474956298, 0.0666666666667};
94 const int numComp = arg.getDataPointSize();
95 const double volume_product = 0.25*m_dx[0]*m_dx[1];
96 for (int ei = 0; ei < m_NE[1]; ++ei) {
97 for (int ej = 0; ej < m_NE[0]; ++ej) {
98 const double *e = arg.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
99 double result = 0;
100 for (int comp = 0; comp < numComp; ++comp) {
101 for (int i = 0; i < 6; ++i) {
102 for (int j = 0; j < 6; ++j) {
103 result += weights[i] * weights[j] * e[INDEX3(comp,i,j,numComp,6)];
104 }
105 }
106 integrals[comp] += result;
107 }
108 }
109 }
110 for (int comp = 0; comp < numComp; ++comp) {
111 integrals[comp] *= volume_product;
112 }
113 }
114
115 void Rectangle::integral_order6(std::vector<double>& integrals, const escript::Data& arg) const {
116 const double weights[] = {0.047619047619, 0.276826047362, 0.43174538121, 0.487619047619, 0.43174538121, 0.276826047362, 0.047619047619};
117 const int numComp = arg.getDataPointSize();
118 const double volume_product = 0.25*m_dx[0]*m_dx[1];
119 for (int ei = 0; ei < m_NE[1]; ++ei) {
120 for (int ej = 0; ej < m_NE[0]; ++ej) {
121 const double *e = arg.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
122 double result = 0;
123 for (int comp = 0; comp < numComp; ++comp) {
124 for (int i = 0; i < 7; ++i) {
125 for (int j = 0; j < 7; ++j) {
126 result += weights[i] * weights[j] * e[INDEX3(comp,i,j,numComp,7)];
127 }
128 }
129 integrals[comp] += result;
130 }
131 }
132 }
133 for (int comp = 0; comp < numComp; ++comp) {
134 integrals[comp] *= volume_product;
135 }
136 }
137
138 void Rectangle::integral_order7(std::vector<double>& integrals, const escript::Data& arg) const {
139 const double weights[] = {0.0357142857143, 0.210704227144, 0.341122692484, 0.412458794659, 0.412458794659, 0.341122692484, 0.210704227144, 0.0357142857143};
140 const int numComp = arg.getDataPointSize();
141 const double volume_product = 0.25*m_dx[0]*m_dx[1];
142 for (int ei = 0; ei < m_NE[1]; ++ei) {
143 for (int ej = 0; ej < m_NE[0]; ++ej) {
144 const double *e = arg.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
145 double result = 0;
146 for (int comp = 0; comp < numComp; ++comp) {
147 for (int i = 0; i < 8; ++i) {
148 for (int j = 0; j < 8; ++j) {
149 result += weights[i] * weights[j] * e[INDEX3(comp,i,j,numComp,8)];
150 }
151 }
152 integrals[comp] += result;
153 }
154 }
155 }
156 for (int comp = 0; comp < numComp; ++comp) {
157 integrals[comp] *= volume_product;
158 }
159 }
160
161 void Rectangle::integral_order8(std::vector<double>& integrals, const escript::Data& arg) const {
162 const double weights[] = {0.0277777777778, 0.165495361561, 0.2745387125, 0.346428510973, 0.371519274376, 0.346428510973, 0.2745387125, 0.165495361561, 0.0277777777778};
163 const int numComp = arg.getDataPointSize();
164 const double volume_product = 0.25*m_dx[0]*m_dx[1];
165 for (int ei = 0; ei < m_NE[1]; ++ei) {
166 for (int ej = 0; ej < m_NE[0]; ++ej) {
167 const double *e = arg.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
168 double result = 0;
169 for (int comp = 0; comp < numComp; ++comp) {
170 for (int i = 0; i < 9; ++i) {
171 for (int j = 0; j < 9; ++j) {
172 result += weights[i] * weights[j] * e[INDEX3(comp,i,j,numComp,9)];
173 }
174 }
175 integrals[comp] += result;
176 }
177 }
178 }
179 for (int comp = 0; comp < numComp; ++comp) {
180 integrals[comp] *= volume_product;
181 }
182 }
183
184 void Rectangle::integral_order9(std::vector<double>& integrals, const escript::Data& arg) const {
185 const double weights[] = {0.0222222222222, 0.133305990851, 0.224889342063, 0.29204268368, 0.327539761184, 0.327539761184, 0.29204268368, 0.224889342063, 0.133305990851, 0.0222222222222};
186 const int numComp = arg.getDataPointSize();
187 const double volume_product = 0.25*m_dx[0]*m_dx[1];
188 for (int ei = 0; ei < m_NE[1]; ++ei) {
189 for (int ej = 0; ej < m_NE[0]; ++ej) {
190 const double *e = arg.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
191 double result = 0;
192 for (int comp = 0; comp < numComp; ++comp) {
193 for (int i = 0; i < 10; ++i) {
194 for (int j = 0; j < 10; ++j) {
195 result += weights[i] * weights[j] * e[INDEX3(comp,i,j,numComp,10)];
196 }
197 }
198 integrals[comp] += result;
199 }
200 }
201 }
202 for (int comp = 0; comp < numComp; ++comp) {
203 integrals[comp] *= volume_product;
204 }
205 }
206
207 void Rectangle::integral_order10(std::vector<double>& integrals, const escript::Data& arg) const {
208 const double weights[] = {0.0181818181818, 0.109612273267, 0.18716988178, 0.248048104264, 0.286879124779, 0.300217595456, 0.286879124779, 0.248048104264, 0.18716988178, 0.109612273267, 0.0181818181818};
209 const int numComp = arg.getDataPointSize();
210 const double volume_product = 0.25*m_dx[0]*m_dx[1];
211 for (int ei = 0; ei < m_NE[1]; ++ei) {
212 for (int ej = 0; ej < m_NE[0]; ++ej) {
213 const double *e = arg.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
214 double result = 0;
215 for (int comp = 0; comp < numComp; ++comp) {
216 for (int i = 0; i < 11; ++i) {
217 for (int j = 0; j < 11; ++j) {
218 result += weights[i] * weights[j] * e[INDEX3(comp,i,j,numComp,11)];
219 }
220 }
221 integrals[comp] += result;
222 }
223 }
224 }
225 for (int comp = 0; comp < numComp; ++comp) {
226 integrals[comp] *= volume_product;
227 }
228 }
229
230 }

  ViewVC Help
Powered by ViewVC 1.1.26