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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4692 - (show annotations)
Wed Feb 19 05:26:55 2014 UTC (5 years, 6 months ago) by sshaw
File size: 46126 byte(s)
fixed duplication for 2D constant data in wave assembler
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2014 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 #include <ripley/WaveAssembler2D.h>
17
18 using namespace std;
19
20 namespace ripley {
21
22 WaveAssembler2D::WaveAssembler2D(Rectangle *dom, double *m_dx, dim_t *m_NX, dim_t *m_NE,
23 dim_t *m_NN, std::map<std::string, escript::Data> c)
24 : AbstractAssembler() {
25 domain = dom;
26 this->m_dx = m_dx;
27 this->m_NX = m_NX;
28 this->m_NE = m_NE;
29 this->m_NN = m_NN;
30 isHTI = isVTI = false;
31 std::map<std::string, escript::Data>::iterator a = c.find("c12"),
32 b = c.find("c23");
33 if (c.find("c11") == c.end()
34 || c.find("c13") == c.end() || c.find("c33") == c.end()
35 || c.find("c44") == c.end() || c.find("c66") == c.end()
36 || (a == c.end() && b == c.end()))
37 throw RipleyException("required constants missing for WaveAssembler");
38
39 if (a != c.end() && b != c.end()) {
40 throw RipleyException("WaveAssembler2D() doesn't support general"
41 " form waves (yet)");
42 } else if (a == c.end()) {
43 c23 = b->second;
44 isHTI = true;
45 } else if (b == c.end()) {
46 c12 = a->second;
47 isVTI = true;
48 } //final else case taken care of with the missing constants above
49 c11 = c.find("c11")->second,
50 c13 = c.find("c13")->second, c33 = c.find("c33")->second,
51 c44 = c.find("c44")->second, c66 = c.find("c66")->second;
52 }
53
54 void WaveAssembler2D::assemblePDESystem(Paso_SystemMatrix* mat,
55 escript::Data& rhs, map<string, escript::Data> coefs) const
56 {
57 const escript::Data A = unpackData("A", coefs), B = unpackData("B", coefs),
58 C = unpackData("C", coefs), D = unpackData("D", coefs),
59 Y = unpackData("Y", coefs), du = unpackData("du", coefs);
60 dim_t numEq, numComp;
61 if (!mat)
62 numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
63 else {
64 numEq=mat->logical_row_block_size;
65 numComp=mat->logical_col_block_size;
66 }
67 const double SQRT3 = 1.73205080756887719318;
68 const double w1 = 1.0/24;
69 const double w5 = -SQRT3/24 + 1.0/12;
70 const double w2 = -SQRT3/24 - 1.0/12;
71 const double w19 = -m_dx[0]/12;
72 const double w11 = w19*(SQRT3 + 3)/12;
73 const double w14 = w19*(-SQRT3 + 3)/12;
74 const double w16 = w19*(5*SQRT3 + 9)/12;
75 const double w17 = w19*(-5*SQRT3 + 9)/12;
76 const double w27 = w19*(-SQRT3 - 3)/2;
77 const double w28 = w19*(SQRT3 - 3)/2;
78 const double w18 = -m_dx[1]/12;
79 const double w10 = w18*(SQRT3 + 3)/12;
80 const double w15 = w18*(-SQRT3 + 3)/12;
81 const double w12 = w18*(5*SQRT3 + 9)/12;
82 const double w13 = w18*(-5*SQRT3 + 9)/12;
83 const double w25 = w18*(-SQRT3 - 3)/2;
84 const double w26 = w18*(SQRT3 - 3)/2;
85 const double w22 = m_dx[0]*m_dx[1]/144;
86 const double w20 = w22*(SQRT3 + 2);
87 const double w21 = w22*(-SQRT3 + 2);
88 const double w23 = w22*(4*SQRT3 + 7);
89 const double w24 = w22*(-4*SQRT3 + 7);
90 const double w3 = m_dx[0]/(24*m_dx[1]);
91 const double w7 = w3*(SQRT3 + 2);
92 const double w8 = w3*(-SQRT3 + 2);
93 const double w6 = -m_dx[1]/(24*m_dx[0]);
94 const double w0 = w6*(SQRT3 + 2);
95 const double w4 = w6*(-SQRT3 + 2);
96
97 rhs.requireWrite();
98 #pragma omp parallel
99 {
100 for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
101 #pragma omp for
102 for (index_t k1=k1_0; k1 < m_NE[1]; k1+=2) {
103 for (index_t k0=0; k0 < m_NE[0]; ++k0) {
104 bool addEM_S=false;
105 bool addEM_F=false;
106 vector<double> EM_S(4*4*numEq*numComp, 0);
107 vector<double> EM_F(4*numEq, 0);
108 const index_t e = k0 + m_NE[0]*k1;
109 ///////////////
110 // process A //
111 ///////////////
112 if (!A.isEmpty()) {
113 addEM_S = true;
114 const double* A_p = A.getSampleDataRO(e);
115 if (A.actsExpanded()) {
116 for (index_t k=0; k<numEq; k++) {
117 for (index_t m=0; m<numComp; m++) {
118 const double A_00_0 = A_p[INDEX5(k,0,m,0,0,numEq,2,numComp,2)];
119 const double A_01_0 = A_p[INDEX5(k,0,m,1,0,numEq,2,numComp,2)];
120 const double A_10_0 = A_p[INDEX5(k,1,m,0,0,numEq,2,numComp,2)];
121 const double A_11_0 = A_p[INDEX5(k,1,m,1,0,numEq,2,numComp,2)];
122 const double A_00_1 = A_p[INDEX5(k,0,m,0,1,numEq,2,numComp,2)];
123 const double A_01_1 = A_p[INDEX5(k,0,m,1,1,numEq,2,numComp,2)];
124 const double A_10_1 = A_p[INDEX5(k,1,m,0,1,numEq,2,numComp,2)];
125 const double A_11_1 = A_p[INDEX5(k,1,m,1,1,numEq,2,numComp,2)];
126 const double A_00_2 = A_p[INDEX5(k,0,m,0,2,numEq,2,numComp,2)];
127 const double A_01_2 = A_p[INDEX5(k,0,m,1,2,numEq,2,numComp,2)];
128 const double A_10_2 = A_p[INDEX5(k,1,m,0,2,numEq,2,numComp,2)];
129 const double A_11_2 = A_p[INDEX5(k,1,m,1,2,numEq,2,numComp,2)];
130 const double A_00_3 = A_p[INDEX5(k,0,m,0,3,numEq,2,numComp,2)];
131 const double A_01_3 = A_p[INDEX5(k,0,m,1,3,numEq,2,numComp,2)];
132 const double A_10_3 = A_p[INDEX5(k,1,m,0,3,numEq,2,numComp,2)];
133 const double A_11_3 = A_p[INDEX5(k,1,m,1,3,numEq,2,numComp,2)];
134 const double tmp0 = w3*(A_11_0 + A_11_1 + A_11_2 + A_11_3);
135 const double tmp1 = w1*(A_01_0 + A_01_3 - A_10_1 - A_10_2);
136 const double tmp2 = w4*(A_00_2 + A_00_3);
137 const double tmp3 = w0*(A_00_0 + A_00_1);
138 const double tmp4 = w5*(A_01_2 - A_10_3);
139 const double tmp5 = w2*(-A_01_1 + A_10_0);
140 const double tmp6 = w5*(A_01_3 + A_10_0);
141 const double tmp7 = w3*(-A_11_0 - A_11_1 - A_11_2 - A_11_3);
142 const double tmp8 = w6*(A_00_0 + A_00_1 + A_00_2 + A_00_3);
143 const double tmp9 = w1*(A_01_1 + A_01_2 + A_10_1 + A_10_2);
144 const double tmp10 = w2*(-A_01_0 - A_10_3);
145 const double tmp11 = w4*(A_00_0 + A_00_1);
146 const double tmp12 = w0*(A_00_2 + A_00_3);
147 const double tmp13 = w5*(A_01_1 - A_10_0);
148 const double tmp14 = w2*(-A_01_2 + A_10_3);
149 const double tmp15 = w7*(A_11_0 + A_11_2);
150 const double tmp16 = w4*(-A_00_2 - A_00_3);
151 const double tmp17 = w0*(-A_00_0 - A_00_1);
152 const double tmp18 = w5*(A_01_3 + A_10_3);
153 const double tmp19 = w8*(A_11_1 + A_11_3);
154 const double tmp20 = w2*(-A_01_0 - A_10_0);
155 const double tmp21 = w7*(A_11_1 + A_11_3);
156 const double tmp22 = w4*(-A_00_0 - A_00_1);
157 const double tmp23 = w0*(-A_00_2 - A_00_3);
158 const double tmp24 = w5*(A_01_0 + A_10_0);
159 const double tmp25 = w8*(A_11_0 + A_11_2);
160 const double tmp26 = w2*(-A_01_3 - A_10_3);
161 const double tmp27 = w5*(-A_01_1 - A_10_2);
162 const double tmp28 = w1*(-A_01_0 - A_01_3 - A_10_0 - A_10_3);
163 const double tmp29 = w2*(A_01_2 + A_10_1);
164 const double tmp30 = w7*(-A_11_1 - A_11_3);
165 const double tmp31 = w1*(-A_01_1 - A_01_2 + A_10_0 + A_10_3);
166 const double tmp32 = w5*(-A_01_0 + A_10_2);
167 const double tmp33 = w8*(-A_11_0 - A_11_2);
168 const double tmp34 = w6*(-A_00_0 - A_00_1 - A_00_2 - A_00_3);
169 const double tmp35 = w2*(A_01_3 - A_10_1);
170 const double tmp36 = w5*(A_01_0 + A_10_3);
171 const double tmp37 = w2*(-A_01_3 - A_10_0);
172 const double tmp38 = w7*(-A_11_0 - A_11_2);
173 const double tmp39 = w5*(-A_01_3 + A_10_1);
174 const double tmp40 = w8*(-A_11_1 - A_11_3);
175 const double tmp41 = w2*(A_01_0 - A_10_2);
176 const double tmp42 = w5*(A_01_1 - A_10_3);
177 const double tmp43 = w2*(-A_01_2 + A_10_0);
178 const double tmp44 = w5*(A_01_2 - A_10_0);
179 const double tmp45 = w2*(-A_01_1 + A_10_3);
180 const double tmp46 = w5*(-A_01_0 + A_10_1);
181 const double tmp47 = w2*(A_01_3 - A_10_2);
182 const double tmp48 = w5*(-A_01_1 - A_10_1);
183 const double tmp49 = w2*(A_01_2 + A_10_2);
184 const double tmp50 = w5*(-A_01_3 + A_10_2);
185 const double tmp51 = w2*(A_01_0 - A_10_1);
186 const double tmp52 = w5*(-A_01_2 - A_10_1);
187 const double tmp53 = w2*(A_01_1 + A_10_2);
188 const double tmp54 = w5*(-A_01_2 - A_10_2);
189 const double tmp55 = w2*(A_01_1 + A_10_1);
190 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp15 + tmp16 + tmp17 + tmp18 + tmp19 + tmp20 + tmp9;
191 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5;
192 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp31 + tmp34 + tmp38 + tmp39 + tmp40 + tmp41;
193 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp28 + tmp52 + tmp53 + tmp7 + tmp8;
194 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0 + tmp2 + tmp3 + tmp31 + tmp50 + tmp51;
195 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp16 + tmp17 + tmp21 + tmp25 + tmp28 + tmp54 + tmp55;
196 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10 + tmp6 + tmp7 + tmp8 + tmp9;
197 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1 + tmp30 + tmp33 + tmp34 + tmp44 + tmp45;
198 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1 + tmp34 + tmp38 + tmp40 + tmp42 + tmp43;
199 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp36 + tmp37 + tmp7 + tmp8 + tmp9;
200 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp15 + tmp19 + tmp22 + tmp23 + tmp28 + tmp48 + tmp49;
201 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0 + tmp11 + tmp12 + tmp31 + tmp46 + tmp47;
202 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp27 + tmp28 + tmp29 + tmp7 + tmp8;
203 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp30 + tmp31 + tmp32 + tmp33 + tmp34 + tmp35;
204 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0 + tmp1 + tmp11 + tmp12 + tmp13 + tmp14;
205 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp21 + tmp22 + tmp23 + tmp24 + tmp25 + tmp26 + tmp9;
206 }
207 }
208 } else { // constant data
209 for (index_t k=0; k<numEq; k++) {
210 for (index_t m=0; m<numComp; m++) {
211 const double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
212 const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
213 const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
214 const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
215 const double tmp0 = 6*w1*(A_01 - A_10);
216 const double tmp1 = 6*w1*(A_01 + A_10);
217 const double tmp2 = 6*w1*(-A_01 - A_10);
218 const double tmp3 = 6*w1*(-A_01 + A_10);
219 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp1;
220 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp0;
221 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp3;
222 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp2;
223 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp3;
224 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp2;
225 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp1;
226 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp0;
227 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp0;
228 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp1;
229 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp2;
230 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp3;
231 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp2;
232 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp3;
233 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp0;
234 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp1;
235 }
236 }
237 }
238 }
239 ///////////////
240 // process B //
241 ///////////////
242 if (!B.isEmpty()) {
243 addEM_S=true;
244 const double* B_p=B.getSampleDataRO(e);
245 if (B.actsExpanded()) {
246 for (index_t k=0; k<numEq; k++) {
247 for (index_t m=0; m<numComp; m++) {
248 const double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];
249 const double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];
250 const double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];
251 const double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];
252 const double B_0_2 = B_p[INDEX4(k,0,m,2, numEq,2,numComp)];
253 const double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];
254 const double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];
255 const double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];
256 const double tmp0 = w11*(B_1_0 + B_1_1);
257 const double tmp1 = w14*(B_1_2 + B_1_3);
258 const double tmp2 = w15*(-B_0_1 - B_0_3);
259 const double tmp3 = w10*(-B_0_0 - B_0_2);
260 const double tmp4 = w11*(B_1_2 + B_1_3);
261 const double tmp5 = w14*(B_1_0 + B_1_1);
262 const double tmp6 = w11*(-B_1_2 - B_1_3);
263 const double tmp7 = w14*(-B_1_0 - B_1_1);
264 const double tmp8 = w11*(-B_1_0 - B_1_1);
265 const double tmp9 = w14*(-B_1_2 - B_1_3);
266 const double tmp10 = w10*(-B_0_1 - B_0_3);
267 const double tmp11 = w15*(-B_0_0 - B_0_2);
268 const double tmp12 = w15*(B_0_0 + B_0_2);
269 const double tmp13 = w10*(B_0_1 + B_0_3);
270 const double tmp14 = w10*(B_0_0 + B_0_2);
271 const double tmp15 = w15*(B_0_1 + B_0_3);
272 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=B_0_0*w12 + B_0_1*w10 + B_0_2*w15 + B_0_3*w13 + B_1_0*w16 + B_1_1*w14 + B_1_2*w11 + B_1_3*w17;
273 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=B_0_0*w10 + B_0_1*w12 + B_0_2*w13 + B_0_3*w15 + tmp0 + tmp1;
274 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=B_1_0*w11 + B_1_1*w17 + B_1_2*w16 + B_1_3*w14 + tmp14 + tmp15;
275 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp12 + tmp13 + tmp4 + tmp5;
276 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-B_0_0*w12 - B_0_1*w10 - B_0_2*w15 - B_0_3*w13 + tmp0 + tmp1;
277 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-B_0_0*w10 - B_0_1*w12 - B_0_2*w13 - B_0_3*w15 + B_1_0*w14 + B_1_1*w16 + B_1_2*w17 + B_1_3*w11;
278 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
279 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=B_1_0*w17 + B_1_1*w11 + B_1_2*w14 + B_1_3*w16 + tmp10 + tmp11;
280 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-B_1_0*w16 - B_1_1*w14 - B_1_2*w11 - B_1_3*w17 + tmp14 + tmp15;
281 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
282 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=B_0_0*w15 + B_0_1*w13 + B_0_2*w12 + B_0_3*w10 - B_1_0*w11 - B_1_1*w17 - B_1_2*w16 - B_1_3*w14;
283 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=B_0_0*w13 + B_0_1*w15 + B_0_2*w10 + B_0_3*w12 + tmp6 + tmp7;
284 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2 + tmp3 + tmp8 + tmp9;
285 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-B_1_0*w14 - B_1_1*w16 - B_1_2*w17 - B_1_3*w11 + tmp10 + tmp11;
286 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-B_0_0*w15 - B_0_1*w13 - B_0_2*w12 - B_0_3*w10 + tmp6 + tmp7;
287 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-B_0_0*w13 - B_0_1*w15 - B_0_2*w10 - B_0_3*w12 - B_1_0*w17 - B_1_1*w11 - B_1_2*w14 - B_1_3*w16;
288 }
289 }
290 } else { // constant data
291 for (index_t k=0; k<numEq; k++) {
292 for (index_t m=0; m<numComp; m++) {
293 const double wB0 = B_p[INDEX3(k,0,m,numEq,2)]*w18;
294 const double wB1 = B_p[INDEX3(k,1,m,numEq,2)]*w19;
295 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+= 2*wB0 + 2*wB1;
296 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+= 2*wB0 + wB1;
297 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+= wB0 + 2*wB1;
298 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+= wB0 + wB1;
299 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-2*wB0 + wB1;
300 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-2*wB0 + 2*wB1;
301 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+= -wB0 + wB1;
302 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+= -wB0 + 2*wB1;
303 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+= wB0 - 2*wB1;
304 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+= wB0 - wB1;
305 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+= 2*wB0 - 2*wB1;
306 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+= 2*wB0 - wB1;
307 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+= -wB0 - wB1;
308 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+= -wB0 - 2*wB1;
309 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-2*wB0 - wB1;
310 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-2*wB0 - 2*wB1;
311 }
312 }
313 }
314 }
315 ///////////////
316 // process C //
317 ///////////////
318 if (!C.isEmpty()) {
319 addEM_S=true;
320 const double* C_p=C.getSampleDataRO(e);
321 if (C.actsExpanded()) {
322 for (index_t k=0; k<numEq; k++) {
323 for (index_t m=0; m<numComp; m++) {
324 const double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];
325 const double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];
326 const double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];
327 const double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];
328 const double C_0_2 = C_p[INDEX4(k,m,0, 2, numEq,numComp,2)];
329 const double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];
330 const double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];
331 const double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];
332 const double tmp0 = w11*(C_1_0 + C_1_1);
333 const double tmp1 = w14*(C_1_2 + C_1_3);
334 const double tmp2 = w15*(C_0_0 + C_0_2);
335 const double tmp3 = w10*(C_0_1 + C_0_3);
336 const double tmp4 = w11*(-C_1_0 - C_1_1);
337 const double tmp5 = w14*(-C_1_2 - C_1_3);
338 const double tmp6 = w11*(-C_1_2 - C_1_3);
339 const double tmp7 = w14*(-C_1_0 - C_1_1);
340 const double tmp8 = w11*(C_1_2 + C_1_3);
341 const double tmp9 = w14*(C_1_0 + C_1_1);
342 const double tmp10 = w10*(-C_0_1 - C_0_3);
343 const double tmp11 = w15*(-C_0_0 - C_0_2);
344 const double tmp12 = w15*(-C_0_1 - C_0_3);
345 const double tmp13 = w10*(-C_0_0 - C_0_2);
346 const double tmp14 = w10*(C_0_0 + C_0_2);
347 const double tmp15 = w15*(C_0_1 + C_0_3);
348 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=C_0_0*w12 + C_0_1*w10 + C_0_2*w15 + C_0_3*w13 + C_1_0*w16 + C_1_1*w14 + C_1_2*w11 + C_1_3*w17;
349 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-C_0_0*w12 - C_0_1*w10 - C_0_2*w15 - C_0_3*w13 + tmp0 + tmp1;
350 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-C_1_0*w16 - C_1_1*w14 - C_1_2*w11 - C_1_3*w17 + tmp14 + tmp15;
351 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp12 + tmp13 + tmp4 + tmp5;
352 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=C_0_0*w10 + C_0_1*w12 + C_0_2*w13 + C_0_3*w15 + tmp0 + tmp1;
353 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-C_0_0*w10 - C_0_1*w12 - C_0_2*w13 - C_0_3*w15 + C_1_0*w14 + C_1_1*w16 + C_1_2*w17 + C_1_3*w11;
354 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
355 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-C_1_0*w14 - C_1_1*w16 - C_1_2*w17 - C_1_3*w11 + tmp10 + tmp11;
356 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=C_1_0*w11 + C_1_1*w17 + C_1_2*w16 + C_1_3*w14 + tmp14 + tmp15;
357 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
358 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=C_0_0*w15 + C_0_1*w13 + C_0_2*w12 + C_0_3*w10 - C_1_0*w11 - C_1_1*w17 - C_1_2*w16 - C_1_3*w14;
359 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-C_0_0*w15 - C_0_1*w13 - C_0_2*w12 - C_0_3*w10 + tmp6 + tmp7;
360 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2 + tmp3 + tmp8 + tmp9;
361 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=C_1_0*w17 + C_1_1*w11 + C_1_2*w14 + C_1_3*w16 + tmp10 + tmp11;
362 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=C_0_0*w13 + C_0_1*w15 + C_0_2*w10 + C_0_3*w12 + tmp6 + tmp7;
363 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-C_0_0*w13 - C_0_1*w15 - C_0_2*w10 - C_0_3*w12 - C_1_0*w17 - C_1_1*w11 - C_1_2*w14 - C_1_3*w16;
364 }
365 }
366 } else { // constant data
367 for (index_t k=0; k<numEq; k++) {
368 for (index_t m=0; m<numComp; m++) {
369 const double wC0 = C_p[INDEX3(k,m,0,numEq,numComp)]*w18;
370 const double wC1 = C_p[INDEX3(k,m,1,numEq,numComp)]*w19;
371 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+= 2*wC0 + 2*wC1;
372 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-2*wC0 + wC1;
373 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+= wC0 - 2*wC1;
374 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+= -wC0 - wC1;
375 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+= 2*wC0 + wC1;
376 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-2*wC0 + 2*wC1;
377 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+= wC0 - wC1;
378 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+= -wC0 - 2*wC1;
379 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+= wC0 + 2*wC1;
380 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+= -wC0 + wC1;
381 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+= 2*wC0 - 2*wC1;
382 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-2*wC0 - wC1;
383 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+= wC0 + wC1;
384 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+= -wC0 + 2*wC1;
385 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+= 2*wC0 - wC1;
386 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-2*wC0 - 2*wC1;
387 }
388 }
389 }
390 }
391 ///////////////
392 // process D //
393 ///////////////
394 if (!D.isEmpty()) {
395 addEM_S=true;
396 const double* D_p=D.getSampleDataRO(e);
397 if (D.actsExpanded()) {
398 for (index_t k=0; k<numEq; k++) {
399 for (index_t m=0; m<numComp; m++) {
400 const double D_0 = D_p[INDEX3(k,m,0,numEq,numComp)];
401 const double D_1 = D_p[INDEX3(k,m,1,numEq,numComp)];
402 const double D_2 = D_p[INDEX3(k,m,2,numEq,numComp)];
403 const double D_3 = D_p[INDEX3(k,m,3,numEq,numComp)];
404 const double tmp0 = w21*(D_2 + D_3);
405 const double tmp1 = w20*(D_0 + D_1);
406 const double tmp2 = w22*(D_0 + D_1 + D_2 + D_3);
407 const double tmp3 = w21*(D_0 + D_1);
408 const double tmp4 = w20*(D_2 + D_3);
409 const double tmp5 = w22*(D_1 + D_2);
410 const double tmp6 = w21*(D_0 + D_2);
411 const double tmp7 = w20*(D_1 + D_3);
412 const double tmp8 = w21*(D_1 + D_3);
413 const double tmp9 = w20*(D_0 + D_2);
414 const double tmp10 = w22*(D_0 + D_3);
415 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=D_0*w23 + D_3*w24 + tmp5;
416 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0 + tmp1;
417 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp8 + tmp9;
418 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp2;
419 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0 + tmp1;
420 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=D_1*w23 + D_2*w24 + tmp10;
421 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2;
422 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp6 + tmp7;
423 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp8 + tmp9;
424 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp2;
425 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=D_1*w24 + D_2*w23 + tmp10;
426 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp3 + tmp4;
427 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2;
428 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp6 + tmp7;
429 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3 + tmp4;
430 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=D_0*w24 + D_3*w23 + tmp5;
431 }
432 }
433 } else { // constant data
434 for (index_t k=0; k<numEq; k++) {
435 for (index_t m=0; m<numComp; m++) {
436 const double D_0 = D_p[INDEX2(k, m, numEq)];
437 EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=16*D_0*w22;
438 EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=8*D_0*w22;
439 EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=8*D_0*w22;
440 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=4*D_0*w22;
441 EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=8*D_0*w22;
442 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=16*D_0*w22;
443 EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=4*D_0*w22;
444 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=8*D_0*w22;
445 EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=8*D_0*w22;
446 EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=4*D_0*w22;
447 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=16*D_0*w22;
448 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=8*D_0*w22;
449 EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=4*D_0*w22;
450 EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=8*D_0*w22;
451 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=8*D_0*w22;
452 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=16*D_0*w22;
453 }
454 }
455 }
456 }
457 ///////////////
458 // process X //
459 ///////////////
460 if (!du.isEmpty()) {
461 addEM_F=true;
462 const double *du_p = du.getSampleDataRO(e);
463 const double *c11_p = c11.getSampleDataRO(e),
464 *c13_p = c13.getSampleDataRO(e),
465 *c33_p = c33.getSampleDataRO(e);
466 if (du.actsExpanded()) {
467 double X_00_0, X_00_1, X_00_2, X_00_3,
468 X_11_0, X_11_1, X_11_2, X_11_3,
469 X_01_0 = -(du_p[INDEX3(1,0,0,numEq,2)] + du_p[INDEX3(0,1,0,numEq,2)]),
470 X_01_1 = -(du_p[INDEX3(1,0,1,numEq,2)] + du_p[INDEX3(0,1,1,numEq,2)]),
471 X_01_2 = -(du_p[INDEX3(1,0,2,numEq,2)] + du_p[INDEX3(0,1,2,numEq,2)]),
472 X_01_3 = -(du_p[INDEX3(1,0,3,numEq,2)] + du_p[INDEX3(0,1,3,numEq,2)]);
473 if (isVTI) {
474 const double *c44_p = c44.getSampleDataRO(e);
475 X_00_0 = -(du_p[INDEX3(0,0,0,numEq,2)] * c11_p[0]
476 + du_p[INDEX3(1,1,0,numEq,2)] * c13_p[0]);
477 X_00_1 = -(du_p[INDEX3(0,0,1,numEq,2)] * c11_p[1]
478 + du_p[INDEX3(1,1,1,numEq,2)] * c13_p[1]);
479 X_00_2 = -(du_p[INDEX3(0,0,2,numEq,2)] * c11_p[2]
480 + du_p[INDEX3(1,1,2,numEq,2)] * c13_p[2]);
481 X_00_3 = -(du_p[INDEX3(0,0,3,numEq,2)] * c11_p[3]
482 + du_p[INDEX3(1,1,3,numEq,2)] * c13_p[3]);
483 X_11_0 = -(du_p[INDEX3(0,0,0,numEq,2)] * c13_p[0]
484 + du_p[INDEX3(1,1,0,numEq,2)] * c33_p[0]);
485 X_11_1 = -(du_p[INDEX3(0,0,1,numEq,2)] * c13_p[1]
486 + du_p[INDEX3(1,1,1,numEq,2)] * c33_p[1]);
487 X_11_2 = -(du_p[INDEX3(0,0,2,numEq,2)] * c13_p[2]
488 + du_p[INDEX3(1,1,2,numEq,2)] * c33_p[2]);
489 X_11_3 = -(du_p[INDEX3(0,0,3,numEq,2)] * c13_p[3]
490 + du_p[INDEX3(1,1,3,numEq,2)] * c33_p[3]);
491 X_01_0 *= c44_p[0];
492 X_01_1 *= c44_p[1];
493 X_01_2 *= c44_p[2];
494 X_01_3 *= c44_p[3];
495 } else if (isHTI) {
496 const double *c66_p = c66.getSampleDataRO(e);
497 X_00_0 = -(du_p[INDEX3(0,0,0,numEq,2)] * c11_p[0]
498 + du_p[INDEX3(1,1,0,numEq,2)] * c13_p[0]);
499 X_00_1 = -(du_p[INDEX3(0,0,1,numEq,2)] * c11_p[1]
500 + du_p[INDEX3(1,1,1,numEq,2)] * c13_p[1]);
501 X_00_2 = -(du_p[INDEX3(0,0,2,numEq,2)] * c11_p[2]
502 + du_p[INDEX3(1,1,2,numEq,2)] * c13_p[2]);
503 X_00_3 = -(du_p[INDEX3(0,0,3,numEq,2)] * c11_p[3]
504 + du_p[INDEX3(1,1,3,numEq,2)] * c13_p[3]);
505 X_11_0 = -(du_p[INDEX3(0,0,0,numEq,2)] * c13_p[0]
506 + du_p[INDEX3(1,1,0,numEq,2)] * c33_p[0]);
507 X_11_1 = -(du_p[INDEX3(0,0,1,numEq,2)] * c13_p[1]
508 + du_p[INDEX3(1,1,1,numEq,2)] * c33_p[1]);
509 X_11_2 = -(du_p[INDEX3(0,0,2,numEq,2)] * c13_p[2]
510 + du_p[INDEX3(1,1,2,numEq,2)] * c33_p[2]);
511 X_11_3 = -(du_p[INDEX3(0,0,3,numEq,2)] * c13_p[3]
512 + du_p[INDEX3(1,1,3,numEq,2)] * c33_p[3]);
513 X_01_0 *= c66_p[0];
514 X_01_1 *= c66_p[1];
515 X_01_2 *= c66_p[2];
516 X_01_3 *= c66_p[3];
517 } else {
518 throw RipleyException("General form solutions"
519 " not yet implemented in WaveAssembler");
520 }
521
522 const double X_10_0 = X_01_0, X_10_1 = X_01_1, X_10_2 = X_01_2, X_10_3 = X_01_3;
523
524 const double Atmp0 = 6*w15*(X_00_2 + X_00_3);
525 const double Atmp1 = 6*w10*(X_00_0 + X_00_1);
526 const double Atmp2 = 6*w11*(X_01_0 + X_01_2);
527 const double Atmp3 = 6*w14*(X_01_1 + X_01_3);
528 const double Atmp4 = 6*w11*(X_01_1 + X_01_3);
529 const double Atmp5 = w25*(X_00_0 + X_00_1);
530 const double Atmp6 = w26*(X_00_2 + X_00_3);
531 const double Atmp7 = 6*w14*(X_01_0 + X_01_2);
532 const double Atmp8 = w27*(X_01_0 + X_01_2);
533 const double Atmp9 = w28*(X_01_1 + X_01_3);
534 const double Atmp10 = w25*(-X_00_2 - X_00_3);
535 const double Atmp11 = w26*(-X_00_0 - X_00_1);
536 const double Atmp12 = w27*(X_01_1 + X_01_3);
537 const double Atmp13 = w28*(X_01_0 + X_01_2);
538 const double Atmp14 = w25*(X_00_2 + X_00_3);
539 const double Atmp15 = w26*(X_00_0 + X_00_1);
540 EM_F[INDEX2(0,0,numEq)]+=Atmp0 + Atmp1 + Atmp2 + Atmp3;
541 EM_F[INDEX2(0,1,numEq)]+=Atmp4 + Atmp5 + Atmp6 + Atmp7;
542 EM_F[INDEX2(0,2,numEq)]+=Atmp10 + Atmp11 + Atmp8 + Atmp9;
543 EM_F[INDEX2(0,3,numEq)]+=Atmp12 + Atmp13 + Atmp14 + Atmp15;
544 const double Btmp0 = 6*w15*(X_10_2 + X_10_3);
545 const double Btmp1 = 6*w10*(X_10_0 + X_10_1);
546 const double Btmp2 = 6*w11*(X_11_0 + X_11_2);
547 const double Btmp3 = 6*w14*(X_11_1 + X_11_3);
548 const double Btmp4 = 6*w11*(X_11_1 + X_11_3);
549 const double Btmp5 = w25*(X_10_0 + X_10_1);
550 const double Btmp6 = w26*(X_10_2 + X_10_3);
551 const double Btmp7 = 6*w14*(X_11_0 + X_11_2);
552 const double Btmp8 = w27*(X_11_0 + X_11_2);
553 const double Btmp9 = w28*(X_11_1 + X_11_3);
554 const double Btmp10 = w25*(-X_10_2 - X_10_3);
555 const double Btmp11 = w26*(-X_10_0 - X_10_1);
556 const double Btmp12 = w27*(X_11_1 + X_11_3);
557 const double Btmp13 = w28*(X_11_0 + X_11_2);
558 const double Btmp14 = w25*(X_10_2 + X_10_3);
559 const double Btmp15 = w26*(X_10_0 + X_10_1);
560 EM_F[INDEX2(1,0,numEq)]+=Btmp0 + Btmp1 + Btmp2 + Btmp3;
561 EM_F[INDEX2(1,1,numEq)]+=Btmp4 + Btmp5 + Btmp6 + Btmp7;
562 EM_F[INDEX2(1,2,numEq)]+=Btmp10 + Btmp11 + Btmp8 + Btmp9;
563 EM_F[INDEX2(1,3,numEq)]+=Btmp12 + Btmp13 + Btmp14 + Btmp15;
564
565 } else { // constant data
566 double wX_00, wX_01, wX_10, wX_11;
567 if (isVTI) {
568 const double *c44_p = c44.getSampleDataRO(e);
569 wX_00 = -(du_p[INDEX2(0,0,numEq)] * c11_p[0]
570 + du_p[INDEX2(1,1,numEq)] * c13_p[0])*w18;
571 wX_01 = -(c44_p[0] *
572 (du_p[INDEX2(1,0,numEq)] + du_p[INDEX2(0,1,numEq)]))*w19;
573 wX_10 = -(c44_p[0] *
574 (du_p[INDEX2(1,0,numEq)] + du_p[INDEX2(0,1,numEq)]))*w18;
575 wX_11 = -(du_p[INDEX2(0,0,numEq)] * c13_p[0]
576 + du_p[INDEX2(1,1,numEq)] * c33_p[0])*w19;
577 } else if (isHTI) {
578 const double *c66_p = c66.getSampleDataRO(e);
579 wX_00 = -(du_p[INDEX2(0,0,numEq)] * c11_p[0]
580 + du_p[INDEX2(1,1,numEq)] * c13_p[0])*w18;
581 wX_01 = -(c66_p[0] *
582 (du_p[INDEX2(1,0,numEq)] + du_p[INDEX2(0,1,numEq)]))*w19;
583 wX_10 = -(c66_p[0] *
584 (du_p[INDEX2(1,0,numEq)] + du_p[INDEX2(0,1,numEq)]))*w18;
585 wX_11 = -(du_p[INDEX2(0,0,numEq)] * c13_p[0]
586 + du_p[INDEX2(1,1,numEq)] * c33_p[0])*w19;
587 } else {
588 throw RipleyException("General form solutions"
589 " not yet implemented in WaveAssembler");
590 }
591 EM_F[INDEX2(0,0,numEq)]+= wX_00 + wX_01;
592 EM_F[INDEX2(0,1,numEq)]+=-wX_00 + wX_01;
593 EM_F[INDEX2(0,2,numEq)]+= wX_00 - wX_01;
594 EM_F[INDEX2(0,3,numEq)]+=-wX_00 - wX_01;
595 EM_F[INDEX2(1,0,numEq)]+= wX_10 + wX_11;
596 EM_F[INDEX2(1,1,numEq)]+=-wX_10 + wX_11;
597 EM_F[INDEX2(1,2,numEq)]+= wX_10 - wX_11;
598 EM_F[INDEX2(1,3,numEq)]+=-wX_10 - wX_11;
599 }
600 }
601 ///////////////
602 // process Y //
603 ///////////////
604 if (!Y.isEmpty()) {
605 addEM_F=true;
606 const double* Y_p=Y.getSampleDataRO(e);
607 if (Y.actsExpanded()) {
608 for (index_t k=0; k<numEq; k++) {
609 const double Y_0 = Y_p[INDEX2(k, 0, numEq)];
610 const double Y_1 = Y_p[INDEX2(k, 1, numEq)];
611 const double Y_2 = Y_p[INDEX2(k, 2, numEq)];
612 const double Y_3 = Y_p[INDEX2(k, 3, numEq)];
613 const double tmp0 = 6*w22*(Y_1 + Y_2);
614 const double tmp1 = 6*w22*(Y_0 + Y_3);
615 EM_F[INDEX2(k,0,numEq)]+=6*Y_0*w20 + 6*Y_3*w21 + tmp0;
616 EM_F[INDEX2(k,1,numEq)]+=6*Y_1*w20 + 6*Y_2*w21 + tmp1;
617 EM_F[INDEX2(k,2,numEq)]+=6*Y_1*w21 + 6*Y_2*w20 + tmp1;
618 EM_F[INDEX2(k,3,numEq)]+=6*Y_0*w21 + 6*Y_3*w20 + tmp0;
619 }
620 } else { // constant data
621 for (index_t k=0; k<numEq; k++) {
622 EM_F[INDEX2(k,0,numEq)]+=36*Y_p[k]*w22;
623 EM_F[INDEX2(k,1,numEq)]+=36*Y_p[k]*w22;
624 EM_F[INDEX2(k,2,numEq)]+=36*Y_p[k]*w22;
625 EM_F[INDEX2(k,3,numEq)]+=36*Y_p[k]*w22;
626 }
627 }
628 }
629
630 // add to matrix (if addEM_S) and RHS (if addEM_F)
631 const index_t firstNode=m_NN[0]*k1+k0;
632 domain->addToMatrixAndRHS(mat, rhs, EM_S, EM_F, addEM_S,
633 addEM_F, firstNode, numEq, numComp);
634 } // end k0 loop
635 } // end k1 loop
636 } // end of colouring
637 } // end of parallel region
638 }
639
640 }
641

  ViewVC Help
Powered by ViewVC 1.1.26