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

Contents of /trunk/speckley/src/RectangleReductions.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: 8848 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
17 #include <speckley/Rectangle.h>
18
19 namespace speckley {
20 void Rectangle::reduction_order2(const escript::Data& in, escript::Data& out) const {
21 const double weights[] = {0.333333333333, 1.33333333333, 0.333333333333};
22 const int numComp = in.getDataPointSize();
23 for (int ei = 0; ei < m_NE[1]; ++ei) {
24 for (int ej = 0; ej < m_NE[0]; ++ej) {
25 const double *e_in = in.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
26 double *e_out = out.getSampleDataRW(INDEX2(ej,ei,m_NE[0]));
27 for (int comp = 0; comp < numComp; ++comp) {
28 double result = 0;
29 for (int i = 0; i < 3; ++i) {
30 for (int j = 0; j < 3; ++j) {
31 result += weights[i] * weights[j] * e_in[INDEX3(comp,j,i,numComp,3)];
32 }
33 }
34 e_out[comp] += result / 4.;
35 }
36 }
37 }
38 }
39
40 void Rectangle::reduction_order3(const escript::Data& in, escript::Data& out) const {
41 const double weights[] = {0.166666666667, 0.833333333333, 0.833333333333, 0.166666666667};
42 const int numComp = in.getDataPointSize();
43 for (int ei = 0; ei < m_NE[1]; ++ei) {
44 for (int ej = 0; ej < m_NE[0]; ++ej) {
45 const double *e_in = in.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
46 double *e_out = out.getSampleDataRW(INDEX2(ej,ei,m_NE[0]));
47 for (int comp = 0; comp < numComp; ++comp) {
48 double result = 0;
49 for (int i = 0; i < 4; ++i) {
50 for (int j = 0; j < 4; ++j) {
51 result += weights[i] * weights[j] * e_in[INDEX3(comp,j,i,numComp,4)];
52 }
53 }
54 e_out[comp] += result / 4.;
55 }
56 }
57 }
58 }
59
60 void Rectangle::reduction_order4(const escript::Data& in, escript::Data& out) const {
61 const double weights[] = {0.1, 0.544444444444, 0.711111111111, 0.544444444444, 0.1};
62 const int numComp = in.getDataPointSize();
63 for (int ei = 0; ei < m_NE[1]; ++ei) {
64 for (int ej = 0; ej < m_NE[0]; ++ej) {
65 const double *e_in = in.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
66 double *e_out = out.getSampleDataRW(INDEX2(ej,ei,m_NE[0]));
67 for (int comp = 0; comp < numComp; ++comp) {
68 double result = 0;
69 for (int i = 0; i < 5; ++i) {
70 for (int j = 0; j < 5; ++j) {
71 result += weights[i] * weights[j] * e_in[INDEX3(comp,j,i,numComp,5)];
72 }
73 }
74 e_out[comp] += result / 4.;
75 }
76 }
77 }
78 }
79
80 void Rectangle::reduction_order5(const escript::Data& in, escript::Data& out) const {
81 const double weights[] = {0.0666666666667, 0.378474956298, 0.554858377035, 0.554858377035, 0.378474956298, 0.0666666666667};
82 const int numComp = in.getDataPointSize();
83 for (int ei = 0; ei < m_NE[1]; ++ei) {
84 for (int ej = 0; ej < m_NE[0]; ++ej) {
85 const double *e_in = in.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
86 double *e_out = out.getSampleDataRW(INDEX2(ej,ei,m_NE[0]));
87 for (int comp = 0; comp < numComp; ++comp) {
88 double result = 0;
89 for (int i = 0; i < 6; ++i) {
90 for (int j = 0; j < 6; ++j) {
91 result += weights[i] * weights[j] * e_in[INDEX3(comp,j,i,numComp,6)];
92 }
93 }
94 e_out[comp] += result / 4.;
95 }
96 }
97 }
98 }
99
100 void Rectangle::reduction_order6(const escript::Data& in, escript::Data& out) const {
101 const double weights[] = {0.047619047619, 0.276826047362, 0.43174538121, 0.487619047619, 0.43174538121, 0.276826047362, 0.047619047619};
102 const int numComp = in.getDataPointSize();
103 for (int ei = 0; ei < m_NE[1]; ++ei) {
104 for (int ej = 0; ej < m_NE[0]; ++ej) {
105 const double *e_in = in.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
106 double *e_out = out.getSampleDataRW(INDEX2(ej,ei,m_NE[0]));
107 for (int comp = 0; comp < numComp; ++comp) {
108 double result = 0;
109 for (int i = 0; i < 7; ++i) {
110 for (int j = 0; j < 7; ++j) {
111 result += weights[i] * weights[j] * e_in[INDEX3(comp,j,i,numComp,7)];
112 }
113 }
114 e_out[comp] += result / 4.;
115 }
116 }
117 }
118 }
119
120 void Rectangle::reduction_order7(const escript::Data& in, escript::Data& out) const {
121 const double weights[] = {0.0357142857143, 0.210704227144, 0.341122692484, 0.412458794659, 0.412458794659, 0.341122692484, 0.210704227144, 0.0357142857143};
122 const int numComp = in.getDataPointSize();
123 for (int ei = 0; ei < m_NE[1]; ++ei) {
124 for (int ej = 0; ej < m_NE[0]; ++ej) {
125 const double *e_in = in.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
126 double *e_out = out.getSampleDataRW(INDEX2(ej,ei,m_NE[0]));
127 for (int comp = 0; comp < numComp; ++comp) {
128 double result = 0;
129 for (int i = 0; i < 8; ++i) {
130 for (int j = 0; j < 8; ++j) {
131 result += weights[i] * weights[j] * e_in[INDEX3(comp,j,i,numComp,8)];
132 }
133 }
134 e_out[comp] += result / 4.;
135 }
136 }
137 }
138 }
139
140 void Rectangle::reduction_order8(const escript::Data& in, escript::Data& out) const {
141 const double weights[] = {0.0277777777778, 0.165495361561, 0.2745387125, 0.346428510973, 0.371519274376, 0.346428510973, 0.2745387125, 0.165495361561, 0.0277777777778};
142 const int numComp = in.getDataPointSize();
143 for (int ei = 0; ei < m_NE[1]; ++ei) {
144 for (int ej = 0; ej < m_NE[0]; ++ej) {
145 const double *e_in = in.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
146 double *e_out = out.getSampleDataRW(INDEX2(ej,ei,m_NE[0]));
147 for (int comp = 0; comp < numComp; ++comp) {
148 double result = 0;
149 for (int i = 0; i < 9; ++i) {
150 for (int j = 0; j < 9; ++j) {
151 result += weights[i] * weights[j] * e_in[INDEX3(comp,j,i,numComp,9)];
152 }
153 }
154 e_out[comp] += result / 4.;
155 }
156 }
157 }
158 }
159
160 void Rectangle::reduction_order9(const escript::Data& in, escript::Data& out) const {
161 const double weights[] = {0.0222222222222, 0.133305990851, 0.224889342063, 0.29204268368, 0.327539761184, 0.327539761184, 0.29204268368, 0.224889342063, 0.133305990851, 0.0222222222222};
162 const int numComp = in.getDataPointSize();
163 for (int ei = 0; ei < m_NE[1]; ++ei) {
164 for (int ej = 0; ej < m_NE[0]; ++ej) {
165 const double *e_in = in.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
166 double *e_out = out.getSampleDataRW(INDEX2(ej,ei,m_NE[0]));
167 for (int comp = 0; comp < numComp; ++comp) {
168 double result = 0;
169 for (int i = 0; i < 10; ++i) {
170 for (int j = 0; j < 10; ++j) {
171 result += weights[i] * weights[j] * e_in[INDEX3(comp,j,i,numComp,10)];
172 }
173 }
174 e_out[comp] += result / 4.;
175 }
176 }
177 }
178 }
179
180 void Rectangle::reduction_order10(const escript::Data& in, escript::Data& out) const {
181 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};
182 const int numComp = in.getDataPointSize();
183 for (int ei = 0; ei < m_NE[1]; ++ei) {
184 for (int ej = 0; ej < m_NE[0]; ++ej) {
185 const double *e_in = in.getSampleDataRO(INDEX2(ej,ei,m_NE[0]));
186 double *e_out = out.getSampleDataRW(INDEX2(ej,ei,m_NE[0]));
187 for (int comp = 0; comp < numComp; ++comp) {
188 double result = 0;
189 for (int i = 0; i < 11; ++i) {
190 for (int j = 0; j < 11; ++j) {
191 result += weights[i] * weights[j] * e_in[INDEX3(comp,j,i,numComp,11)];
192 }
193 }
194 e_out[comp] += result / 4.;
195 }
196 }
197 }
198 }
199
200 }

  ViewVC Help
Powered by ViewVC 1.1.26