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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6112 - (show annotations)
Thu Mar 31 09:40:10 2016 UTC (3 years, 3 months ago) by jfenwick
File size: 26625 byte(s)
Relicense all the things!


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
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/WaveAssembler2D.h>
18 #include <speckley/domainhelpers.h>
19
20 #include <escript/index.h>
21
22 const double all_weights[][11] = {
23 {0.333333333333, 1.33333333333, 0.333333333333},
24 {0.166666666667, 0.833333333333, 0.833333333333, 0.166666666667},
25 {0.1, 0.544444444444, 0.711111111111, 0.544444444444, 0.1},
26 {0.0666666666667, 0.378474956298, 0.554858377035, 0.554858377035, 0.378474956298, 0.0666666666667},
27 {0.047619047619, 0.276826047362, 0.43174538121, 0.487619047619, 0.43174538121, 0.276826047362, 0.047619047619},
28 {0.0357142857143, 0.210704227144, 0.341122692484, 0.412458794659, 0.412458794659, 0.341122692484,0.210704227144, 0.0357142857143},
29 {0.0277777777778, 0.165495361561, 0.2745387125, 0.346428510973, 0.371519274376, 0.346428510973, 0.2745387125, 0.165495361561, 0.0277777777778},
30 {0.0222222222222, 0.133305990851, 0.224889342063, 0.29204268368, 0.327539761184, 0.327539761184, 0.29204268368, 0.224889342063, 0.133305990851, 0.0222222222222},
31 {0.0181818181818, 0.109612273267, 0.18716988178, 0.248048104264, 0.286879124779, 0.300217595456, 0.286879124779, 0.248048104264, 0.18716988178, 0.109612273267, 0.0181818181818}
32 };
33
34 const double all_lagrange_derivs[][11][11] = {
35 { // order 2
36 {-1.50000000000000, -0.500000000000000, 0.500000000000000},
37 {2.00000000000000, 0, -2.00000000000000},
38 {-0.500000000000000, 0.500000000000000, 1.50000000000000}
39 }, { //order 3
40 {-3.00000000000000, -0.809016994374948, 0.309016994374948, -0.500000000000000},
41 {4.04508497187474, 4.44089209850063e-16, -1.11803398874990, 1.54508497187474},
42 {-1.54508497187474, 1.11803398874989, 2.22044604925031e-16, -4.04508497187474},
43 {0.500000000000000, -0.309016994374947, 0.809016994374948, 3.00000000000000}
44 }, { //order 4
45 {-4.99999999999999, -1.24099025303098, 0.374999999999999, -0.259009746969017, 0.499999999999999},
46 {6.75650248872424, -6.66133814775094e-15, -1.33658457769545, 0.763762615825974, -1.41016417794243},
47 {-2.66666666666667, 1.74574312188794, 1.44328993201270e-15, -1.74574312188794, 2.66666666666667},
48 {1.41016417794243, -0.763762615825974, 1.33658457769545, 1.66533453693773e-15, -6.75650248872424},
49 {-0.500000000000001, 0.259009746969017, -0.375000000000000, 1.24099025303098, 5.00000000000000}
50 }, { //order 5
51 {-7.50000000000002, -1.78636494833911, 0.484951047853572, -0.269700610832040, 0.237781177984232, -0.500000000000002},
52 {10.1414159363197, 2.13162820728030e-14, -1.72125695283023, 0.786356672223240, -0.653547507429800, 1.34991331419049},
53 {-4.03618727030532, 2.52342677742945, -4.66293670342566e-15, -1.75296196636786, 1.15282815853593, -2.24468464817616},
54 {2.24468464817616, -1.15282815853593, 1.75296196636787, -1.77635683940025e-15, -2.52342677742946, 4.03618727030535},
55 {-1.34991331419048, 0.653547507429800, -0.786356672223242, 1.72125695283023, 2.22044604925031e-15, -10.1414159363197},
56 {0.499999999999998, -0.237781177984231, 0.269700610832039, -0.484951047853569, 1.78636494833909, 7.50000000000000}
57 }, { //order 6
58 {-10.5000000000000, -2.44292601424426, 0.625256665515336, -0.312499999999997, 0.226099400942572, -0.226611870395444, 0.500000000000001},
59 {14.2015766029198, -4.17443857259059e-14, -2.21580428316997, 0.907544471268819, -0.616390835517577, 0.602247179635785, -1.31737343570244},
60 {-5.66898522554555, 3.45582821429430, 3.10862446895044e-15, -2.00696924058875, 1.06644190400637, -0.961339797288714, 2.04996481307676},
61 {3.20000000000003, -1.59860668809837, 2.26669808708600, 1.33226762955019e-15, -2.26669808708599, 1.59860668809837, -3.20000000000003},
62 {-2.04996481307676, 0.961339797288717, -1.06644190400638, 2.00696924058876, -1.77635683940025e-14, -3.45582821429431, 5.66898522554558},
63 {1.31737343570245, -0.602247179635788, 0.616390835517580, -0.907544471268822, 2.21580428316998, 6.76125821996720e-14, -14.2015766029198},
64 {-0.500000000000000, 0.226611870395444, -0.226099400942572, 0.312499999999997, -0.625256665515335, 2.44292601424425, 10.5000000000000}
65 }, { //order 7
66 {-13.9999999999999, -3.20991570300295, 0.792476681320508, -0.372150435728592, 0.243330712723790, -0.203284568900591, 0.219957514771299, -0.499999999999980},
67 {18.9375986071174, -6.23945339839338e-14, -2.80647579473643, 1.07894468879045, -0.661157350900312, 0.537039586157660, -0.573565414940254, 1.29768738832019},
68 {-7.56928981934855, 4.54358506456659, 1.17683640610267e-14, -2.37818723351551, 1.13535801688112, -0.845022556506511, 0.869448098331479, -1.94165942554406},
69 {4.29790816426521, -2.11206121431455, 2.87551740597250, 3.77475828372553e-15, -2.38892435915824, 1.37278583180603, -1.29423205091348, 2.81018898925786},
70 {-2.81018898925796, 1.29423205091350, -1.37278583180602, 2.38892435915823, 7.77156117237610e-15, -2.87551740597249, 2.11206121431450, -4.29790816426503},
71 {1.94165942554413, -0.869448098331493, 0.845022556506508, -1.13535801688111, 2.37818723351550, -2.44249065417534e-14, -4.54358506456649, 7.56928981934828},
72 {-1.29768738832026, 0.573565414940274, -0.537039586157669, 0.661157350900321, -1.07894468879047, 2.80647579473648, -1.71085368094737e-13, -18.9375986071174},
73 {0.500000000000021, -0.219957514771312, 0.203284568900599, -0.243330712723800, 0.372150435728609, -0.792476681320546, 3.20991570300311, 14.0000000000002}
74 }, { //order 8
75 {-18.0000000000010, -4.08701370203454, 0.985360090074639, -0.444613449281139, 0.273437500000029, -0.207734512035617, 0.189655591978376, -0.215654018702531, 0.500000000000095},
76 {24.3497451715930, 1.34892097491957e-12, -3.48835875343438, 1.28796075006388, -0.741782397916244, 0.547300160534042, -0.492350938315503, 0.555704981283736, -1.28483063269969},
77 {-9.73870165721010, 5.78680581663678, -3.21298543326520e-13, -2.83445891207935, 1.26941308635811, -0.855726185092640, 0.738349277190360, -0.816756381741392, 1.87444087344708},
78 {5.54496390694879, -2.69606544031400, 3.57668094012577, -3.10862446895044e-15, -2.65931021757391, 1.37696489376050, -1.07980381128263, 1.14565373845518, -2.59074567655957},
79 {-3.65714285714248, 1.66522164500537, -1.71783215719513, 2.85191596846290, -1.06581410364015e-14, -2.85191596846287, 1.71783215719506, -1.66522164500546, 3.65714285714316},
80 {2.59074567655910, -1.14565373845513, 1.07980381128268, -1.37696489376052, 2.65931021757393, -3.71924713249427e-14, -3.57668094012566, 2.69606544031420, -5.54496390694988},
81 {-1.87444087344680, 0.816756381741386, -0.738349277190418, 0.855726185092682, -1.26941308635816, 2.83445891207943, 9.27036225562006e-14, -5.78680581663756, 9.73870165721225},
82 {1.28483063269940, -0.555704981283693, 0.492350938315507, -0.547300160534031, 0.741782397916225, -1.28796075006385, 3.48835875343430, 5.71542813077031e-13, -24.3497451715928},
83 {-0.499999999999903, 0.215654018702479, -0.189655591978347, 0.207734512035579, -0.273437499999975, 0.444613449281046, -0.985360090074401, 4.08701370203326, 17.9999999999994}
84 }, { //order 9
85 {-22.4999999999988, -5.07406470297709, 1.20335199285206, -0.528369376820220, 0.312047255608382, -0.223527944742433, 0.186645789393719, -0.180786585489230, 0.212702758009187, -0.500000000000077},
86 {30.4381450292820, -1.52677870346452e-12, -4.25929735496529, 1.52990263818163, -0.845813573406436, 0.588082143045176, -0.483462326333953, 0.464274958908154, -0.543753738235757, 1.27595483609299},
87 {-12.1779467074315, 7.18550286970643, 3.27293747659496e-13, -3.36412586829791, 1.44485031560171, -0.916555180336469, 0.721237312721631, -0.676797087196100, 0.783239293138005, -1.82956393190377},
88 {6.94378848513465, -3.35166386274684, 4.36867455701003, 3.17523785042795e-14, -3.02021795819936, 1.46805550939000, -1.04618936550250, 0.936603213139437, -1.05915446364554, 2.45288417544331},
89 {-4.59935476110357, 2.07820799403643, -2.10435017941307, 3.38731810120242, 2.22044604925031e-16, -3.02518848775198, 1.64649408398706, -1.33491548387823, 1.44494844875159, -3.29464303375000},
90 {3.29464303374949, -1.44494844875146, 1.33491548387820, -1.64649408398705, 3.02518848775197, 6.66133814775094e-15, -3.38731810120245, 2.10435017941312, -2.07820799403661, 4.59935476110425},
91 {-2.45288417544291, 1.05915446364544, -0.936603213139410, 1.04618936550249, -1.46805550938999, 3.02021795819934, -1.55431223447522e-15, -4.36867455701010, 3.35166386274711, -6.94378848513561},
92 {1.82956393190345, -0.783239293137921, 0.676797087196070, -0.721237312721610, 0.916555180336448, -1.44485031560168, 3.36412586829786, -2.10609307771392e-13, -7.18550286970689, 12.1779467074328},
93 {-1.27595483609268, 0.543753738235664, -0.464274958908104, 0.483462326333909, -0.588082143045126, 0.845813573406365, -1.52990263818151, 4.25929735496501, 2.64499533386697e-12, -30.4381450292815},
94 {0.499999999999919, -0.212702758009135, 0.180786585489197, -0.186645789393687, 0.223527944742396, -0.312047255608330, 0.528369376820133, -1.20335199285186, 5.07406470297626, 22.4999999999976}
95 }, { //order 10
96 {-27.4999999999896, -6.17098569730879, 1.44617248279108, -0.622725214738251, 0.357476373116252, -0.246093749999837, 0.194287668796289, -0.172970108511720, 0.174657862947473, -0.210587346312973, 0.500000000000200},
97 {37.2028673819635, -1.45093936865237e-11, -5.11821182477732, 1.80264679987985, -0.968487138802308, 0.646939963832195, -0.502643313012917, 0.443395636437341, -0.445313527290911, 0.535331085929629, -1.26956267628665},
98 {-14.8873962951139, 8.73967005352577, 4.30699920173083e-12, -3.96199657704633, 1.65273663422738, -1.00650540857701, 0.747734824688410, -0.643586209360019, 0.637362056418227, -0.760400982244270, 1.79798803579712},
99 {8.49561949463899, -4.07931619371279, 5.25066175545241, -4.06341627012807e-13, -3.45061471617355, 1.60812790372552, -1.07998725049569, 0.884587314556481, -0.852916813558380, 1.00338624297419, -2.35976991309107},
100 {-5.64038799768972, 2.53474117868649, -2.53318340860873, 3.99079578802906, 9.19264664389630e-14, -3.30517685337838, 1.69057057046883, -1.24905529156928, 1.14606853428002, -1.31552671443958, 3.06553920088515},
101 {4.06349206349476, -1.77190705526895, 1.61441910793459, -1.94634945457108, 3.45885134807703, 3.86357612569554e-14, -3.45885134807708, 1.94634945457117, -1.61441910793481, 1.77190705526946, -4.06349206349634},
102 {-3.06553920088389, 1.31552671443917, -1.14606853427984, 1.24905529156920, -1.69057057046877, 3.30517685337831, -4.28546087505310e-14, -3.99079578802918, 2.53318340860904, -2.53474117868717, 5.64038799769177},
103 {2.35976991309003, -1.00338624297385, 0.852916813558223, -0.884587314556402, 1.07998725049562, -1.60812790372545, 3.45061471617347, 5.73430192218893e-13, -5.25066175545293, 4.07931619371373, -8.49561949464169},
104 {-1.79798803579616, 0.760400982243941, -0.637362056418051, 0.643586209359903, -0.747734824688294, 1.00650540857687, -1.65273663422718, 3.96199657704598, -3.31623617455534e-12, -8.73967005352657, 14.8873962951164},
105 {1.26956267628575, -0.535331085929307, 0.445313527290713, -0.443395636437185, 0.502643313012754, -0.646939963831994, 0.968487138802021, -1.80264679987934, 5.11821182477613, 1.64738223062955e-11, -37.2028673819613},
106 {-0.499999999999790, 0.210587346312823, -0.174657862947375, 0.172970108511639, -0.194287668796203, 0.246093749999731, -0.357476373116102, 0.622725214737994, -1.44617248279054, 6.17098569730708, 27.4999999999864}
107 }
108 };
109
110 using escript::AbstractSystemMatrix;
111 using escript::Data;
112
113 namespace speckley {
114
115 void WaveAssembler2D::collateFunctionSpaceTypes(std::vector<int>& fsTypes,
116 const DataMap& coefs) const
117 {
118 if (isNotEmpty("D", coefs))
119 fsTypes.push_back(coefs.find("D")->second.getFunctionSpace().getTypeCode());
120 if (isNotEmpty("Y", coefs))
121 fsTypes.push_back(coefs.find("Y")->second.getFunctionSpace().getTypeCode());
122 if (isNotEmpty("du", coefs))
123 fsTypes.push_back(coefs.find("du")->second.getFunctionSpace().getTypeCode());
124 }
125
126 void WaveAssembler2D::assemblePDESingle(AbstractSystemMatrix* mat,
127 Data& rhs, const DataMap& coefs) const
128 {
129 const Data& A = unpackData("A", coefs);
130 const Data& B = unpackData("B", coefs);
131 const Data& C = unpackData("C", coefs);
132 const Data& D = unpackData("D", coefs);
133 const Data& du = unpackData("du", coefs);
134 const Data& Y = unpackData("Y", coefs);
135 assemblePDESingle(mat, rhs, A, B, C, D, du, Y);
136
137 }
138
139 void WaveAssembler2D::assemblePDEBoundarySingle(AbstractSystemMatrix* mat,
140 Data& rhs, const DataMap& coefs) const
141 {
142 const Data& d = unpackData("d", coefs);
143 const Data& y = unpackData("y", coefs);
144 assemblePDEBoundarySingle(mat, rhs, d, y);
145 }
146
147 void WaveAssembler2D::assemblePDESingleReduced(AbstractSystemMatrix* mat,
148 Data& rhs, const DataMap& coefs) const
149 {
150 const Data& A = unpackData("A", coefs);
151 const Data& B = unpackData("B", coefs);
152 const Data& C = unpackData("C", coefs);
153 const Data& D = unpackData("D", coefs);
154 const Data& du = unpackData("du", coefs);
155 const Data& Y = unpackData("Y", coefs);
156 assemblePDESingleReduced(mat, rhs, A, B, C, D, du, Y);
157 }
158
159 void WaveAssembler2D::assemblePDEBoundarySingleReduced(
160 AbstractSystemMatrix* mat, Data& rhs,
161 const DataMap& coefs) const
162 {
163 const Data& d = unpackData("d", coefs);
164 const Data& y = unpackData("y", coefs);
165 assemblePDEBoundarySingleReduced(mat, rhs, d, y);
166 }
167
168 void WaveAssembler2D::assemblePDESystem(AbstractSystemMatrix* mat,
169 Data& rhs, const DataMap& coefs) const
170 {
171 if (!unpackData("X", coefs).isEmpty())
172 throw SpeckleyException("Wave assembler does not support X");
173 const Data& A = unpackData("A", coefs);
174 const Data& B = unpackData("B", coefs);
175 const Data& C = unpackData("C", coefs);
176 const Data& D = unpackData("D", coefs);
177 const Data& du = unpackData("du", coefs);
178 const Data& Y = unpackData("Y", coefs);
179 assemblePDESystem(mat, rhs, A, B, C, D, du, Y);
180 }
181
182 void WaveAssembler2D::assemblePDEBoundarySystem(AbstractSystemMatrix* mat,
183 Data& rhs, const DataMap& coefs) const
184 {
185 const Data& d = unpackData("d", coefs);
186 const Data& y = unpackData("y", coefs);
187 assemblePDEBoundarySystem(mat, rhs, d, y);
188 }
189
190 void WaveAssembler2D::assemblePDESystemReduced(AbstractSystemMatrix* mat,
191 Data& rhs, const DataMap& coefs) const
192 {
193 const Data& A = unpackData("A", coefs);
194 const Data& B = unpackData("B", coefs);
195 const Data& C = unpackData("C", coefs);
196 const Data& D = unpackData("D", coefs);
197 const Data& du = unpackData("du", coefs);
198 const Data& Y = unpackData("Y", coefs);
199 assemblePDESystemReduced(mat, rhs, A, B, C, D, du, Y);
200 }
201
202 void WaveAssembler2D::assemblePDEBoundarySystemReduced(
203 AbstractSystemMatrix* mat, Data& rhs,
204 const DataMap& coefs) const
205 {
206 const Data& d = unpackData("d", coefs);
207 const Data& y = unpackData("y", coefs);
208 assemblePDEBoundarySystemReduced(mat, rhs, d, y);
209 }
210
211 void WaveAssembler2D::assemblePDESystem(AbstractSystemMatrix* mat,
212 Data& rhs, const Data& A, const Data& B,
213 const Data& C, const Data& D,
214 const Data& du, const Data& Y) const
215 {
216 if (!(A.isEmpty() && B.isEmpty() && C.isEmpty()))
217 throw SpeckleyException("Speckley does not support PDEs using A, B or C");
218 int order = domain->m_order;
219 const double *weights = all_weights[order-2];
220 const double volume_product = m_dx[0]*m_dx[1]/4.;
221 const int NE0 = m_NE[0];
222 const int NE1 = m_NE[1];
223 const int quads = order + 1;
224 const int max_x = m_NN[0];
225 dim_t numComp;
226 if (!mat)
227 numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
228 else {
229 numComp=mat->getColumnBlockSize();
230 }
231
232 rhs.requireWrite();
233 int size = 1;
234 if (!Y.isEmpty()) {
235 size = Y.getDataPointSize();
236 }
237 const index_t y_indices[2] = {0,size-1};
238 size = 1;
239 if (!D.isEmpty()) {
240 size = D.getDataPointSize();
241 }
242 const index_t d_indices[2] = {0,size-1};
243
244 if (!D.isEmpty() && (!du.isEmpty() || !Y.isEmpty()))
245 throw SpeckleyException("Speckley does not support adding left and right sides concurrently");
246
247 for (dim_t colouring = 0; colouring < 2; colouring++) {
248 #pragma omp parallel for
249 for (dim_t ey = colouring; ey < NE1; ey += 2) {
250 for (dim_t ex = 0; ex < NE0; ex++) {
251 const index_t e_index = ex + ey*NE0; //element number for Elements
252 const index_t start = ex*order + ey*max_x*order; //node start for Nodes
253
254 if (!D.isEmpty()) {
255 const double *e = D.getSampleDataRO(e_index);
256 if (D.actsExpanded()) {
257 for (index_t comp = 0; comp < numComp; comp++) {
258 for (short qy = 0; qy < quads; qy++) {
259 for (short qx = 0; qx < quads; qx++) {
260 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
261 /* D */ out[comp] += volume_product * weights[qx] * weights[qy]
262 * e[INDEX3(comp, qx, qy, numComp, quads)]; //D(i),qx,qy
263 }
264 }
265 }
266 } else { //constant
267 for (index_t comp = 0; comp < numComp; comp++) {
268 for (short qy = 0; qy < quads; qy++) {
269 for (short qx = 0; qx < quads; qx++) {
270 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
271 /* D */ out[comp] += volume_product * weights[qx] * weights[qy]
272 * e[d_indices[comp]]; //D(i)
273 }
274 }
275 }
276 }
277 }
278
279 if (!du.isEmpty()) {
280 const double *du_p = du.getSampleDataRO(e_index);
281 const double c11_v = -c11.getSampleDataRO(e_index)[0];
282 const double c13_v = -c13.getSampleDataRO(e_index)[0];
283 const double c33_v = -c33.getSampleDataRO(e_index)[0];
284 const double diagonal = (
285 isHTI ? -c66.getSampleDataRO(e_index)[0]
286 : -c44.getSampleDataRO(e_index)[0]);
287
288
289 if (du.actsExpanded()) {
290 for (short qy = 0; qy < quads; qy++) { //a
291 for (short qx = 0; qx < quads; qx++) { //b
292 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
293 double res_a = 0, res_b = 0; //comp 0
294 double res_c = 0, res_d = 0; //comp 1
295 for (short k = 0; k < quads; k++) {
296 #define SI(_x_,_y_,_qx_,_qy_) INDEX4((_x_),(_y_),(_qx_),(_qy_),numComp,2,quads)
297 res_a += weights[k] * all_lagrange_derivs[order-2][qx][k]
298 * (c11_v*du_p[SI(0,0,k,qy)] + c13_v*du_p[SI(1,1,k,qy)]);
299 res_b += weights[k] * all_lagrange_derivs[order-2][qy][k]
300 * diagonal*( du_p[SI(0,1,qx,k)] + du_p[SI(1,0,qx,k)] );
301
302 res_c += weights[k] * all_lagrange_derivs[order-2][qx][k]
303 * diagonal*(du_p[SI(0,1,k,qy)] + du_p[SI(1,0,k,qy)]);
304 res_d += weights[k] * all_lagrange_derivs[order-2][qy][k]
305 * (c13_v*du_p[SI(0,0,qx,k)] + c33_v*du_p[SI(1,1,qx,k)]);
306 #undef SI
307 }
308 out[0] += 2 * volume_product
309 * (res_a * weights[qy] / m_dx[0] + res_b * weights[qx] / m_dx[1]);
310 out[1] += 2 * volume_product
311 * (res_c * weights[qy] / m_dx[0] + res_d * weights[qx] / m_dx[1]);
312 }
313 }
314 } else { //constant
315 for (short qy = 0; qy < quads; qy++) {
316 for (short qx = 0; qx < quads; qx++) {
317 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
318 double res_a = 0, res_b = 0; //comp 0
319 double res_c = 0, res_d = 0; //comp 1
320 for (short k = 0; k < quads; k++) {
321 res_a += weights[k] * all_lagrange_derivs[order-2][qx][k]
322 * (c11_v*du_p[0] + c13_v*du_p[2]);
323 res_b += weights[k] * all_lagrange_derivs[order-2][qy][k]
324 * diagonal*(du_p[1] + du_p[2]);
325 res_c += weights[k] * all_lagrange_derivs[order-2][qx][k]
326 * diagonal*(du_p[1] + du_p[2]);
327 res_d += weights[k] * all_lagrange_derivs[order-2][qy][k]
328 * (c13_v*du_p[0] + c33_v*du_p[3]);
329 }
330 out[0] += 2 * volume_product
331 * (res_a * weights[qy] / m_dx[0] + res_b * weights[qx] / m_dx[1]);
332 out[1] += 2 * volume_product
333 * (res_c * weights[qy] / m_dx[0] + res_d * weights[qx] / m_dx[1]);
334 }
335 }
336 }
337 }
338
339 if (!Y.isEmpty()) {
340 const double *e = Y.getSampleDataRO(e_index);
341 if (Y.actsExpanded()) {
342 for (index_t comp = 0; comp < numComp; comp++) {
343 for (short qy = 0; qy < quads; qy++) {
344 for (short qx = 0; qx < quads; qx++) {
345 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
346 /* Y */ out[comp] += volume_product * weights[qx] * weights[qy]
347 * e[INDEX3(comp, qx, qy, numComp, quads)]; //Y(i),qx,qy
348 }
349 }
350 }
351 } else { //constant
352 for (index_t comp = 0; comp < numComp; comp++) {
353 for (short qy = 0; qy < quads; qy++) {
354 for (short qx = 0; qx < quads; qx++) {
355 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
356 /* Y */ out[comp] += volume_product * weights[qx] * weights[qy]
357 * e[y_indices[comp]]; //Y(i)
358 }
359 }
360 }
361 }
362 }
363 }
364 }
365 }
366 }
367
368 void WaveAssembler2D::assemblePDESystemReduced(AbstractSystemMatrix* mat,
369 Data& rhs, const Data& A, const Data& B,
370 const Data& C, const Data& D,
371 const Data& du, const Data& Y) const
372 {
373 throw SpeckleyException("Speckley does not support reduced functionspaces");
374
375 }
376
377 void WaveAssembler2D::assemblePDEBoundarySingle(AbstractSystemMatrix* mat,
378 Data& rhs, const Data& d, const Data& y) const
379 {
380 throw SpeckleyException("Speckley does not support boundary functionspaces");
381 }
382
383 void WaveAssembler2D::assemblePDEBoundarySingleReduced(
384 AbstractSystemMatrix* mat, Data& rhs,
385 const Data& d, const Data& y) const
386 {
387 throw SpeckleyException("Speckley does not support reduced functionspaces");
388 }
389
390 void WaveAssembler2D::assemblePDEBoundarySystem(AbstractSystemMatrix* mat,
391 Data& rhs, const Data& d, const Data& y) const
392 {
393 throw SpeckleyException("Speckley does not support boundary functionspaces");
394 }
395
396 void WaveAssembler2D::assemblePDEBoundarySystemReduced(
397 AbstractSystemMatrix* mat, Data& rhs,
398 const Data& d, const Data& y) const
399 {
400 throw SpeckleyException("Speckley does not support reduced functionspaces");
401 }
402
403 //protected
404 void WaveAssembler2D::assemblePDESingle(AbstractSystemMatrix *mat,
405 escript::Data& rhs, const escript::Data& A, const escript::Data& B,
406 const escript::Data& C, const escript::Data& D,
407 const escript::Data& du, const escript::Data& Y) const
408 {
409 throw SpeckleyException("WaveAssembler does not support assemblePDESingle()");
410 }
411
412 //protected
413 void WaveAssembler2D::assemblePDESingleReduced(AbstractSystemMatrix* mat,
414 Data& rhs, const Data& A, const Data& B,
415 const Data& C, const Data& D,
416 const Data& du, const Data& Y) const
417 {
418 throw SpeckleyException("Speckley does not support reduced functionspaces");
419 }
420
421
422 } // namespace speckley
423

  ViewVC Help
Powered by ViewVC 1.1.26