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

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 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 #define ESNEEDPYTHON 17 #include "esysUtils/first.h" 18 19 #include 20 #include 21 22 namespace speckley { 23 void Rectangle::integral_order2(std::vector& 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& 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& 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& 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& 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& 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& 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& 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& 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 }