# Contents of /trunk/speckley/src/RectangleReductions.cpp

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 6 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 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 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 }