/[escript]/branches/diaplayground/ripley/src/Brick.cpp
ViewVC logotype

Diff of /branches/diaplayground/ripley/src/Brick.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3761 by caltinay, Mon Jan 9 06:50:33 2012 UTC revision 3762 by caltinay, Tue Jan 10 00:01:45 2012 UTC
# Line 5329  void Brick::assemblePDESingle(Paso_Syste Line 5329  void Brick::assemblePDESingle(Paso_Syste
5329  }  }
5330    
5331  //protected  //protected
5332    void Brick::assemblePDESingleReduced(Paso_SystemMatrix* mat,
5333            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
5334            const escript::Data& C, const escript::Data& D,
5335            const escript::Data& X, const escript::Data& Y,
5336            const escript::Data& d, const escript::Data& y) const
5337    {
5338        const double h0 = m_l0/m_gNE0;
5339        const double h1 = m_l1/m_gNE1;
5340        const double h2 = m_l2/m_gNE2;
5341        /* GENERATOR SNIP_PDE_SINGLE_REDUCED_PRE TOP */
5342        const double w0 = 0.0625*h1*h2/h0;
5343        const double w1 = 0.0625*h2;
5344        const double w2 = -0.0625*h1;
5345        const double w3 = 0.0625*h0*h2/h1;
5346        const double w4 = -0.0625*h0;
5347        const double w5 = 0.0625*h1;
5348        const double w6 = 0.0625*h0;
5349        const double w7 = -0.0625*h0*h1/h2;
5350        const double w8 = -0.0625*h1*h2/h0;
5351        const double w9 = -0.0625*h2;
5352        const double w10 = -0.0625*h0*h2/h1;
5353        const double w11 = 0.0625*h0*h1/h2;
5354        const double w12 = 0.03125*h1*h2;
5355        const double w13 = 0.03125*h0*h2;
5356        const double w14 = 0.03125*h0*h1;
5357        const double w15 = -0.03125*h1*h2;
5358        const double w16 = -0.03125*h0*h2;
5359        const double w17 = -0.03125*h0*h1;
5360        const double w18 = 0.015625*h0*h1*h2;
5361        const double w19 = -0.25*h1*h2;
5362        const double w20 = -0.25*h0*h2;
5363        const double w21 = -0.25*h0*h1;
5364        const double w22 = 0.25*h1*h2;
5365        const double w23 = 0.25*h0*h2;
5366        const double w24 = 0.25*h0*h1;
5367        const double w25 = 0.125*h0*h1*h2;
5368        /* GENERATOR SNIP_PDE_SINGLE_REDUCED_PRE BOTTOM */
5369    
5370        rhs.requireWrite();
5371    #pragma omp parallel
5372        {
5373            for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
5374    #pragma omp for
5375                for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
5376                    for (index_t k1=0; k1<m_NE1; ++k1) {
5377                        for (index_t k0=0; k0<m_NE0; ++k0)  {
5378                            bool add_EM_S=false;
5379                            bool add_EM_F=false;
5380                            vector<double> EM_S(8*8, 0);
5381                            vector<double> EM_F(8, 0);
5382                            const index_t e = k0 + m_NE0*k1 + m_NE0*m_NE1*k2;
5383                            /* GENERATOR SNIP_PDE_SINGLE_REDUCED TOP */
5384                            ///////////////
5385                            // process A //
5386                            ///////////////
5387                            if (!A.isEmpty()) {
5388                                add_EM_S=true;
5389                                const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
5390                                const register double A_00 = A_p[INDEX2(0,0,3)];
5391                                const register double A_01 = A_p[INDEX2(0,1,3)];
5392                                const register double A_02 = A_p[INDEX2(0,2,3)];
5393                                const register double A_10 = A_p[INDEX2(1,0,3)];
5394                                const register double A_11 = A_p[INDEX2(1,1,3)];
5395                                const register double A_12 = A_p[INDEX2(1,2,3)];
5396                                const register double A_20 = A_p[INDEX2(2,0,3)];
5397                                const register double A_21 = A_p[INDEX2(2,1,3)];
5398                                const register double A_22 = A_p[INDEX2(2,2,3)];
5399                                const register double tmp0_0 = A_01 + A_10;
5400                                const register double tmp1_0 = A_02 + A_20;
5401                                const register double tmp2_0 = A_12 + A_21;
5402                                const register double tmp3_1 = A_22*w7;
5403                                const register double tmp10_1 = A_11*w10;
5404                                const register double tmp21_1 = A_02*w5;
5405                                const register double tmp2_1 = A_00*w0;
5406                                const register double tmp23_1 = tmp2_0*w6;
5407                                const register double tmp19_1 = A_20*w2;
5408                                const register double tmp4_1 = A_11*w3;
5409                                const register double tmp22_1 = tmp1_0*w5;
5410                                const register double tmp13_1 = A_21*w4;
5411                                const register double tmp5_1 = A_21*w6;
5412                                const register double tmp8_1 = A_00*w8;
5413                                const register double tmp7_1 = A_20*w5;
5414                                const register double tmp18_1 = tmp2_0*w4;
5415                                const register double tmp6_1 = A_02*w2;
5416                                const register double tmp9_1 = A_22*w11;
5417                                const register double tmp15_1 = tmp1_0*w2;
5418                                const register double tmp12_1 = A_01*w1;
5419                                const register double tmp0_1 = tmp0_0*w1;
5420                                const register double tmp20_1 = A_01*w9;
5421                                const register double tmp14_1 = A_12*w6;
5422                                const register double tmp1_1 = A_12*w4;
5423                                const register double tmp16_1 = A_10*w9;
5424                                const register double tmp11_1 = tmp0_0*w9;
5425                                const register double tmp17_1 = A_10*w1;
5426                                EM_S[INDEX2(7,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
5427                                EM_S[INDEX2(4,7,8)]+=tmp10_1 + tmp11_1 + tmp1_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
5428                                EM_S[INDEX2(1,3,8)]+=tmp10_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1 + tmp16_1 + tmp2_1 + tmp9_1;
5429                                EM_S[INDEX2(6,4,8)]+=tmp10_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1 + tmp16_1 + tmp2_1 + tmp9_1;
5430                                EM_S[INDEX2(3,0,8)]+=tmp10_1 + tmp11_1 + tmp1_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
5431                                EM_S[INDEX2(5,4,8)]+=tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp4_1 + tmp8_1 + tmp9_1;
5432                                EM_S[INDEX2(0,7,8)]+=tmp10_1 + tmp11_1 + tmp15_1 + tmp18_1 + tmp3_1 + tmp8_1;
5433                                EM_S[INDEX2(5,6,8)]+=tmp0_1 + tmp10_1 + tmp19_1 + tmp1_1 + tmp21_1 + tmp5_1 + tmp8_1 + tmp9_1;
5434                                EM_S[INDEX2(2,6,8)]+=tmp11_1 + tmp13_1 + tmp14_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp6_1 + tmp7_1;
5435                                EM_S[INDEX2(1,6,8)]+=tmp0_1 + tmp10_1 + tmp18_1 + tmp22_1 + tmp3_1 + tmp8_1;
5436                                EM_S[INDEX2(5,1,8)]+=tmp11_1 + tmp13_1 + tmp14_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp6_1 + tmp7_1;
5437                                EM_S[INDEX2(3,7,8)]+=tmp0_1 + tmp13_1 + tmp14_1 + tmp19_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1;
5438                                EM_S[INDEX2(2,5,8)]+=tmp0_1 + tmp10_1 + tmp15_1 + tmp23_1 + tmp3_1 + tmp8_1;
5439                                EM_S[INDEX2(0,3,8)]+=tmp10_1 + tmp11_1 + tmp13_1 + tmp14_1 + tmp19_1 + tmp21_1 + tmp8_1 + tmp9_1;
5440                                EM_S[INDEX2(7,2,8)]+=tmp12_1 + tmp15_1 + tmp16_1 + tmp1_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp8_1;
5441                                EM_S[INDEX2(4,0,8)]+=tmp0_1 + tmp13_1 + tmp14_1 + tmp19_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1;
5442                                EM_S[INDEX2(1,2,8)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp14_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
5443                                EM_S[INDEX2(6,7,8)]+=tmp17_1 + tmp20_1 + tmp23_1 + tmp4_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
5444                                EM_S[INDEX2(3,3,8)]+=tmp0_1 + tmp15_1 + tmp18_1 + tmp2_1 + tmp4_1 + tmp9_1;
5445                                EM_S[INDEX2(2,0,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp1_1 + tmp22_1 + tmp2_1 + tmp5_1 + tmp9_1;
5446                                EM_S[INDEX2(7,6,8)]+=tmp12_1 + tmp16_1 + tmp19_1 + tmp21_1 + tmp23_1 + tmp4_1 + tmp8_1 + tmp9_1;
5447                                EM_S[INDEX2(4,4,8)]+=tmp0_1 + tmp15_1 + tmp18_1 + tmp2_1 + tmp4_1 + tmp9_1;
5448                                EM_S[INDEX2(6,3,8)]+=tmp17_1 + tmp1_1 + tmp20_1 + tmp22_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp8_1;
5449                                EM_S[INDEX2(1,5,8)]+=tmp11_1 + tmp19_1 + tmp1_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
5450                                EM_S[INDEX2(3,6,8)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp16_1 + tmp22_1 + tmp3_1 + tmp4_1 + tmp8_1;
5451                                EM_S[INDEX2(2,2,8)]+=tmp11_1 + tmp18_1 + tmp22_1 + tmp2_1 + tmp4_1 + tmp9_1;
5452                                EM_S[INDEX2(7,7,8)]+=tmp0_1 + tmp22_1 + tmp23_1 + tmp2_1 + tmp4_1 + tmp9_1;
5453                                EM_S[INDEX2(5,7,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp1_1 + tmp22_1 + tmp2_1 + tmp5_1 + tmp9_1;
5454                                EM_S[INDEX2(5,3,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp23_1 + tmp2_1 + tmp3_1 + tmp6_1 + tmp7_1;
5455                                EM_S[INDEX2(4,1,8)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp16_1 + tmp22_1 + tmp3_1 + tmp4_1 + tmp8_1;
5456                                EM_S[INDEX2(1,1,8)]+=tmp11_1 + tmp15_1 + tmp23_1 + tmp2_1 + tmp4_1 + tmp9_1;
5457                                EM_S[INDEX2(2,7,8)]+=tmp13_1 + tmp14_1 + tmp15_1 + tmp17_1 + tmp20_1 + tmp3_1 + tmp4_1 + tmp8_1;
5458                                EM_S[INDEX2(3,2,8)]+=tmp12_1 + tmp16_1 + tmp18_1 + tmp4_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
5459                                EM_S[INDEX2(0,0,8)]+=tmp0_1 + tmp22_1 + tmp23_1 + tmp2_1 + tmp4_1 + tmp9_1;
5460                                EM_S[INDEX2(6,6,8)]+=tmp11_1 + tmp15_1 + tmp23_1 + tmp2_1 + tmp4_1 + tmp9_1;
5461                                EM_S[INDEX2(5,0,8)]+=tmp13_1 + tmp14_1 + tmp15_1 + tmp17_1 + tmp20_1 + tmp3_1 + tmp4_1 + tmp8_1;
5462                                EM_S[INDEX2(7,1,8)]+=tmp10_1 + tmp17_1 + tmp18_1 + tmp20_1 + tmp2_1 + tmp3_1 + tmp6_1 + tmp7_1;
5463                                EM_S[INDEX2(4,5,8)]+=tmp12_1 + tmp16_1 + tmp18_1 + tmp4_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
5464                                EM_S[INDEX2(0,4,8)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
5465                                EM_S[INDEX2(5,5,8)]+=tmp11_1 + tmp18_1 + tmp22_1 + tmp2_1 + tmp4_1 + tmp9_1;
5466                                EM_S[INDEX2(1,4,8)]+=tmp17_1 + tmp1_1 + tmp20_1 + tmp22_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp8_1;
5467                                EM_S[INDEX2(6,0,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp18_1 + tmp19_1 + tmp21_1 + tmp2_1 + tmp3_1;
5468                                EM_S[INDEX2(7,5,8)]+=tmp10_1 + tmp13_1 + tmp14_1 + tmp17_1 + tmp20_1 + tmp22_1 + tmp2_1 + tmp9_1;
5469                                EM_S[INDEX2(2,3,8)]+=tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp4_1 + tmp8_1 + tmp9_1;
5470                                EM_S[INDEX2(2,1,8)]+=tmp0_1 + tmp10_1 + tmp19_1 + tmp1_1 + tmp21_1 + tmp5_1 + tmp8_1 + tmp9_1;
5471                                EM_S[INDEX2(4,2,8)]+=tmp10_1 + tmp17_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp23_1 + tmp2_1 + tmp3_1;
5472                                EM_S[INDEX2(1,0,8)]+=tmp17_1 + tmp20_1 + tmp23_1 + tmp4_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
5473                                EM_S[INDEX2(6,5,8)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp14_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
5474                                EM_S[INDEX2(3,5,8)]+=tmp10_1 + tmp17_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp23_1 + tmp2_1 + tmp3_1;
5475                                EM_S[INDEX2(0,1,8)]+=tmp12_1 + tmp16_1 + tmp19_1 + tmp21_1 + tmp23_1 + tmp4_1 + tmp8_1 + tmp9_1;
5476                                EM_S[INDEX2(7,0,8)]+=tmp10_1 + tmp11_1 + tmp15_1 + tmp18_1 + tmp3_1 + tmp8_1;
5477                                EM_S[INDEX2(4,6,8)]+=tmp10_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp20_1 + tmp2_1 + tmp5_1 + tmp9_1;
5478                                EM_S[INDEX2(5,2,8)]+=tmp0_1 + tmp10_1 + tmp15_1 + tmp23_1 + tmp3_1 + tmp8_1;
5479                                EM_S[INDEX2(6,1,8)]+=tmp0_1 + tmp10_1 + tmp18_1 + tmp22_1 + tmp3_1 + tmp8_1;
5480                                EM_S[INDEX2(3,1,8)]+=tmp10_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp20_1 + tmp2_1 + tmp5_1 + tmp9_1;
5481                                EM_S[INDEX2(0,2,8)]+=tmp10_1 + tmp13_1 + tmp14_1 + tmp17_1 + tmp20_1 + tmp22_1 + tmp2_1 + tmp9_1;
5482                                EM_S[INDEX2(7,4,8)]+=tmp10_1 + tmp11_1 + tmp13_1 + tmp14_1 + tmp19_1 + tmp21_1 + tmp8_1 + tmp9_1;
5483                                EM_S[INDEX2(0,6,8)]+=tmp10_1 + tmp17_1 + tmp18_1 + tmp20_1 + tmp2_1 + tmp3_1 + tmp6_1 + tmp7_1;
5484                                EM_S[INDEX2(6,2,8)]+=tmp11_1 + tmp19_1 + tmp1_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
5485                                EM_S[INDEX2(4,3,8)]+=tmp10_1 + tmp11_1 + tmp22_1 + tmp23_1 + tmp3_1 + tmp8_1;
5486                                EM_S[INDEX2(1,7,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp18_1 + tmp19_1 + tmp21_1 + tmp2_1 + tmp3_1;
5487                                EM_S[INDEX2(0,5,8)]+=tmp12_1 + tmp15_1 + tmp16_1 + tmp1_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp8_1;
5488                                EM_S[INDEX2(3,4,8)]+=tmp10_1 + tmp11_1 + tmp22_1 + tmp23_1 + tmp3_1 + tmp8_1;
5489                                EM_S[INDEX2(2,4,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp23_1 + tmp2_1 + tmp3_1 + tmp6_1 + tmp7_1;
5490                            }
5491                            ///////////////
5492                            // process B //
5493                            ///////////////
5494                            if (!B.isEmpty()) {
5495                                add_EM_S=true;
5496                                const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
5497                                const register double B_0 = B_p[0];
5498                                const register double B_1 = B_p[1];
5499                                const register double B_2 = B_p[2];
5500                                const register double tmp4_1 = B_0*w15;
5501                                const register double tmp3_1 = B_1*w16;
5502                                const register double tmp2_1 = B_0*w12;
5503                                const register double tmp5_1 = B_2*w17;
5504                                const register double tmp1_1 = B_2*w14;
5505                                const register double tmp0_1 = B_1*w13;
5506                                EM_S[INDEX2(7,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5507                                EM_S[INDEX2(4,7,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
5508                                EM_S[INDEX2(1,3,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
5509                                EM_S[INDEX2(6,4,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
5510                                EM_S[INDEX2(3,0,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
5511                                EM_S[INDEX2(5,4,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
5512                                EM_S[INDEX2(0,7,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5513                                EM_S[INDEX2(5,6,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
5514                                EM_S[INDEX2(2,6,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
5515                                EM_S[INDEX2(1,6,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
5516                                EM_S[INDEX2(5,1,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
5517                                EM_S[INDEX2(3,7,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
5518                                EM_S[INDEX2(2,5,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
5519                                EM_S[INDEX2(0,3,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5520                                EM_S[INDEX2(7,2,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5521                                EM_S[INDEX2(4,0,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
5522                                EM_S[INDEX2(1,2,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
5523                                EM_S[INDEX2(6,7,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
5524                                EM_S[INDEX2(3,3,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
5525                                EM_S[INDEX2(2,0,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
5526                                EM_S[INDEX2(7,6,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5527                                EM_S[INDEX2(4,4,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
5528                                EM_S[INDEX2(6,3,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
5529                                EM_S[INDEX2(1,5,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
5530                                EM_S[INDEX2(3,6,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
5531                                EM_S[INDEX2(2,2,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
5532                                EM_S[INDEX2(7,7,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5533                                EM_S[INDEX2(5,7,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
5534                                EM_S[INDEX2(5,3,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
5535                                EM_S[INDEX2(4,1,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
5536                                EM_S[INDEX2(1,1,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
5537                                EM_S[INDEX2(2,7,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
5538                                EM_S[INDEX2(3,2,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
5539                                EM_S[INDEX2(0,0,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5540                                EM_S[INDEX2(6,6,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
5541                                EM_S[INDEX2(5,0,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
5542                                EM_S[INDEX2(7,1,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5543                                EM_S[INDEX2(4,5,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
5544                                EM_S[INDEX2(0,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5545                                EM_S[INDEX2(5,5,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
5546                                EM_S[INDEX2(1,4,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
5547                                EM_S[INDEX2(6,0,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
5548                                EM_S[INDEX2(7,5,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5549                                EM_S[INDEX2(2,3,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
5550                                EM_S[INDEX2(2,1,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
5551                                EM_S[INDEX2(4,2,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
5552                                EM_S[INDEX2(1,0,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
5553                                EM_S[INDEX2(6,5,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
5554                                EM_S[INDEX2(3,5,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
5555                                EM_S[INDEX2(0,1,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5556                                EM_S[INDEX2(7,0,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5557                                EM_S[INDEX2(4,6,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
5558                                EM_S[INDEX2(5,2,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
5559                                EM_S[INDEX2(6,1,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
5560                                EM_S[INDEX2(3,1,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
5561                                EM_S[INDEX2(0,2,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5562                                EM_S[INDEX2(7,4,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5563                                EM_S[INDEX2(0,6,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5564                                EM_S[INDEX2(6,2,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
5565                                EM_S[INDEX2(4,3,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
5566                                EM_S[INDEX2(1,7,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
5567                                EM_S[INDEX2(0,5,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5568                                EM_S[INDEX2(3,4,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
5569                                EM_S[INDEX2(2,4,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
5570                            }
5571                            ///////////////
5572                            // process C //
5573                            ///////////////
5574                            if (!C.isEmpty()) {
5575                                add_EM_S=true;
5576                                const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
5577                                const register double C_0 = C_p[0];
5578                                const register double C_1 = C_p[1];
5579                                const register double C_2 = C_p[2];
5580                                const register double tmp5_1 = C_0*w15;
5581                                const register double tmp2_1 = C_0*w12;
5582                                const register double tmp4_1 = C_1*w16;
5583                                const register double tmp1_1 = C_2*w17;
5584                                const register double tmp3_1 = C_2*w14;
5585                                const register double tmp0_1 = C_1*w13;
5586                                EM_S[INDEX2(7,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5587                                EM_S[INDEX2(4,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
5588                                EM_S[INDEX2(1,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5589                                EM_S[INDEX2(6,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5590                                EM_S[INDEX2(3,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
5591                                EM_S[INDEX2(5,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5592                                EM_S[INDEX2(0,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
5593                                EM_S[INDEX2(5,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
5594                                EM_S[INDEX2(2,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
5595                                EM_S[INDEX2(1,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
5596                                EM_S[INDEX2(5,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
5597                                EM_S[INDEX2(3,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
5598                                EM_S[INDEX2(2,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
5599                                EM_S[INDEX2(0,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5600                                EM_S[INDEX2(7,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
5601                                EM_S[INDEX2(4,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
5602                                EM_S[INDEX2(1,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
5603                                EM_S[INDEX2(6,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
5604                                EM_S[INDEX2(3,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5605                                EM_S[INDEX2(2,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
5606                                EM_S[INDEX2(7,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
5607                                EM_S[INDEX2(4,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5608                                EM_S[INDEX2(6,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5609                                EM_S[INDEX2(1,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
5610                                EM_S[INDEX2(3,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
5611                                EM_S[INDEX2(2,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
5612                                EM_S[INDEX2(7,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
5613                                EM_S[INDEX2(5,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
5614                                EM_S[INDEX2(5,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5615                                EM_S[INDEX2(4,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
5616                                EM_S[INDEX2(1,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
5617                                EM_S[INDEX2(2,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
5618                                EM_S[INDEX2(3,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
5619                                EM_S[INDEX2(0,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
5620                                EM_S[INDEX2(6,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
5621                                EM_S[INDEX2(5,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
5622                                EM_S[INDEX2(7,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
5623                                EM_S[INDEX2(4,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
5624                                EM_S[INDEX2(0,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5625                                EM_S[INDEX2(5,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
5626                                EM_S[INDEX2(1,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5627                                EM_S[INDEX2(6,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
5628                                EM_S[INDEX2(7,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
5629                                EM_S[INDEX2(2,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5630                                EM_S[INDEX2(2,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
5631                                EM_S[INDEX2(4,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
5632                                EM_S[INDEX2(1,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
5633                                EM_S[INDEX2(6,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
5634                                EM_S[INDEX2(3,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
5635                                EM_S[INDEX2(0,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
5636                                EM_S[INDEX2(7,0,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
5637                                EM_S[INDEX2(4,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
5638                                EM_S[INDEX2(5,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
5639                                EM_S[INDEX2(6,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
5640                                EM_S[INDEX2(3,1,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
5641                                EM_S[INDEX2(0,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
5642                                EM_S[INDEX2(7,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5643                                EM_S[INDEX2(0,6,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
5644                                EM_S[INDEX2(6,2,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
5645                                EM_S[INDEX2(4,3,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
5646                                EM_S[INDEX2(1,7,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
5647                                EM_S[INDEX2(0,5,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
5648                                EM_S[INDEX2(3,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5649                                EM_S[INDEX2(2,4,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
5650                            }
5651                            ///////////////
5652                            // process D //
5653                            ///////////////
5654                            if (!D.isEmpty()) {
5655                                add_EM_S=true;
5656                                const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
5657                                const register double D_0 = D_p[0];
5658                                const register double tmp0_1 = D_0*w18;
5659                                EM_S[INDEX2(7,3,8)]+=tmp0_1;
5660                                EM_S[INDEX2(4,7,8)]+=tmp0_1;
5661                                EM_S[INDEX2(1,3,8)]+=tmp0_1;
5662                                EM_S[INDEX2(6,4,8)]+=tmp0_1;
5663                                EM_S[INDEX2(3,0,8)]+=tmp0_1;
5664                                EM_S[INDEX2(5,4,8)]+=tmp0_1;
5665                                EM_S[INDEX2(0,7,8)]+=tmp0_1;
5666                                EM_S[INDEX2(5,6,8)]+=tmp0_1;
5667                                EM_S[INDEX2(2,6,8)]+=tmp0_1;
5668                                EM_S[INDEX2(1,6,8)]+=tmp0_1;
5669                                EM_S[INDEX2(5,1,8)]+=tmp0_1;
5670                                EM_S[INDEX2(3,7,8)]+=tmp0_1;
5671                                EM_S[INDEX2(2,5,8)]+=tmp0_1;
5672                                EM_S[INDEX2(0,3,8)]+=tmp0_1;
5673                                EM_S[INDEX2(7,2,8)]+=tmp0_1;
5674                                EM_S[INDEX2(4,0,8)]+=tmp0_1;
5675                                EM_S[INDEX2(1,2,8)]+=tmp0_1;
5676                                EM_S[INDEX2(6,7,8)]+=tmp0_1;
5677                                EM_S[INDEX2(3,3,8)]+=tmp0_1;
5678                                EM_S[INDEX2(2,0,8)]+=tmp0_1;
5679                                EM_S[INDEX2(7,6,8)]+=tmp0_1;
5680                                EM_S[INDEX2(4,4,8)]+=tmp0_1;
5681                                EM_S[INDEX2(6,3,8)]+=tmp0_1;
5682                                EM_S[INDEX2(1,5,8)]+=tmp0_1;
5683                                EM_S[INDEX2(3,6,8)]+=tmp0_1;
5684                                EM_S[INDEX2(2,2,8)]+=tmp0_1;
5685                                EM_S[INDEX2(7,7,8)]+=tmp0_1;
5686                                EM_S[INDEX2(5,7,8)]+=tmp0_1;
5687                                EM_S[INDEX2(5,3,8)]+=tmp0_1;
5688                                EM_S[INDEX2(4,1,8)]+=tmp0_1;
5689                                EM_S[INDEX2(1,1,8)]+=tmp0_1;
5690                                EM_S[INDEX2(2,7,8)]+=tmp0_1;
5691                                EM_S[INDEX2(3,2,8)]+=tmp0_1;
5692                                EM_S[INDEX2(0,0,8)]+=tmp0_1;
5693                                EM_S[INDEX2(6,6,8)]+=tmp0_1;
5694                                EM_S[INDEX2(5,0,8)]+=tmp0_1;
5695                                EM_S[INDEX2(7,1,8)]+=tmp0_1;
5696                                EM_S[INDEX2(4,5,8)]+=tmp0_1;
5697                                EM_S[INDEX2(0,4,8)]+=tmp0_1;
5698                                EM_S[INDEX2(5,5,8)]+=tmp0_1;
5699                                EM_S[INDEX2(1,4,8)]+=tmp0_1;
5700                                EM_S[INDEX2(6,0,8)]+=tmp0_1;
5701                                EM_S[INDEX2(7,5,8)]+=tmp0_1;
5702                                EM_S[INDEX2(2,3,8)]+=tmp0_1;
5703                                EM_S[INDEX2(2,1,8)]+=tmp0_1;
5704                                EM_S[INDEX2(4,2,8)]+=tmp0_1;
5705                                EM_S[INDEX2(1,0,8)]+=tmp0_1;
5706                                EM_S[INDEX2(6,5,8)]+=tmp0_1;
5707                                EM_S[INDEX2(3,5,8)]+=tmp0_1;
5708                                EM_S[INDEX2(0,1,8)]+=tmp0_1;
5709                                EM_S[INDEX2(7,0,8)]+=tmp0_1;
5710                                EM_S[INDEX2(4,6,8)]+=tmp0_1;
5711                                EM_S[INDEX2(5,2,8)]+=tmp0_1;
5712                                EM_S[INDEX2(6,1,8)]+=tmp0_1;
5713                                EM_S[INDEX2(3,1,8)]+=tmp0_1;
5714                                EM_S[INDEX2(0,2,8)]+=tmp0_1;
5715                                EM_S[INDEX2(7,4,8)]+=tmp0_1;
5716                                EM_S[INDEX2(0,6,8)]+=tmp0_1;
5717                                EM_S[INDEX2(6,2,8)]+=tmp0_1;
5718                                EM_S[INDEX2(4,3,8)]+=tmp0_1;
5719                                EM_S[INDEX2(1,7,8)]+=tmp0_1;
5720                                EM_S[INDEX2(0,5,8)]+=tmp0_1;
5721                                EM_S[INDEX2(3,4,8)]+=tmp0_1;
5722                                EM_S[INDEX2(2,4,8)]+=tmp0_1;
5723                            }
5724                            ///////////////
5725                            // process X //
5726                            ///////////////
5727                            if (!X.isEmpty()) {
5728                                add_EM_F=true;
5729                                const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
5730                                const register double X_0 = X_p[0];
5731                                const register double X_1 = X_p[1];
5732                                const register double X_2 = X_p[2];
5733                                const register double tmp1_1 = X_0*w19;
5734                                const register double tmp2_1 = X_1*w20;
5735                                const register double tmp3_1 = X_0*w22;
5736                                const register double tmp4_1 = X_1*w23;
5737                                const register double tmp5_1 = X_2*w24;
5738                                const register double tmp0_1 = X_2*w21;
5739                                EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
5740                                EM_F[1]+=tmp0_1 + tmp2_1 + tmp3_1;
5741                                EM_F[2]+=tmp0_1 + tmp1_1 + tmp4_1;
5742                                EM_F[3]+=tmp0_1 + tmp3_1 + tmp4_1;
5743                                EM_F[4]+=tmp1_1 + tmp2_1 + tmp5_1;
5744                                EM_F[5]+=tmp2_1 + tmp3_1 + tmp5_1;
5745                                EM_F[6]+=tmp1_1 + tmp4_1 + tmp5_1;
5746                                EM_F[7]+=tmp3_1 + tmp4_1 + tmp5_1;
5747                            }
5748                            ///////////////
5749                            // process Y //
5750                            ///////////////
5751                            if (!Y.isEmpty()) {
5752                                add_EM_F=true;
5753                                const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
5754                                const register double Y_0 = Y_p[0];
5755                                const register double tmp0_1 = Y_0*w25;
5756                                EM_F[0]+=tmp0_1;
5757                                EM_F[1]+=tmp0_1;
5758                                EM_F[2]+=tmp0_1;
5759                                EM_F[3]+=tmp0_1;
5760                                EM_F[4]+=tmp0_1;
5761                                EM_F[5]+=tmp0_1;
5762                                EM_F[6]+=tmp0_1;
5763                                EM_F[7]+=tmp0_1;
5764                            }
5765                            /* GENERATOR SNIP_PDE_SINGLE_REDUCED BOTTOM */
5766    
5767                            // add to matrix (if add_EM_S) and RHS (if add_EM_F)
5768                            const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;
5769                            IndexVector rowIndex;
5770                            rowIndex.push_back(m_dofMap[firstNode]);
5771                            rowIndex.push_back(m_dofMap[firstNode+1]);
5772                            rowIndex.push_back(m_dofMap[firstNode+m_N0]);
5773                            rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
5774                            rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1]);
5775                            rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1+1]);
5776                            rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);
5777                            rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);
5778                            if (add_EM_F) {
5779                                //cout << "-- AddtoRHS -- " << endl;
5780                                double *F_p=rhs.getSampleDataRW(0);
5781                                for (index_t i=0; i<rowIndex.size(); i++) {
5782                                    if (rowIndex[i]<getNumDOF()) {
5783                                        F_p[rowIndex[i]]+=EM_F[i];
5784                                        //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;
5785                                    }
5786                                }
5787                                //cout << "---"<<endl;
5788                            }
5789                            if (add_EM_S) {
5790                                //cout << "-- AddtoSystem -- " << endl;
5791                                addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
5792                            }
5793                        } // end k0 loop
5794                    } // end k1 loop
5795                } // end k2 loop
5796            } // end of colouring
5797        } // end of parallel region
5798    }
5799    
5800    //protected
5801  void Brick::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,  void Brick::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
5802          const escript::Data& A, const escript::Data& B,          const escript::Data& A, const escript::Data& B,
5803          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
# Line 5346  void Brick::assemblePDESystem(Paso_Syste Line 5815  void Brick::assemblePDESystem(Paso_Syste
5815          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
5816      }      }
5817    
5818      /* GENERATOR SNIP_PDE_SYSTEM_PRE TOP */      /*** GENERATOR SNIP_PDE_SYSTEM_PRE TOP */
5819      const double w0 = 0.0009303791403858427308*h1*h2/h0;      const double w0 = 0.0009303791403858427308*h1*h2/h0;
5820      const double w1 = 0.0009303791403858427308*h2;      const double w1 = 0.0009303791403858427308*h2;
5821      const double w10 = 0.012958509748503046158*h0*h2/h1;      const double w10 = 0.012958509748503046158*h0*h2/h1;
# Line 5546  void Brick::assemblePDESystem(Paso_Syste Line 6015  void Brick::assemblePDESystem(Paso_Syste
6015                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
6016                          vector<double> EM_F(8*numEq, 0);                          vector<double> EM_F(8*numEq, 0);
6017                          const index_t e = k0 + m_NE0*k1 + m_NE0*m_NE1*k2;                          const index_t e = k0 + m_NE0*k1 + m_NE0*m_NE1*k2;
6018                          /* GENERATOR SNIP_PDE_SYSTEM TOP */                          /*** GENERATOR SNIP_PDE_SYSTEM TOP */
6019                          ///////////////                          ///////////////
6020                          // process A //                          // process A //
6021                          ///////////////                          ///////////////
# Line 8472  void Brick::assemblePDESystem(Paso_Syste Line 8941  void Brick::assemblePDESystem(Paso_Syste
8941    
8942                          // add to matrix (if add_EM_S) and RHS (if add_EM_F)                          // add to matrix (if add_EM_S) and RHS (if add_EM_F)
8943                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;                          const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;
8944                            IndexVector rowIndex;
8945                            rowIndex.push_back(m_dofMap[firstNode]);
8946                            rowIndex.push_back(m_dofMap[firstNode+1]);
8947                            rowIndex.push_back(m_dofMap[firstNode+m_N0]);
8948                            rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
8949                            rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1]);
8950                            rowIndex.push_back(m_dofMap[firstNode+m_N0*m_N1+1]);
8951                            rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)]);
8952                            rowIndex.push_back(m_dofMap[firstNode+m_N0*(m_N1+1)+1]);
8953                            if (add_EM_F) {
8954                                //cout << "-- AddtoRHS -- " << endl;
8955                                double *F_p=rhs.getSampleDataRW(0);
8956                                for (index_t i=0; i<rowIndex.size(); i++) {
8957                                    if (rowIndex[i]<getNumDOF()) {
8958                                        for (index_t eq=0; eq<numEq; eq++) {
8959                                            F_p[INDEX2(eq,rowIndex[i],numEq)]+=EM_F[INDEX2(eq,i,numEq)];
8960                                            //cout << "F[" << INDEX2(eq,rowIndex[i],numEq) << "]=" << F_p[INDEX2(eq,rowIndex[i],numEq)] << endl;
8961                                        }
8962                                    }
8963                                }
8964                                //cout << "---"<<endl;
8965                            }
8966                            if (add_EM_S) {
8967                                //cout << "-- AddtoSystem -- " << endl;
8968                                addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);
8969                            }
8970                        } // end k0 loop
8971                    } // end k1 loop
8972                } // end k2 loop
8973            } // end of colouring
8974        } // end of parallel region
8975    }
8976    
8977    //protected
8978    void Brick::assemblePDESystemReduced(Paso_SystemMatrix* mat,
8979            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
8980            const escript::Data& C, const escript::Data& D,
8981            const escript::Data& X, const escript::Data& Y,
8982            const escript::Data& d, const escript::Data& y) const
8983    {
8984        const double h0 = m_l0/m_gNE0;
8985        const double h1 = m_l1/m_gNE1;
8986        const double h2 = m_l2/m_gNE2;
8987        dim_t numEq, numComp;
8988        if (!mat)
8989            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
8990        else {
8991            numEq=mat->logical_row_block_size;
8992            numComp=mat->logical_col_block_size;
8993        }
8994    
8995        /* GENERATOR SNIP_PDE_SYSTEM_REDUCED_PRE TOP */
8996        const double w0 = 0.0625*h1*h2/h0;
8997        const double w1 = 0.0625*h2;
8998        const double w10 = -0.0625*h0*h2/h1;
8999        const double w11 = 0.0625*h0*h1/h2;
9000        const double w12 = 0.03125*h1*h2;
9001        const double w13 = 0.03125*h0*h2;
9002        const double w14 = 0.03125*h0*h1;
9003        const double w15 = -0.03125*h1*h2;
9004        const double w16 = -0.03125*h0*h2;
9005        const double w17 = -0.03125*h0*h1;
9006        const double w18 = 0.015625*h0*h1*h2;
9007        const double w19 = -0.25*h1*h2;
9008        const double w2 = -0.0625*h1;
9009        const double w20 = -0.25*h0*h2;
9010        const double w21 = -0.25*h0*h1;
9011        const double w22 = 0.25*h1*h2;
9012        const double w23 = 0.25*h0*h2;
9013        const double w24 = 0.25*h0*h1;
9014        const double w25 = 0.125*h0*h1*h2;
9015        const double w3 = 0.0625*h0*h2/h1;
9016        const double w4 = -0.0625*h0;
9017        const double w5 = 0.0625*h1;
9018        const double w6 = 0.0625*h0;
9019        const double w7 = -0.0625*h0*h1/h2;
9020        const double w8 = -0.0625*h1*h2/h0;
9021        const double w9 = -0.0625*h2;
9022        /* GENERATOR SNIP_PDE_SYSTEM_REDUCED_PRE BOTTOM */
9023    
9024        rhs.requireWrite();
9025    #pragma omp parallel
9026        {
9027            for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
9028    #pragma omp for
9029                for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
9030                    for (index_t k1=0; k1<m_NE1; ++k1) {
9031                        for (index_t k0=0; k0<m_NE0; ++k0)  {
9032                            bool add_EM_S=false;
9033                            bool add_EM_F=false;
9034                            vector<double> EM_S(8*8*numEq*numComp, 0);
9035                            vector<double> EM_F(8*numEq, 0);
9036                            const index_t e = k0 + m_NE0*k1 + m_NE0*m_NE1*k2;
9037                            /* GENERATOR SNIP_PDE_SYSTEM_REDUCED TOP */
9038                            ///////////////
9039                            // process A //
9040                            ///////////////
9041                            if (!A.isEmpty()) {
9042                                add_EM_S=true;
9043                                const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
9044                                for (index_t k=0; k<numEq; k++) {
9045                                    for (index_t m=0; m<numComp; m++) {
9046                                        const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,3, numComp)];
9047                                        const register double A_01 = A_p[INDEX4(k,0,m,1, numEq,3, numComp)];
9048                                        const register double A_02 = A_p[INDEX4(k,0,m,2, numEq,3, numComp)];
9049                                        const register double A_10 = A_p[INDEX4(k,1,m,0, numEq,3, numComp)];
9050                                        const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,3, numComp)];
9051                                        const register double A_12 = A_p[INDEX4(k,1,m,2, numEq,3, numComp)];
9052                                        const register double A_20 = A_p[INDEX4(k,2,m,0, numEq,3, numComp)];
9053                                        const register double A_21 = A_p[INDEX4(k,2,m,1, numEq,3, numComp)];
9054                                        const register double A_22 = A_p[INDEX4(k,2,m,2, numEq,3, numComp)];
9055                                        const register double tmp0_0 = A_01 + A_10;
9056                                        const register double tmp1_0 = A_02 + A_20;
9057                                        const register double tmp2_0 = A_12 + A_21;
9058                                        const register double tmp3_1 = A_22*w7;
9059                                        const register double tmp10_1 = A_11*w10;
9060                                        const register double tmp21_1 = A_02*w5;
9061                                        const register double tmp2_1 = A_00*w0;
9062                                        const register double tmp23_1 = tmp2_0*w6;
9063                                        const register double tmp19_1 = A_20*w2;
9064                                        const register double tmp4_1 = A_11*w3;
9065                                        const register double tmp22_1 = tmp1_0*w5;
9066                                        const register double tmp13_1 = A_21*w4;
9067                                        const register double tmp5_1 = A_21*w6;
9068                                        const register double tmp8_1 = A_00*w8;
9069                                        const register double tmp7_1 = A_20*w5;
9070                                        const register double tmp18_1 = tmp2_0*w4;
9071                                        const register double tmp6_1 = A_02*w2;
9072                                        const register double tmp9_1 = A_22*w11;
9073                                        const register double tmp15_1 = tmp1_0*w2;
9074                                        const register double tmp12_1 = A_01*w1;
9075                                        const register double tmp0_1 = tmp0_0*w1;
9076                                        const register double tmp20_1 = A_01*w9;
9077                                        const register double tmp14_1 = A_12*w6;
9078                                        const register double tmp1_1 = A_12*w4;
9079                                        const register double tmp16_1 = A_10*w9;
9080                                        const register double tmp11_1 = tmp0_0*w9;
9081                                        const register double tmp17_1 = A_10*w1;
9082                                        EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
9083                                        EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp1_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
9084                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp10_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1 + tmp16_1 + tmp2_1 + tmp9_1;
9085                                        EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp10_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1 + tmp16_1 + tmp2_1 + tmp9_1;
9086                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp1_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
9087                                        EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp4_1 + tmp8_1 + tmp9_1;
9088                                        EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp15_1 + tmp18_1 + tmp3_1 + tmp8_1;
9089                                        EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0_1 + tmp10_1 + tmp19_1 + tmp1_1 + tmp21_1 + tmp5_1 + tmp8_1 + tmp9_1;
9090                                        EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp11_1 + tmp13_1 + tmp14_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp6_1 + tmp7_1;
9091                                        EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=tmp0_1 + tmp10_1 + tmp18_1 + tmp22_1 + tmp3_1 + tmp8_1;
9092                                        EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp11_1 + tmp13_1 + tmp14_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp6_1 + tmp7_1;
9093                                        EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1 + tmp13_1 + tmp14_1 + tmp19_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1;
9094                                        EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=tmp0_1 + tmp10_1 + tmp15_1 + tmp23_1 + tmp3_1 + tmp8_1;
9095                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp13_1 + tmp14_1 + tmp19_1 + tmp21_1 + tmp8_1 + tmp9_1;
9096                                        EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp12_1 + tmp15_1 + tmp16_1 + tmp1_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp8_1;
9097                                        EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0_1 + tmp13_1 + tmp14_1 + tmp19_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1;
9098                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp14_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
9099                                        EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp17_1 + tmp20_1 + tmp23_1 + tmp4_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
9100                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0_1 + tmp15_1 + tmp18_1 + tmp2_1 + tmp4_1 + tmp9_1;
9101                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp1_1 + tmp22_1 + tmp2_1 + tmp5_1 + tmp9_1;
9102                                        EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp12_1 + tmp16_1 + tmp19_1 + tmp21_1 + tmp23_1 + tmp4_1 + tmp8_1 + tmp9_1;
9103                                        EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0_1 + tmp15_1 + tmp18_1 + tmp2_1 + tmp4_1 + tmp9_1;
9104                                        EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp17_1 + tmp1_1 + tmp20_1 + tmp22_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp8_1;
9105                                        EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp11_1 + tmp19_1 + tmp1_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
9106                                        EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp16_1 + tmp22_1 + tmp3_1 + tmp4_1 + tmp8_1;
9107                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp11_1 + tmp18_1 + tmp22_1 + tmp2_1 + tmp4_1 + tmp9_1;
9108                                        EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0_1 + tmp22_1 + tmp23_1 + tmp2_1 + tmp4_1 + tmp9_1;
9109                                        EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp1_1 + tmp22_1 + tmp2_1 + tmp5_1 + tmp9_1;
9110                                        EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp23_1 + tmp2_1 + tmp3_1 + tmp6_1 + tmp7_1;
9111                                        EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp16_1 + tmp22_1 + tmp3_1 + tmp4_1 + tmp8_1;
9112                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp11_1 + tmp15_1 + tmp23_1 + tmp2_1 + tmp4_1 + tmp9_1;
9113                                        EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp13_1 + tmp14_1 + tmp15_1 + tmp17_1 + tmp20_1 + tmp3_1 + tmp4_1 + tmp8_1;
9114                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp12_1 + tmp16_1 + tmp18_1 + tmp4_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
9115                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0_1 + tmp22_1 + tmp23_1 + tmp2_1 + tmp4_1 + tmp9_1;
9116                                        EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp11_1 + tmp15_1 + tmp23_1 + tmp2_1 + tmp4_1 + tmp9_1;
9117                                        EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp13_1 + tmp14_1 + tmp15_1 + tmp17_1 + tmp20_1 + tmp3_1 + tmp4_1 + tmp8_1;
9118                                        EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp10_1 + tmp17_1 + tmp18_1 + tmp20_1 + tmp2_1 + tmp3_1 + tmp6_1 + tmp7_1;
9119                                        EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp12_1 + tmp16_1 + tmp18_1 + tmp4_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
9120                                        EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
9121                                        EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp11_1 + tmp18_1 + tmp22_1 + tmp2_1 + tmp4_1 + tmp9_1;
9122                                        EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp17_1 + tmp1_1 + tmp20_1 + tmp22_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp8_1;
9123                                        EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp18_1 + tmp19_1 + tmp21_1 + tmp2_1 + tmp3_1;
9124                                        EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp10_1 + tmp13_1 + tmp14_1 + tmp17_1 + tmp20_1 + tmp22_1 + tmp2_1 + tmp9_1;
9125                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp4_1 + tmp8_1 + tmp9_1;
9126                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp0_1 + tmp10_1 + tmp19_1 + tmp1_1 + tmp21_1 + tmp5_1 + tmp8_1 + tmp9_1;
9127                                        EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp10_1 + tmp17_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp23_1 + tmp2_1 + tmp3_1;
9128                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp17_1 + tmp20_1 + tmp23_1 + tmp4_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
9129                                        EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp14_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
9130                                        EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp10_1 + tmp17_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp23_1 + tmp2_1 + tmp3_1;
9131                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp12_1 + tmp16_1 + tmp19_1 + tmp21_1 + tmp23_1 + tmp4_1 + tmp8_1 + tmp9_1;
9132                                        EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp15_1 + tmp18_1 + tmp3_1 + tmp8_1;
9133                                        EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp10_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp20_1 + tmp2_1 + tmp5_1 + tmp9_1;
9134                                        EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=tmp0_1 + tmp10_1 + tmp15_1 + tmp23_1 + tmp3_1 + tmp8_1;
9135                                        EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=tmp0_1 + tmp10_1 + tmp18_1 + tmp22_1 + tmp3_1 + tmp8_1;
9136                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp10_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp20_1 + tmp2_1 + tmp5_1 + tmp9_1;
9137                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp10_1 + tmp13_1 + tmp14_1 + tmp17_1 + tmp20_1 + tmp22_1 + tmp2_1 + tmp9_1;
9138                                        EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp13_1 + tmp14_1 + tmp19_1 + tmp21_1 + tmp8_1 + tmp9_1;
9139                                        EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp10_1 + tmp17_1 + tmp18_1 + tmp20_1 + tmp2_1 + tmp3_1 + tmp6_1 + tmp7_1;
9140                                        EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp11_1 + tmp19_1 + tmp1_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
9141                                        EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp22_1 + tmp23_1 + tmp3_1 + tmp8_1;
9142                                        EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp18_1 + tmp19_1 + tmp21_1 + tmp2_1 + tmp3_1;
9143                                        EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp12_1 + tmp15_1 + tmp16_1 + tmp1_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp8_1;
9144                                        EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=tmp10_1 + tmp11_1 + tmp22_1 + tmp23_1 + tmp3_1 + tmp8_1;
9145                                        EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp10_1 + tmp12_1 + tmp16_1 + tmp23_1 + tmp2_1 + tmp3_1 + tmp6_1 + tmp7_1;
9146                                    }
9147                                }
9148                            }
9149                            ///////////////
9150                            // process B //
9151                            ///////////////
9152                            if (!B.isEmpty()) {
9153                                add_EM_S=true;
9154                                const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
9155                                for (index_t k=0; k<numEq; k++) {
9156                                    for (index_t m=0; m<numComp; m++) {
9157                                        const register double B_0 = B_p[INDEX3(k,0,m, numEq, 3)];
9158                                        const register double B_1 = B_p[INDEX3(k,1,m, numEq, 3)];
9159                                        const register double B_2 = B_p[INDEX3(k,2,m, numEq, 3)];
9160                                        const register double tmp4_1 = B_0*w15;
9161                                        const register double tmp3_1 = B_1*w16;
9162                                        const register double tmp2_1 = B_0*w12;
9163                                        const register double tmp5_1 = B_2*w17;
9164                                        const register double tmp1_1 = B_2*w14;
9165                                        const register double tmp0_1 = B_1*w13;
9166                                        EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9167                                        EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
9168                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
9169                                        EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
9170                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
9171                                        EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
9172                                        EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9173                                        EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
9174                                        EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
9175                                        EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
9176                                        EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
9177                                        EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
9178                                        EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
9179                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9180                                        EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9181                                        EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
9182                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
9183                                        EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
9184                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
9185                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
9186                                        EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9187                                        EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
9188                                        EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
9189                                        EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
9190                                        EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
9191                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
9192                                        EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9193                                        EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
9194                                        EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
9195                                        EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
9196                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
9197                                        EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
9198                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
9199                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9200                                        EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
9201                                        EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
9202                                        EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9203                                        EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
9204                                        EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9205                                        EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
9206                                        EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
9207                                        EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
9208                                        EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9209                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
9210                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
9211                                        EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
9212                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
9213                                        EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
9214                                        EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
9215                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9216                                        EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9217                                        EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
9218                                        EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp3_1;
9219                                        EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
9220                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
9221                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9222                                        EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9223                                        EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9224                                        EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp4_1;
9225                                        EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=tmp1_1 + tmp3_1 + tmp4_1;
9226                                        EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp5_1;
9227                                        EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9228                                        EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp5_1;
9229                                        EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp0_1 + tmp4_1 + tmp5_1;
9230                                    }
9231                                }
9232                            }
9233                            ///////////////
9234                            // process C //
9235                            ///////////////
9236                            if (!C.isEmpty()) {
9237                                add_EM_S=true;
9238                                const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
9239                                for (index_t k=0; k<numEq; k++) {
9240                                    for (index_t m=0; m<numComp; m++) {
9241                                        const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
9242                                        const register double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
9243                                        const register double C_2 = C_p[INDEX3(k, m, 2, numEq, numComp)];
9244                                        const register double tmp5_1 = C_0*w15;
9245                                        const register double tmp2_1 = C_0*w12;
9246                                        const register double tmp4_1 = C_1*w16;
9247                                        const register double tmp1_1 = C_2*w17;
9248                                        const register double tmp3_1 = C_2*w14;
9249                                        const register double tmp0_1 = C_1*w13;
9250                                        EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9251                                        EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
9252                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9253                                        EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9254                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
9255                                        EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9256                                        EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
9257                                        EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
9258                                        EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
9259                                        EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
9260                                        EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
9261                                        EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
9262                                        EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
9263                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9264                                        EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
9265                                        EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
9266                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
9267                                        EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
9268                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9269                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
9270                                        EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
9271                                        EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9272                                        EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9273                                        EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
9274                                        EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
9275                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
9276                                        EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
9277                                        EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
9278                                        EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9279                                        EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
9280                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
9281                                        EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
9282                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
9283                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
9284                                        EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
9285                                        EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
9286                                        EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
9287                                        EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
9288                                        EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9289                                        EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
9290                                        EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9291                                        EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
9292                                        EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
9293                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9294                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
9295                                        EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
9296                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
9297                                        EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
9298                                        EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
9299                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
9300                                        EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=tmp1_1 + tmp4_1 + tmp5_1;
9301                                        EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
9302                                        EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
9303                                        EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
9304                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp1_1 + tmp2_1 + tmp4_1;
9305                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
9306                                        EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9307                                        EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp0_1 + tmp3_1 + tmp5_1;
9308                                        EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp5_1;
9309                                        EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=tmp0_1 + tmp1_1 + tmp2_1;
9310                                        EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp0_1 + tmp2_1 + tmp3_1;
9311                                        EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp2_1 + tmp3_1 + tmp4_1;
9312                                        EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9313                                        EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp3_1 + tmp4_1 + tmp5_1;
9314                                    }
9315                                }
9316                            }
9317                            ///////////////
9318                            // process D //
9319                            ///////////////
9320                            if (!D.isEmpty()) {
9321                                add_EM_S=true;
9322                                const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
9323                                for (index_t k=0; k<numEq; k++) {
9324                                    for (index_t m=0; m<numComp; m++) {
9325                                        const register double D_0 = D_p[INDEX2(k, m, numEq)];
9326                                        const register double tmp0_1 = D_0*w18;
9327                                        EM_S[INDEX4(k,m,7,3,numEq,numComp,8)]+=tmp0_1;
9328                                        EM_S[INDEX4(k,m,4,7,numEq,numComp,8)]+=tmp0_1;
9329                                        EM_S[INDEX4(k,m,1,3,numEq,numComp,8)]+=tmp0_1;
9330                                        EM_S[INDEX4(k,m,6,4,numEq,numComp,8)]+=tmp0_1;
9331                                        EM_S[INDEX4(k,m,3,0,numEq,numComp,8)]+=tmp0_1;
9332                                        EM_S[INDEX4(k,m,5,4,numEq,numComp,8)]+=tmp0_1;
9333                                        EM_S[INDEX4(k,m,0,7,numEq,numComp,8)]+=tmp0_1;
9334                                        EM_S[INDEX4(k,m,5,6,numEq,numComp,8)]+=tmp0_1;
9335                                        EM_S[INDEX4(k,m,2,6,numEq,numComp,8)]+=tmp0_1;
9336                                        EM_S[INDEX4(k,m,1,6,numEq,numComp,8)]+=tmp0_1;
9337                                        EM_S[INDEX4(k,m,5,1,numEq,numComp,8)]+=tmp0_1;
9338                                        EM_S[INDEX4(k,m,3,7,numEq,numComp,8)]+=tmp0_1;
9339                                        EM_S[INDEX4(k,m,2,5,numEq,numComp,8)]+=tmp0_1;
9340                                        EM_S[INDEX4(k,m,0,3,numEq,numComp,8)]+=tmp0_1;
9341                                        EM_S[INDEX4(k,m,7,2,numEq,numComp,8)]+=tmp0_1;
9342                                        EM_S[INDEX4(k,m,4,0,numEq,numComp,8)]+=tmp0_1;
9343                                        EM_S[INDEX4(k,m,1,2,numEq,numComp,8)]+=tmp0_1;
9344                                        EM_S[INDEX4(k,m,6,7,numEq,numComp,8)]+=tmp0_1;
9345                                        EM_S[INDEX4(k,m,3,3,numEq,numComp,8)]+=tmp0_1;
9346                                        EM_S[INDEX4(k,m,2,0,numEq,numComp,8)]+=tmp0_1;
9347                                        EM_S[INDEX4(k,m,7,6,numEq,numComp,8)]+=tmp0_1;
9348                                        EM_S[INDEX4(k,m,4,4,numEq,numComp,8)]+=tmp0_1;
9349                                        EM_S[INDEX4(k,m,6,3,numEq,numComp,8)]+=tmp0_1;
9350                                        EM_S[INDEX4(k,m,1,5,numEq,numComp,8)]+=tmp0_1;
9351                                        EM_S[INDEX4(k,m,3,6,numEq,numComp,8)]+=tmp0_1;
9352                                        EM_S[INDEX4(k,m,2,2,numEq,numComp,8)]+=tmp0_1;
9353                                        EM_S[INDEX4(k,m,7,7,numEq,numComp,8)]+=tmp0_1;
9354                                        EM_S[INDEX4(k,m,5,7,numEq,numComp,8)]+=tmp0_1;
9355                                        EM_S[INDEX4(k,m,5,3,numEq,numComp,8)]+=tmp0_1;
9356                                        EM_S[INDEX4(k,m,4,1,numEq,numComp,8)]+=tmp0_1;
9357                                        EM_S[INDEX4(k,m,1,1,numEq,numComp,8)]+=tmp0_1;
9358                                        EM_S[INDEX4(k,m,2,7,numEq,numComp,8)]+=tmp0_1;
9359                                        EM_S[INDEX4(k,m,3,2,numEq,numComp,8)]+=tmp0_1;
9360                                        EM_S[INDEX4(k,m,0,0,numEq,numComp,8)]+=tmp0_1;
9361                                        EM_S[INDEX4(k,m,6,6,numEq,numComp,8)]+=tmp0_1;
9362                                        EM_S[INDEX4(k,m,5,0,numEq,numComp,8)]+=tmp0_1;
9363                                        EM_S[INDEX4(k,m,7,1,numEq,numComp,8)]+=tmp0_1;
9364                                        EM_S[INDEX4(k,m,4,5,numEq,numComp,8)]+=tmp0_1;
9365                                        EM_S[INDEX4(k,m,0,4,numEq,numComp,8)]+=tmp0_1;
9366                                        EM_S[INDEX4(k,m,5,5,numEq,numComp,8)]+=tmp0_1;
9367                                        EM_S[INDEX4(k,m,1,4,numEq,numComp,8)]+=tmp0_1;
9368                                        EM_S[INDEX4(k,m,6,0,numEq,numComp,8)]+=tmp0_1;
9369                                        EM_S[INDEX4(k,m,7,5,numEq,numComp,8)]+=tmp0_1;
9370                                        EM_S[INDEX4(k,m,2,3,numEq,numComp,8)]+=tmp0_1;
9371                                        EM_S[INDEX4(k,m,2,1,numEq,numComp,8)]+=tmp0_1;
9372                                        EM_S[INDEX4(k,m,4,2,numEq,numComp,8)]+=tmp0_1;
9373                                        EM_S[INDEX4(k,m,1,0,numEq,numComp,8)]+=tmp0_1;
9374                                        EM_S[INDEX4(k,m,6,5,numEq,numComp,8)]+=tmp0_1;
9375                                        EM_S[INDEX4(k,m,3,5,numEq,numComp,8)]+=tmp0_1;
9376                                        EM_S[INDEX4(k,m,0,1,numEq,numComp,8)]+=tmp0_1;
9377                                        EM_S[INDEX4(k,m,7,0,numEq,numComp,8)]+=tmp0_1;
9378                                        EM_S[INDEX4(k,m,4,6,numEq,numComp,8)]+=tmp0_1;
9379                                        EM_S[INDEX4(k,m,5,2,numEq,numComp,8)]+=tmp0_1;
9380                                        EM_S[INDEX4(k,m,6,1,numEq,numComp,8)]+=tmp0_1;
9381                                        EM_S[INDEX4(k,m,3,1,numEq,numComp,8)]+=tmp0_1;
9382                                        EM_S[INDEX4(k,m,0,2,numEq,numComp,8)]+=tmp0_1;
9383                                        EM_S[INDEX4(k,m,7,4,numEq,numComp,8)]+=tmp0_1;
9384                                        EM_S[INDEX4(k,m,0,6,numEq,numComp,8)]+=tmp0_1;
9385                                        EM_S[INDEX4(k,m,6,2,numEq,numComp,8)]+=tmp0_1;
9386                                        EM_S[INDEX4(k,m,4,3,numEq,numComp,8)]+=tmp0_1;
9387                                        EM_S[INDEX4(k,m,1,7,numEq,numComp,8)]+=tmp0_1;
9388                                        EM_S[INDEX4(k,m,0,5,numEq,numComp,8)]+=tmp0_1;
9389                                        EM_S[INDEX4(k,m,3,4,numEq,numComp,8)]+=tmp0_1;
9390                                        EM_S[INDEX4(k,m,2,4,numEq,numComp,8)]+=tmp0_1;
9391                                    }
9392                                }
9393                            }
9394                            ///////////////
9395                            // process X //
9396                            ///////////////
9397                            if (!X.isEmpty()) {
9398                                add_EM_F=true;
9399                                const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
9400                                for (index_t k=0; k<numEq; k++) {
9401                                    const register double X_0 = X_p[INDEX2(k, 0, numEq)];
9402                                    const register double X_1 = X_p[INDEX2(k, 1, numEq)];
9403                                    const register double X_2 = X_p[INDEX2(k, 2, numEq)];
9404                                    const register double tmp1_1 = X_0*w19;
9405                                    const register double tmp2_1 = X_1*w20;
9406                                    const register double tmp3_1 = X_0*w22;
9407                                    const register double tmp4_1 = X_1*w23;
9408                                    const register double tmp5_1 = X_2*w24;
9409                                    const register double tmp0_1 = X_2*w21;
9410                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
9411                                    EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1 + tmp3_1;
9412                                    EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp1_1 + tmp4_1;
9413                                    EM_F[INDEX2(k,3,numEq)]+=tmp0_1 + tmp3_1 + tmp4_1;
9414                                    EM_F[INDEX2(k,4,numEq)]+=tmp1_1 + tmp2_1 + tmp5_1;
9415                                    EM_F[INDEX2(k,5,numEq)]+=tmp2_1 + tmp3_1 + tmp5_1;
9416                                    EM_F[INDEX2(k,6,numEq)]+=tmp1_1 + tmp4_1 + tmp5_1;
9417                                    EM_F[INDEX2(k,7,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;
9418                                }
9419                            }
9420                            ///////////////
9421                            // process Y //
9422                            ///////////////
9423                            if (!Y.isEmpty()) {
9424                                add_EM_F=true;
9425                                const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
9426                                for (index_t k=0; k<numEq; k++) {
9427                                    const register double Y_0 = Y_p[k];
9428                                    const register double tmp0_1 = Y_0*w25;
9429                                    EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
9430                                    EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
9431                                    EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
9432                                    EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
9433                                    EM_F[INDEX2(k,4,numEq)]+=tmp0_1;
9434                                    EM_F[INDEX2(k,5,numEq)]+=tmp0_1;
9435                                    EM_F[INDEX2(k,6,numEq)]+=tmp0_1;
9436                                    EM_F[INDEX2(k,7,numEq)]+=tmp0_1;
9437                                }
9438                            }
9439                            /* GENERATOR SNIP_PDE_SYSTEM_REDUCED BOTTOM */
9440    
9441                            // add to matrix (if add_EM_S) and RHS (if add_EM_F)
9442                            const index_t firstNode=m_N0*m_N1*k2+m_N0*k1+k0;
9443                          IndexVector rowIndex;                          IndexVector rowIndex;
9444                          rowIndex.push_back(m_dofMap[firstNode]);                          rowIndex.push_back(m_dofMap[firstNode]);
9445                          rowIndex.push_back(m_dofMap[firstNode+1]);                          rowIndex.push_back(m_dofMap[firstNode+1]);

Legend:
Removed from v.3761  
changed lines
  Added in v.3762

  ViewVC Help
Powered by ViewVC 1.1.26