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

Contents of /trunk/speckley/src/DefaultAssembler2D.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: 30223 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/DefaultAssembler2D.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 DefaultAssembler2D::collateFunctionSpaceTypes(std::vector<int>& fsTypes,
116 const DataMap& coefs) const
117 {
118 if (isNotEmpty("A", coefs))
119 fsTypes.push_back(coefs.find("A")->second.getFunctionSpace().getTypeCode());
120 if (isNotEmpty("B", coefs))
121 fsTypes.push_back(coefs.find("B")->second.getFunctionSpace().getTypeCode());
122 if (isNotEmpty("C", coefs))
123 fsTypes.push_back(coefs.find("C")->second.getFunctionSpace().getTypeCode());
124 if (isNotEmpty("D", coefs))
125 fsTypes.push_back(coefs.find("D")->second.getFunctionSpace().getTypeCode());
126 if (isNotEmpty("X", coefs))
127 fsTypes.push_back(coefs.find("X")->second.getFunctionSpace().getTypeCode());
128 if (isNotEmpty("Y", coefs))
129 fsTypes.push_back(coefs.find("Y")->second.getFunctionSpace().getTypeCode());
130 }
131
132 void DefaultAssembler2D::assemblePDESingle(AbstractSystemMatrix* mat,
133 Data& rhs, const DataMap& coefs) const
134 {
135 const Data& A = unpackData("A", coefs);
136 const Data& B = unpackData("B", coefs);
137 const Data& C = unpackData("C", coefs);
138 const Data& D = unpackData("D", coefs);
139 const Data& X = unpackData("X", coefs);
140 const Data& Y = unpackData("Y", coefs);
141 assemblePDESingle(mat, rhs, A, B, C, D, X, Y);
142
143 }
144
145 void DefaultAssembler2D::assemblePDEBoundarySingle(AbstractSystemMatrix* mat,
146 Data& rhs, const DataMap& coefs) const
147 {
148 const Data& d = unpackData("d", coefs);
149 const Data& y = unpackData("y", coefs);
150 assemblePDEBoundarySingle(mat, rhs, d, y);
151 }
152
153 void DefaultAssembler2D::assemblePDESingleReduced(AbstractSystemMatrix* mat,
154 Data& rhs, const DataMap& coefs) const
155 {
156 const Data& A = unpackData("A", coefs);
157 const Data& B = unpackData("B", coefs);
158 const Data& C = unpackData("C", coefs);
159 const Data& D = unpackData("D", coefs);
160 const Data& X = unpackData("X", coefs);
161 const Data& Y = unpackData("Y", coefs);
162 assemblePDESingleReduced(mat, rhs, A, B, C, D, X, Y);
163 }
164
165 void DefaultAssembler2D::assemblePDEBoundarySingleReduced(
166 AbstractSystemMatrix* mat, Data& rhs,
167 const DataMap& coefs) const
168 {
169 const Data& d = unpackData("d", coefs);
170 const Data& y = unpackData("y", coefs);
171 assemblePDEBoundarySingleReduced(mat, rhs, d, y);
172 }
173
174 void DefaultAssembler2D::assemblePDESystem(AbstractSystemMatrix* mat,
175 Data& rhs, const DataMap& coefs) const
176 {
177 const Data& A = unpackData("A", coefs);
178 const Data& B = unpackData("B", coefs);
179 const Data& C = unpackData("C", coefs);
180 const Data& D = unpackData("D", coefs);
181 const Data& X = unpackData("X", coefs);
182 const Data& Y = unpackData("Y", coefs);
183 assemblePDESystem(mat, rhs, A, B, C, D, X, Y);
184 }
185
186 void DefaultAssembler2D::assemblePDEBoundarySystem(AbstractSystemMatrix* mat,
187 Data& rhs, const DataMap& coefs) const
188 {
189 const Data& d = unpackData("d", coefs);
190 const Data& y = unpackData("y", coefs);
191 assemblePDEBoundarySystem(mat, rhs, d, y);
192 }
193
194 void DefaultAssembler2D::assemblePDESystemReduced(AbstractSystemMatrix* mat,
195 Data& rhs, const DataMap& coefs) const
196 {
197 const Data& A = unpackData("A", coefs);
198 const Data& B = unpackData("B", coefs);
199 const Data& C = unpackData("C", coefs);
200 const Data& D = unpackData("D", coefs);
201 const Data& X = unpackData("X", coefs);
202 const Data& Y = unpackData("Y", coefs);
203 assemblePDESystemReduced(mat, rhs, A, B, C, D, X, Y);
204 }
205
206 void DefaultAssembler2D::assemblePDEBoundarySystemReduced(
207 AbstractSystemMatrix* mat, Data& rhs,
208 const DataMap& coefs) const
209 {
210 const Data& d = unpackData("d", coefs);
211 const Data& y = unpackData("y", coefs);
212 assemblePDEBoundarySystemReduced(mat, rhs, d, y);
213 }
214
215 void DefaultAssembler2D::assemblePDESystem(AbstractSystemMatrix* mat,
216 Data& rhs, const Data& A, const Data& B,
217 const Data& C, const Data& D,
218 const Data& X, const Data& Y) const
219 {
220 if (!(A.isEmpty() && B.isEmpty() && C.isEmpty()))
221 throw SpeckleyException("Speckley does not support PDEs using A, B or C");
222 int order = domain->m_order;
223 const double *weights = all_weights[order-2];
224 const double volume_product = m_dx[0]*m_dx[1]/4.;
225 const int NE0 = m_NE[0];
226 const int NE1 = m_NE[1];
227 const int quads = order + 1;
228 const int max_x = m_NN[0];
229 dim_t numComp;
230 if (!mat)
231 numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
232 else {
233 numComp=mat->getColumnBlockSize();
234 }
235
236 rhs.requireWrite();
237 int size = 1;
238 if (!Y.isEmpty()) {
239 size = Y.getDataPointSize();
240 }
241 const index_t y_indices[2] = {0,size-1};
242 size = 1;
243 if (!D.isEmpty()) {
244 size = D.getDataPointSize();
245 }
246 const index_t d_indices[2] = {0,size-1};
247
248 if (!D.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
249 throw SpeckleyException("Speckley does not support adding left and right sides concurrently");
250
251 for (dim_t colouring = 0; colouring < 2; colouring++) {
252 #pragma omp parallel for
253 for (dim_t ey = colouring; ey < NE1; ey += 2) {
254 for (dim_t ex = 0; ex < NE0; ex++) {
255 const index_t e_index = ex + ey*NE0; //element number for Elements
256 const index_t start = ex*order + ey*max_x*order; //node start for Nodes
257
258 if (!D.isEmpty()) {
259 const double *e = D.getSampleDataRO(e_index);
260 if (D.actsExpanded()) {
261 for (index_t comp = 0; comp < numComp; comp++) {
262 for (short qy = 0; qy < quads; qy++) {
263 for (short qx = 0; qx < quads; qx++) {
264 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
265 /* D */ out[comp] += volume_product * weights[qx] * weights[qy]
266 * e[INDEX3(comp, qx, qy, numComp, quads)]; //D(i),qx,qy
267 }
268 }
269 }
270 } else { //constant
271 for (index_t comp = 0; comp < numComp; comp++) {
272 for (short qy = 0; qy < quads; qy++) {
273 for (short qx = 0; qx < quads; qx++) {
274 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
275 /* D */ out[comp] += volume_product * weights[qx] * weights[qy]
276 * e[d_indices[comp]]; //D(i)
277 }
278 }
279 }
280 }
281 }
282
283 if (!X.isEmpty()) {
284 const double *e = X.getSampleDataRO(e_index);
285 if (X.actsExpanded()) {
286 for (index_t comp = 0; comp < numComp; comp++) {
287 for (short qy = 0; qy < quads; qy++) {
288 for (short qx = 0; qx < quads; qx++) {
289 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
290 double res_a = 0, res_b = 0;
291 for (short k = 0; k < quads; k++) {
292 res_a += weights[k] * all_lagrange_derivs[order-2][qx][k]
293 * e[INDEX4(comp,0,k,qy,numComp,2,quads)]; //X(i1),k,qy
294 res_b += weights[k] * all_lagrange_derivs[order-2][qy][k]
295 * e[INDEX4(comp,1,qx,k,numComp,2,quads)]; //X(i2),qx,k
296 }
297 /* X */ out[comp] += 2 * volume_product
298 * (res_a * weights[qy] / m_dx[0] + res_b * weights[qx] / m_dx[1]);
299 }
300 }
301 }
302 } else { //constant
303 for (index_t comp = 0; comp < numComp; comp++) {
304 for (short qy = 0; qy < quads; qy++) {
305 for (short qx = 0; qx < quads; qx++) {
306 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
307 double res_a = 0, res_b = 0;
308 for (short k = 0; k < quads; k++) {
309 res_a += weights[k] * all_lagrange_derivs[order-2][qx][k]
310 * e[INDEX2(comp,0,numComp)]; //X(i1)
311 res_b += weights[k] * all_lagrange_derivs[order-2][qy][k]
312 * e[INDEX2(comp,1,numComp)]; //X(i2)
313 }
314 /* X */ out[comp] += 2 * volume_product
315 * (res_a * weights[qy] / m_dx[0] + res_b * weights[qx] / m_dx[1]);
316 }
317 }
318 }
319 }
320 }
321
322 if (!Y.isEmpty()) {
323 const double *e = Y.getSampleDataRO(e_index);
324 if (Y.actsExpanded()) {
325 for (index_t comp = 0; comp < numComp; comp++) {
326 for (short qy = 0; qy < quads; qy++) {
327 for (short qx = 0; qx < quads; qx++) {
328 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
329 /* Y */ out[comp] += volume_product * weights[qx] * weights[qy]
330 * e[INDEX3(comp, qx, qy, numComp, quads)]; //Y(i),qx,qy
331 }
332 }
333 }
334 } else { //constant
335 for (index_t comp = 0; comp < numComp; comp++) {
336 for (short qy = 0; qy < quads; qy++) {
337 for (short qx = 0; qx < quads; qx++) {
338 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
339 /* Y */ out[comp] += volume_product * weights[qx] * weights[qy]
340 * e[y_indices[comp]]; //Y(i)
341 }
342 }
343 }
344 }
345 }
346 }
347 }
348 }
349 }
350
351 void DefaultAssembler2D::assemblePDESystemReduced(AbstractSystemMatrix* mat,
352 Data& rhs, const Data& A, const Data& B,
353 const Data& C, const Data& D,
354 const Data& X, const Data& Y) const
355 {
356 throw SpeckleyException("Speckley does not support reduced functionspaces");
357
358 }
359
360 void DefaultAssembler2D::assemblePDEBoundarySingle(AbstractSystemMatrix* mat,
361 Data& rhs, const Data& d, const Data& y) const
362 {
363 throw SpeckleyException("Speckley does not support boundary functionspaces");
364 }
365
366 void DefaultAssembler2D::assemblePDEBoundarySingleReduced(
367 AbstractSystemMatrix* mat, Data& rhs,
368 const Data& d, const Data& y) const
369 {
370 throw SpeckleyException("Speckley does not support reduced functionspaces");
371 }
372
373 void DefaultAssembler2D::assemblePDEBoundarySystem(AbstractSystemMatrix* mat,
374 Data& rhs, const Data& d, const Data& y) const
375 {
376 throw SpeckleyException("Speckley does not support boundary functionspaces");
377 }
378
379 void DefaultAssembler2D::assemblePDEBoundarySystemReduced(
380 AbstractSystemMatrix* mat, Data& rhs,
381 const Data& d, const Data& y) const
382 {
383 throw SpeckleyException("Speckley does not support reduced functionspaces");
384 }
385
386 //protected
387 void DefaultAssembler2D::assemblePDESingle(AbstractSystemMatrix *mat,
388 escript::Data& rhs, const escript::Data& A, const escript::Data& B,
389 const escript::Data& C, const escript::Data& D,
390 const escript::Data& X, const escript::Data& Y) const
391 {
392 if (!(A.isEmpty() && B.isEmpty() && C.isEmpty()))
393 throw SpeckleyException("Speckley does not support PDEs using A, B or C");
394 int order = domain->m_order;
395 const double *weights = all_weights[order-2];
396 const double volume_product = m_dx[0]*m_dx[1]/4.;
397 const int NE0 = m_NE[0];
398 const int NE1 = m_NE[1];
399 const int quads = order + 1;
400 const int max_x = m_NN[0];
401 rhs.requireWrite();
402
403
404 if (!D.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
405 throw SpeckleyException("Speckley does not support adding left and right sides concurrently");
406
407 for (dim_t colouring = 0; colouring < 2; colouring++) {
408 #pragma omp parallel for
409 for (dim_t ey = colouring; ey < NE1; ey += 2) {
410 for (dim_t ex = 0; ex < NE0; ex++) {
411 const index_t e_index = INDEX2(ex,ey,NE0); //element number for Elements
412 const index_t start = order * (INDEX2(ex, ey, max_x)); //node start for Nodes
413
414 if (!D.isEmpty()) {
415 const double *e = D.getSampleDataRO(e_index);
416 if (D.actsExpanded()) {
417 for (short qy = 0; qy < quads; qy++) {
418 for (short qx = 0; qx < quads; qx++) {
419 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
420 /* D */ out[0] += volume_product * weights[qx] * weights[qy]
421 * e[INDEX2(qx, qy, quads)]; //D,qx,qy
422 }
423 }
424 } else { // constant form
425 for (short qy = 0; qy < quads; qy++) {
426 for (short qx = 0; qx < quads; qx++) {
427 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
428 /* D */ out[0] += volume_product * weights[qx] * weights[qy] * e[0]; //D
429 }
430 }
431 }
432 }
433
434 if (!X.isEmpty()) {
435 const double *e = X.getSampleDataRO(e_index);
436 if (X.actsExpanded()) {
437 for (short qy = 0; qy < quads; qy++) {
438 for (short qx = 0; qx < quads; qx++) {
439 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
440 double res_a = 0, res_b = 0;
441 for (short k = 0; k < quads; k++) {
442 res_a += weights[k] * all_lagrange_derivs[order-2][qx][k]
443 * e[INDEX3(0,k,qy,2,quads)]; //X(1),k,qy
444 res_b += weights[k] * all_lagrange_derivs[order-2][qy][k]
445 * e[INDEX3(1,qx,k,2,quads)]; //X(2),qx,k
446 }
447 /* X */ out[0] += 2 * volume_product * (res_a * weights[qy] / m_dx[0] + res_b * weights[qx] / m_dx[1]);
448 }
449 }
450 } else { //constant
451 for (short qy = 0; qy < quads; qy++) {
452 for (short qx = 0; qx < quads; qx++) {
453 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
454 double res_a = 0, res_b = 0;
455 for (short k = 0; k < quads; k++) {
456 res_a += weights[k] * all_lagrange_derivs[order-2][qx][k]
457 * e[0]; //X(1)
458 res_b += weights[k] * all_lagrange_derivs[order-2][qy][k]
459 * e[1]; //X(2)
460 }
461 /* X */ out[0] += 2 * volume_product
462 * (res_a * weights[qy] / m_dx[0] + res_b * weights[qx] / m_dx[1]);
463 }
464 }
465 }
466 }
467
468 if (!Y.isEmpty()) {
469 const double *e = Y.getSampleDataRO(e_index);
470 if (Y.actsExpanded()) {
471 for (short qy = 0; qy < quads; qy++) {
472 for (short qx = 0; qx < quads; qx++) {
473 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
474 /* Y */ out[0] += volume_product * weights[qx] * weights[qy]
475 * e[INDEX2(qx, qy, quads)]; //Y,qx,qy
476 }
477 }
478 } else { // constant form
479 for (short qy = 0; qy < quads; qy++) {
480 for (short qx = 0; qx < quads; qx++) {
481 double *out = rhs.getSampleDataRW(start + INDEX2(qx,qy,max_x));
482 /* Y */ out[0] += volume_product * weights[qx] * weights[qy] * e[0]; //Y
483 }
484 }
485 }
486 }
487 }
488 }
489 }
490 }
491
492 //protected
493 void DefaultAssembler2D::assemblePDESingleReduced(AbstractSystemMatrix* mat,
494 Data& rhs, const Data& A, const Data& B,
495 const Data& C, const Data& D,
496 const Data& X, const Data& Y) const
497 {
498 throw SpeckleyException("Speckley does not support reduced functionspaces");
499 }
500
501
502 } // namespace speckley
503

  ViewVC Help
Powered by ViewVC 1.1.26