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

Diff of /trunk/ripley/src/Rectangle.cpp

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

revision 3769 by caltinay, Tue Jan 17 04:09:55 2012 UTC revision 3770 by gross, Wed Jan 18 01:40:15 2012 UTC
# Line 3385  void Rectangle::assemblePDESystemReduced Line 3385  void Rectangle::assemblePDESystemReduced
3385    
3386  //protected  //protected
3387  void Rectangle::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
3388        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
3389            const escript::Data& C, const escript::Data& D,
3390            const escript::Data& X, const escript::Data& Y) const
3391  {  {
3392        const double h0 = m_l0/m_gNE0;
3393        const double h1 = m_l1/m_gNE1;
3394        dim_t numEq, numComp;
3395        if (!mat)
3396            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
3397        else {
3398            numEq=mat->logical_row_block_size;
3399            numComp=mat->logical_col_block_size;
3400        }
3401      /* GENERATOR SNIP_PDEBC_SINGLE_PRE TOP */
3402      const double w0 = -0.31100423396407310779*h1/pow(h0,2);
3403      const double w1 = -0.39433756729740644113/h0;
3404      const double w10 = 0.083333333333333333333*h1/pow(h0,2);
3405      const double w11 = -0.5/h1;
3406      const double w12 = -0.33333333333333333333*h1/pow(h0,2);
3407      const double w13 = -0.5/h0;
3408      const double w14 = -0.16666666666666666667*h1/pow(h0,2);
3409      const double w15 = 0.5/h0;
3410      const double w16 = 0.33333333333333333333*h1/pow(h0,2);
3411      const double w17 = 1.0/h1;
3412      const double w18 = 0.16666666666666666667*h1/pow(h0,2);
3413      const double w19 = -1.0/h1;
3414      const double w2 = -0.022329099369260225539*h1/pow(h0,2);
3415      const double w20 = 0.083333333333333333333*h1/h0;
3416      const double w21 = 0.022329099369260225539*h1/h0;
3417      const double w22 = 0.31100423396407310779*h1/h0;
3418      const double w23 = -0.31100423396407310779*h1/h0;
3419      const double w24 = -0.39433756729740643311;
3420      const double w25 = -0.022329099369260225539*h1/h0;
3421      const double w26 = -0.10566243270259356689;
3422      const double w27 = -0.083333333333333333333*h1/h0;
3423      const double w28 = 0.39433756729740643311;
3424      const double w29 = 0.10566243270259356689;
3425      const double w3 = -0.10566243270259355887/h0;
3426      const double w30 = 0.16666666666666666667*h1/h0;
3427      const double w31 = 0.33333333333333333333*h1/h0;
3428      const double w32 = -0.33333333333333333333*h1/h0;
3429      const double w33 = -0.50000000000000000000;
3430      const double w34 = -0.16666666666666666667*h1/h0;
3431      const double w35 = 0.50000000000000000000;
3432      const double w36 = 0.31100423396407310779*h1;
3433      const double w37 = 0.022329099369260225539*h1;
3434      const double w38 = 0.083333333333333333333*h1;
3435      const double w39 = 0.33333333333333333333*h1;
3436      const double w4 = -0.083333333333333333333*h1/pow(h0,2);
3437      const double w40 = 0.16666666666666666667*h1;
3438      const double w41 = -0.39433756729740644113*h1/h0;
3439      const double w42 = -0.10566243270259355887*h1/h0;
3440      const double w43 = 0.39433756729740644113*h1/h0;
3441      const double w44 = 0.10566243270259355887*h1/h0;
3442      const double w45 = -0.5*h1/h0;
3443      const double w46 = -1.0000000000000000000;
3444      const double w47 = 0.5*h1/h0;
3445      const double w48 = 1.0000000000000000000;
3446      const double w49 = 0.39433756729740644113*h1;
3447      const double w5 = 0.39433756729740644113/h0;
3448      const double w50 = 0.10566243270259355887*h1;
3449      const double w51 = 0.5*h1;
3450      const double w52 = -0.5/h0;
3451      const double w53 = 0.10566243270259355887/h1;
3452      const double w54 = -0.39433756729740644113/h1;
3453      const double w55 = 0.083333333333333333333*h0/pow(h1,2);
3454      const double w56 = 0.39433756729740644113/h1;
3455      const double w57 = -0.10566243270259355887/h1;
3456      const double w58 = -0.083333333333333333333*h0/pow(h1,2);
3457      const double w59 = 0.5/h0;
3458      const double w6 = 0.10566243270259355887/h0;
3459      const double w60 = 0.31100423396407310779*h0/pow(h1,2);
3460      const double w61 = 0.022329099369260225539*h0/pow(h1,2);
3461      const double w62 = -0.022329099369260225539*h0/pow(h1,2);
3462      const double w63 = -0.31100423396407310779*h0/pow(h1,2);
3463      const double w64 = -1.0/h0;
3464      const double w65 = 0.5/h1;
3465      const double w66 = -0.5/h1;
3466      const double w67 = 0.16666666666666666667*h0/pow(h1,2);
3467      const double w68 = -0.16666666666666666667*h0/pow(h1,2);
3468      const double w69 = 1.0/h0;
3469      const double w7 = 0.31100423396407310779*h1/pow(h0,2);
3470      const double w70 = 0.33333333333333333333*h0/pow(h1,2);
3471      const double w71 = -0.33333333333333333333*h0/pow(h1,2);
3472      const double w72 = -0.083333333333333333333*h0/h1;
3473      const double w73 = -0.31100423396407310779*h0/h1;
3474      const double w74 = -0.022329099369260225539*h0/h1;
3475      const double w75 = 0.083333333333333333333*h0/h1;
3476      const double w76 = 0.022329099369260225539*h0/h1;
3477      const double w77 = 0.31100423396407310779*h0/h1;
3478      const double w78 = -0.16666666666666666667*h0/h1;
3479      const double w79 = -0.33333333333333333333*h0/h1;
3480      const double w8 = 0.5/h1;
3481      const double w80 = 0.16666666666666666667*h0/h1;
3482      const double w81 = 0.33333333333333333333*h0/h1;
3483      const double w82 = 0.083333333333333333333*h0;
3484      const double w83 = 0.31100423396407310779*h0;
3485      const double w84 = 0.022329099369260225539*h0;
3486      const double w85 = 0.16666666666666666667*h0;
3487      const double w86 = 0.33333333333333333333*h0;
3488      const double w87 = -0.39433756729740644113*h0/h1;
3489      const double w88 = -0.10566243270259355887*h0/h1;
3490      const double w89 = 0.39433756729740644113*h0/h1;
3491      const double w9 = 0.022329099369260225539*h1/pow(h0,2);
3492      const double w90 = 0.10566243270259355887*h0/h1;
3493      const double w91 = -0.5*h0/h1;
3494      const double w92 = 0.5*h0/h1;
3495      const double w93 = 0.39433756729740644113*h0;
3496      const double w94 = 0.10566243270259355887*h0;
3497      const double w95 = 0.5*h0;
3498      /* GENERATOR SNIP_PDEBC_SINGLE_PRE BOTTOM */
3499      #pragma omp parallel
3500            {
3501                if (m_faceOffset[0] > -1) {
3502                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3503    #pragma omp for nowait
3504                        for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3505                  /* GENERATOR SNIP_PDEBC_SINGLE_0 TOP */
3506            ///////////////
3507            // process A //
3508            ///////////////
3509            if (!A.isEmpty()) {
3510                add_EM_S=true;
3511                const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
3512                if (A.actsExpanded()) {
3513                    const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
3514                    const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
3515                    const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
3516                    const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
3517                    const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
3518                    const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
3519                    const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
3520                    const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
3521                    const register double tmp1_0 = A_11_0 + A_11_1;
3522                    const register double tmp2_0 = A_01_1 + A_10_1;
3523                    const register double tmp0_0 = A_00_0 + A_00_1;
3524                    const register double tmp3_0 = A_01_0 + A_10_0;
3525                    const register double tmp7_1 = A_00_1*w0;
3526                    const register double tmp16_1 = A_00_0*w9;
3527                    const register double tmp23_1 = A_01_0*w1;
3528                    const register double tmp4_1 = A_01_1*w6;
3529                    const register double tmp13_1 = tmp2_0*w6;
3530                    const register double tmp30_1 = A_10_0*w3;
3531                    const register double tmp2_1 = A_10_0*w1;
3532                    const register double tmp14_1 = A_00_0*w7;
3533                    const register double tmp5_1 = tmp0_0*w4;
3534                    const register double tmp15_1 = tmp3_0*w5;
3535                    const register double tmp22_1 = A_10_0*w5;
3536                    const register double tmp10_1 = A_00_0*w2;
3537                    const register double tmp20_1 = tmp0_0*w10;
3538                    const register double tmp29_1 = tmp2_0*w1;
3539                    const register double tmp0_1 = A_10_1*w3;
3540                    const register double tmp31_1 = A_10_1*w1;
3541                    const register double tmp11_1 = tmp1_0*w8;
3542                    const register double tmp21_1 = A_10_1*w6;
3543                    const register double tmp9_1 = A_01_1*w5;
3544                    const register double tmp12_1 = A_00_1*w9;
3545                    const register double tmp27_1 = A_10_1*w5;
3546                    const register double tmp6_1 = A_01_0*w5;
3547                    const register double tmp8_1 = A_01_0*w6;
3548                    const register double tmp1_1 = A_00_0*w0;
3549                    const register double tmp3_1 = A_00_1*w2;
3550                    const register double tmp19_1 = A_01_1*w1;
3551                    const register double tmp25_1 = A_01_1*w3;
3552                    const register double tmp18_1 = A_01_0*w3;
3553                    const register double tmp24_1 = tmp1_0*w11;
3554                    const register double tmp26_1 = A_10_0*w6;
3555                    const register double tmp17_1 = A_00_1*w7;
3556                    const register double tmp28_1 = tmp3_0*w3;
3557                    EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
3558                    EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
3559                    EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp7_1 + tmp8_1 + tmp9_1;
3560                    EM_S[INDEX2(0,0,4)]+=tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
3561                    EM_S[INDEX2(3,3,4)]+=tmp16_1 + tmp17_1;
3562                    EM_S[INDEX2(3,0,4)]+=tmp18_1 + tmp19_1 + tmp5_1;
3563                    EM_S[INDEX2(3,1,4)]+=tmp20_1;
3564                    EM_S[INDEX2(2,1,4)]+=tmp21_1 + tmp22_1 + tmp5_1;
3565                    EM_S[INDEX2(0,2,4)]+=tmp20_1 + tmp23_1 + tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1;
3566                    EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp20_1 + tmp24_1 + tmp2_1 + tmp8_1 + tmp9_1;
3567                    EM_S[INDEX2(1,3,4)]+=tmp20_1;
3568                    EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp26_1 + tmp27_1 + tmp7_1;
3569                    EM_S[INDEX2(2,2,4)]+=tmp11_1 + tmp16_1 + tmp17_1 + tmp28_1 + tmp29_1;
3570                    EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp23_1 + tmp25_1 + tmp3_1;
3571                    EM_S[INDEX2(0,3,4)]+=tmp30_1 + tmp31_1 + tmp5_1;
3572                    EM_S[INDEX2(1,1,4)]+=tmp12_1 + tmp14_1;
3573                } else { /* constant data */
3574                    const register double A_00 = A_p[INDEX2(0,0,2)];
3575                    const register double A_01 = A_p[INDEX2(0,1,2)];
3576                    const register double A_10 = A_p[INDEX2(1,0,2)];
3577                    const register double A_11 = A_p[INDEX2(1,1,2)];
3578                    const register double tmp0_0 = A_01 + A_10;
3579                    const register double tmp5_1 = A_11*w17;
3580                    const register double tmp8_1 = A_00*w18;
3581                    const register double tmp0_1 = A_10*w13;
3582                    const register double tmp11_1 = tmp0_0*w13;
3583                    const register double tmp3_1 = A_01*w15;
3584                    const register double tmp10_1 = A_11*w19;
3585                    const register double tmp4_1 = A_00*w16;
3586                    const register double tmp2_1 = A_00*w14;
3587                    const register double tmp6_1 = tmp0_0*w15;
3588                    const register double tmp1_1 = A_00*w12;
3589                    const register double tmp7_1 = A_01*w13;
3590                    const register double tmp9_1 = A_10*w15;
3591                    EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
3592                    EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
3593                    EM_S[INDEX2(3,2,4)]+=tmp1_1 + tmp3_1;
3594                    EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
3595                    EM_S[INDEX2(3,3,4)]+=tmp4_1;
3596                    EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp7_1;
3597                    EM_S[INDEX2(3,1,4)]+=tmp8_1;
3598                    EM_S[INDEX2(2,1,4)]+=tmp2_1 + tmp9_1;
3599                    EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp7_1 + tmp8_1 + tmp9_1;
3600                    EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp10_1 + tmp3_1 + tmp8_1;
3601                    EM_S[INDEX2(1,3,4)]+=tmp8_1;
3602                    EM_S[INDEX2(2,3,4)]+=tmp1_1 + tmp9_1;
3603                    EM_S[INDEX2(2,2,4)]+=tmp11_1 + tmp4_1 + tmp5_1;
3604                    EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp7_1;
3605                    EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;
3606                    EM_S[INDEX2(1,1,4)]+=tmp4_1;
3607                }
3608            }
3609            ///////////////
3610            // process B //
3611            ///////////////
3612            if (!B.isEmpty()) {
3613                add_EM_S=true;
3614                const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
3615                if (B.actsExpanded()) {
3616                    const register double B_0_0 = B_p[INDEX2(0,0,2)];
3617                    const register double B_1_0 = B_p[INDEX2(1,0,2)];
3618                    const register double B_0_1 = B_p[INDEX2(0,1,2)];
3619                    const register double B_1_1 = B_p[INDEX2(1,1,2)];
3620                    const register double tmp0_0 = B_0_0 + B_0_1;
3621                    const register double tmp7_1 = tmp0_0*w27;
3622                    const register double tmp5_1 = B_1_0*w24;
3623                    const register double tmp8_1 = B_1_0*w26;
3624                    const register double tmp11_1 = B_1_0*w28;
3625                    const register double tmp10_1 = B_1_1*w29;
3626                    const register double tmp15_1 = B_0_1*w23;
3627                    const register double tmp0_1 = tmp0_0*w20;
3628                    const register double tmp4_1 = B_0_1*w25;
3629                    const register double tmp9_1 = B_1_1*w24;
3630                    const register double tmp16_1 = B_0_0*w22;
3631                    const register double tmp3_1 = B_1_1*w26;
3632                    const register double tmp13_1 = B_1_1*w28;
3633                    const register double tmp12_1 = B_1_0*w29;
3634                    const register double tmp1_1 = B_0_1*w22;
3635                    const register double tmp17_1 = B_0_1*w21;
3636                    const register double tmp6_1 = B_0_0*w23;
3637                    const register double tmp2_1 = B_0_0*w21;
3638                    const register double tmp14_1 = B_0_0*w25;
3639                    EM_S[INDEX2(1,2,4)]+=tmp0_1;
3640                    EM_S[INDEX2(3,2,4)]+=tmp1_1 + tmp2_1;
3641                    EM_S[INDEX2(0,0,4)]+=tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1;
3642                    EM_S[INDEX2(3,0,4)]+=tmp0_1;
3643                    EM_S[INDEX2(0,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
3644                    EM_S[INDEX2(2,0,4)]+=tmp10_1 + tmp11_1 + tmp7_1;
3645                    EM_S[INDEX2(2,2,4)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
3646                    EM_S[INDEX2(1,0,4)]+=tmp16_1 + tmp17_1;
3647                } else { /* constant data */
3648                    const register double B_0 = B_p[0];
3649                    const register double B_1 = B_p[1];
3650                    const register double tmp5_1 = B_1*w35;
3651                    const register double tmp1_1 = B_0*w31;
3652                    const register double tmp3_1 = B_1*w33;
3653                    const register double tmp0_1 = B_0*w30;
3654                    const register double tmp2_1 = B_0*w32;
3655                    const register double tmp4_1 = B_0*w34;
3656                    EM_S[INDEX2(1,2,4)]+=tmp0_1;
3657                    EM_S[INDEX2(3,2,4)]+=tmp1_1;
3658                    EM_S[INDEX2(0,0,4)]+=tmp2_1 + tmp3_1;
3659                    EM_S[INDEX2(3,0,4)]+=tmp0_1;
3660                    EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp4_1;
3661                    EM_S[INDEX2(2,0,4)]+=tmp4_1 + tmp5_1;
3662                    EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp5_1;
3663                    EM_S[INDEX2(1,0,4)]+=tmp1_1;
3664                }
3665            }
3666            ///////////////
3667            // process C //
3668            ///////////////
3669            if (!C.isEmpty()) {
3670                add_EM_S=true;
3671                const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
3672                if (C.actsExpanded()) {
3673                    const register double C_0_0 = C_p[INDEX2(0,0,2)];
3674                    const register double C_1_0 = C_p[INDEX2(1,0,2)];
3675                    const register double C_0_1 = C_p[INDEX2(0,1,2)];
3676                    const register double C_1_1 = C_p[INDEX2(1,1,2)];
3677                    const register double tmp0_0 = C_0_0 + C_0_1;
3678                    const register double tmp14_1 = C_1_0*w29;
3679                    const register double tmp9_1 = tmp0_0*w27;
3680                    const register double tmp11_1 = C_1_1*w24;
3681                    const register double tmp2_1 = C_1_1*w26;
3682                    const register double tmp15_1 = C_1_1*w28;
3683                    const register double tmp12_1 = C_0_1*w22;
3684                    const register double tmp0_1 = C_0_0*w22;
3685                    const register double tmp6_1 = tmp0_0*w20;
3686                    const register double tmp17_1 = C_0_1*w23;
3687                    const register double tmp3_1 = C_0_1*w25;
3688                    const register double tmp16_1 = C_0_0*w25;
3689                    const register double tmp4_1 = C_1_0*w24;
3690                    const register double tmp10_1 = C_1_0*w26;
3691                    const register double tmp8_1 = C_1_0*w28;
3692                    const register double tmp7_1 = C_1_1*w29;
3693                    const register double tmp5_1 = C_0_0*w23;
3694                    const register double tmp1_1 = C_0_1*w21;
3695                    const register double tmp13_1 = C_0_0*w21;
3696                    EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
3697                    EM_S[INDEX2(0,0,4)]+=tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
3698                    EM_S[INDEX2(2,1,4)]+=tmp6_1;
3699                    EM_S[INDEX2(0,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
3700                    EM_S[INDEX2(2,0,4)]+=tmp10_1 + tmp11_1 + tmp9_1;
3701                    EM_S[INDEX2(2,3,4)]+=tmp12_1 + tmp13_1;
3702                    EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1;
3703                    EM_S[INDEX2(0,3,4)]+=tmp6_1;
3704                } else { /* constant data */
3705                    const register double C_0 = C_p[0];
3706                    const register double C_1 = C_p[1];
3707                    const register double tmp3_1 = C_0*w30;
3708                    const register double tmp1_1 = C_0*w32;
3709                    const register double tmp5_1 = C_0*w34;
3710                    const register double tmp2_1 = C_1*w33;
3711                    const register double tmp4_1 = C_1*w35;
3712                    const register double tmp0_1 = C_0*w31;
3713                    EM_S[INDEX2(0,1,4)]+=tmp0_1;
3714                    EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp2_1;
3715                    EM_S[INDEX2(2,1,4)]+=tmp3_1;
3716                    EM_S[INDEX2(0,2,4)]+=tmp4_1 + tmp5_1;
3717                    EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp5_1;
3718                    EM_S[INDEX2(2,3,4)]+=tmp0_1;
3719                    EM_S[INDEX2(2,2,4)]+=tmp1_1 + tmp4_1;
3720                    EM_S[INDEX2(0,3,4)]+=tmp3_1;
3721                }
3722            }
3723            ///////////////
3724            // process D //
3725            ///////////////
3726            if (!D.isEmpty()) {
3727                add_EM_S=true;
3728                const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
3729                if (D.actsExpanded()) {
3730                    const register double D_0 = D_p[0];
3731                    const register double D_1 = D_p[1];
3732                    const register double tmp0_0 = D_0 + D_1;
3733                    const register double tmp1_1 = D_1*w37;
3734                    const register double tmp4_1 = D_0*w37;
3735                    const register double tmp0_1 = D_0*w36;
3736                    const register double tmp3_1 = D_1*w36;
3737                    const register double tmp2_1 = tmp0_0*w38;
3738                    EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
3739                    EM_S[INDEX2(0,2,4)]+=tmp2_1;
3740                    EM_S[INDEX2(2,0,4)]+=tmp2_1;
3741                    EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;
3742                } else { /* constant data */
3743                    const register double D_0 = D_p[0];
3744                    const register double tmp0_1 = D_0*w39;
3745                    const register double tmp1_1 = D_0*w40;
3746                    EM_S[INDEX2(0,0,4)]+=tmp0_1;
3747                    EM_S[INDEX2(0,2,4)]+=tmp1_1;
3748                    EM_S[INDEX2(2,0,4)]+=tmp1_1;
3749                    EM_S[INDEX2(2,2,4)]+=tmp0_1;
3750                }
3751            }
3752            ///////////////
3753            // process X //
3754            ///////////////
3755            if (!X.isEmpty()) {
3756                add_EM_F=true;
3757                const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
3758                if (X.actsExpanded()) {
3759                    const register double X_0_0 = X_p[INDEX2(0,0,2)];
3760                    const register double X_1_0 = X_p[INDEX2(1,0,2)];
3761                    const register double X_0_1 = X_p[INDEX2(0,1,2)];
3762                    const register double X_1_1 = X_p[INDEX2(1,1,2)];
3763                    const register double tmp0_0 = X_1_0 + X_1_1;
3764                    const register double tmp3_1 = X_0_1*w44;
3765                    const register double tmp1_1 = X_0_1*w42;
3766                    const register double tmp8_1 = X_0_0*w44;
3767                    const register double tmp4_1 = X_0_0*w43;
3768                    const register double tmp7_1 = X_0_0*w42;
3769                    const register double tmp5_1 = tmp0_0*w35;
3770                    const register double tmp0_1 = X_0_0*w41;
3771                    const register double tmp9_1 = X_0_1*w43;
3772                    const register double tmp2_1 = tmp0_0*w33;
3773                    const register double tmp6_1 = X_0_1*w41;
3774                    EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
3775                    EM_F[1]+=tmp3_1 + tmp4_1;
3776                    EM_F[2]+=tmp5_1 + tmp6_1 + tmp7_1;
3777                    EM_F[3]+=tmp8_1 + tmp9_1;
3778                } else { /* constant data */
3779                    const register double X_0 = X_p[0];
3780                    const register double X_1 = X_p[1];
3781                    const register double tmp2_1 = X_0*w47;
3782                    const register double tmp1_1 = X_0*w45;
3783                    const register double tmp3_1 = X_1*w48;
3784                    const register double tmp0_1 = X_1*w46;
3785                    EM_F[0]+=tmp0_1 + tmp1_1;
3786                    EM_F[1]+=tmp2_1;
3787                    EM_F[2]+=tmp1_1 + tmp3_1;
3788                    EM_F[3]+=tmp2_1;
3789                }
3790            }
3791            ///////////////
3792            // process Y //
3793            ///////////////
3794            if (!Y.isEmpty()) {
3795                add_EM_F=true;
3796                const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
3797                if (Y.actsExpanded()) {
3798                    const register double Y_0 = Y_p[0];
3799                    const register double Y_1 = Y_p[1];
3800                    const register double tmp0_1 = Y_0*w49;
3801                    const register double tmp1_1 = Y_1*w50;
3802                    const register double tmp3_1 = Y_0*w50;
3803                    const register double tmp2_1 = Y_1*w49;
3804                    EM_F[0]+=tmp0_1 + tmp1_1;
3805                    EM_F[2]+=tmp2_1 + tmp3_1;
3806                } else { /* constant data */
3807                    const register double Y_0 = Y_p[0];
3808                    const register double tmp0_1 = Y_0*w51;
3809                    EM_F[0]+=tmp0_1;
3810                    EM_F[2]+=tmp0_1;
3811                }
3812            }
3813                          /* GENERATOR SNIP_PDEBC_SINGLE_0 BOTTOM */
3814                        }
3815                        // ADD EM_F  and EM_S
3816            } // end colouring
3817                }
3818    
3819                if (m_faceOffset[1] > -1) {
3820                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3821    #pragma omp for nowait
3822                        for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3823                  /* GENERATOR SNIP_PDEBC_SINGLE_1 TOP */
3824            ///////////////
3825            // process A //
3826            ///////////////
3827            if (!A.isEmpty()) {
3828                add_EM_S=true;
3829                const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
3830                if (A.actsExpanded()) {
3831                    const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
3832                    const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
3833                    const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
3834                    const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
3835                    const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
3836                    const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
3837                    const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
3838                    const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
3839                    const register double tmp1_0 = A_11_0 + A_11_1;
3840                    const register double tmp3_0 = A_01_1 + A_10_1;
3841                    const register double tmp0_0 = A_00_0 + A_00_1;
3842                    const register double tmp2_0 = A_01_0 + A_10_0;
3843                    const register double tmp8_1 = A_00_1*w0;
3844                    const register double tmp14_1 = A_00_0*w9;
3845                    const register double tmp29_1 = A_01_0*w1;
3846                    const register double tmp1_1 = A_01_1*w6;
3847                    const register double tmp15_1 = tmp2_0*w6;
3848                    const register double tmp7_1 = A_10_0*w3;
3849                    const register double tmp19_1 = A_10_0*w1;
3850                    const register double tmp12_1 = A_00_0*w7;
3851                    const register double tmp5_1 = tmp0_0*w4;
3852                    const register double tmp17_1 = tmp3_0*w5;
3853                    const register double tmp25_1 = A_10_0*w5;
3854                    const register double tmp10_1 = A_00_0*w2;
3855                    const register double tmp23_1 = tmp0_0*w10;
3856                    const register double tmp31_1 = tmp2_0*w1;
3857                    const register double tmp18_1 = A_10_1*w3;
3858                    const register double tmp9_1 = A_10_1*w1;
3859                    const register double tmp13_1 = tmp1_0*w8;
3860                    const register double tmp24_1 = A_10_1*w6;
3861                    const register double tmp27_1 = A_01_1*w5;
3862                    const register double tmp11_1 = A_00_1*w9;
3863                    const register double tmp6_1 = A_10_1*w5;
3864                    const register double tmp2_1 = A_01_0*w5;
3865                    const register double tmp26_1 = A_01_0*w6;
3866                    const register double tmp0_1 = A_00_0*w0;
3867                    const register double tmp3_1 = A_00_1*w2;
3868                    const register double tmp20_1 = A_01_1*w1;
3869                    const register double tmp28_1 = A_01_1*w3;
3870                    const register double tmp22_1 = A_01_0*w3;
3871                    const register double tmp21_1 = tmp1_0*w11;
3872                    const register double tmp4_1 = A_10_0*w6;
3873                    const register double tmp16_1 = A_00_1*w7;
3874                    const register double tmp30_1 = tmp3_0*w3;
3875                    EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
3876                    EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
3877                    EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp7_1 + tmp8_1 + tmp9_1;
3878                    EM_S[INDEX2(0,0,4)]+=tmp11_1 + tmp12_1;
3879                    EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1;
3880                    EM_S[INDEX2(3,0,4)]+=tmp18_1 + tmp19_1 + tmp5_1;
3881                    EM_S[INDEX2(3,1,4)]+=tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
3882                    EM_S[INDEX2(2,1,4)]+=tmp26_1 + tmp27_1 + tmp5_1;
3883                    EM_S[INDEX2(0,2,4)]+=tmp23_1;
3884                    EM_S[INDEX2(2,0,4)]+=tmp23_1;
3885                    EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp21_1 + tmp23_1 + tmp2_1 + tmp7_1 + tmp9_1;
3886                    EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp20_1 + tmp22_1 + tmp8_1;
3887                    EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp16_1;
3888                    EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp24_1 + tmp25_1 + tmp3_1;
3889                    EM_S[INDEX2(0,3,4)]+=tmp28_1 + tmp29_1 + tmp5_1;
3890                    EM_S[INDEX2(1,1,4)]+=tmp11_1 + tmp12_1 + tmp13_1 + tmp30_1 + tmp31_1;
3891                } else { /* constant data */
3892                    const register double A_00 = A_p[INDEX2(0,0,2)];
3893                    const register double A_01 = A_p[INDEX2(0,1,2)];
3894                    const register double A_10 = A_p[INDEX2(1,0,2)];
3895                    const register double A_11 = A_p[INDEX2(1,1,2)];
3896                    const register double tmp0_0 = A_01 + A_10;
3897                    const register double tmp6_1 = A_11*w17;
3898                    const register double tmp9_1 = A_00*w18;
3899                    const register double tmp4_1 = A_10*w13;
3900                    const register double tmp11_1 = tmp0_0*w13;
3901                    const register double tmp0_1 = A_01*w15;
3902                    const register double tmp10_1 = A_11*w19;
3903                    const register double tmp5_1 = A_00*w16;
3904                    const register double tmp2_1 = A_00*w14;
3905                    const register double tmp7_1 = tmp0_0*w15;
3906                    const register double tmp1_1 = A_00*w12;
3907                    const register double tmp8_1 = A_01*w13;
3908                    const register double tmp3_1 = A_10*w15;
3909                    EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
3910                    EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
3911                    EM_S[INDEX2(3,2,4)]+=tmp1_1 + tmp4_1;
3912                    EM_S[INDEX2(0,0,4)]+=tmp5_1;
3913                    EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
3914                    EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;
3915                    EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp3_1 + tmp8_1 + tmp9_1;
3916                    EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp2_1;
3917                    EM_S[INDEX2(0,2,4)]+=tmp9_1;
3918                    EM_S[INDEX2(2,0,4)]+=tmp9_1;
3919                    EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp10_1 + tmp4_1 + tmp9_1;
3920                    EM_S[INDEX2(2,3,4)]+=tmp1_1 + tmp8_1;
3921                    EM_S[INDEX2(2,2,4)]+=tmp5_1;
3922                    EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
3923                    EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp8_1;
3924                    EM_S[INDEX2(1,1,4)]+=tmp11_1 + tmp5_1 + tmp6_1;
3925                }
3926            }
3927            ///////////////
3928            // process B //
3929            ///////////////
3930            if (!B.isEmpty()) {
3931                add_EM_S=true;
3932                const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
3933                if (B.actsExpanded()) {
3934                    const register double B_0_0 = B_p[INDEX2(0,0,2)];
3935                    const register double B_1_0 = B_p[INDEX2(1,0,2)];
3936                    const register double B_0_1 = B_p[INDEX2(0,1,2)];
3937                    const register double B_1_1 = B_p[INDEX2(1,1,2)];
3938                    const register double tmp0_0 = B_0_0 + B_0_1;
3939                    const register double tmp9_1 = tmp0_0*w27;
3940                    const register double tmp17_1 = B_1_0*w24;
3941                    const register double tmp10_1 = B_1_0*w26;
3942                    const register double tmp8_1 = B_1_0*w28;
3943                    const register double tmp6_1 = B_1_1*w29;
3944                    const register double tmp13_1 = B_0_1*w23;
3945                    const register double tmp7_1 = tmp0_0*w20;
3946                    const register double tmp0_1 = B_0_1*w25;
3947                    const register double tmp11_1 = B_1_1*w24;
3948                    const register double tmp14_1 = B_0_0*w22;
3949                    const register double tmp16_1 = B_1_1*w26;
3950                    const register double tmp4_1 = B_1_1*w28;
3951                    const register double tmp5_1 = B_1_0*w29;
3952                    const register double tmp3_1 = B_0_0*w21;
3953                    const register double tmp15_1 = B_0_1*w21;
3954                    const register double tmp1_1 = B_0_0*w23;
3955                    const register double tmp2_1 = B_0_1*w22;
3956                    const register double tmp12_1 = B_0_0*w25;
3957                    EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
3958                    EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
3959                    EM_S[INDEX2(3,1,4)]+=tmp6_1 + tmp7_1 + tmp8_1;
3960                    EM_S[INDEX2(2,1,4)]+=tmp9_1;
3961                    EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1 + tmp7_1;
3962                    EM_S[INDEX2(2,3,4)]+=tmp12_1 + tmp13_1;
3963                    EM_S[INDEX2(0,3,4)]+=tmp9_1;
3964                    EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1;
3965                } else { /* constant data */
3966                    const register double B_0 = B_p[0];
3967                    const register double B_1 = B_p[1];
3968                    const register double tmp1_1 = B_1*w35;
3969                    const register double tmp2_1 = B_0*w31;
3970                    const register double tmp5_1 = B_1*w33;
3971                    const register double tmp3_1 = B_0*w30;
3972                    const register double tmp0_1 = B_0*w32;
3973                    const register double tmp4_1 = B_0*w34;
3974                    EM_S[INDEX2(0,1,4)]+=tmp0_1;
3975                    EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;
3976                    EM_S[INDEX2(3,1,4)]+=tmp1_1 + tmp3_1;
3977                    EM_S[INDEX2(2,1,4)]+=tmp4_1;
3978                    EM_S[INDEX2(1,3,4)]+=tmp3_1 + tmp5_1;
3979                    EM_S[INDEX2(2,3,4)]+=tmp0_1;
3980                    EM_S[INDEX2(0,3,4)]+=tmp4_1;
3981                    EM_S[INDEX2(1,1,4)]+=tmp2_1 + tmp5_1;
3982                }
3983            }
3984            ///////////////
3985            // process C //
3986            ///////////////
3987            if (!C.isEmpty()) {
3988                add_EM_S=true;
3989                const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
3990                if (C.actsExpanded()) {
3991                    const register double C_0_0 = C_p[INDEX2(0,0,2)];
3992                    const register double C_1_0 = C_p[INDEX2(1,0,2)];
3993                    const register double C_0_1 = C_p[INDEX2(0,1,2)];
3994                    const register double C_1_1 = C_p[INDEX2(1,1,2)];
3995                    const register double tmp0_0 = C_0_0 + C_0_1;
3996                    const register double tmp7_1 = tmp0_0*w20;
3997                    const register double tmp15_1 = C_0_1*w21;
3998                    const register double tmp0_1 = tmp0_0*w27;
3999                    const register double tmp9_1 = C_1_1*w24;
4000                    const register double tmp16_1 = C_1_1*w26;
4001                    const register double tmp5_1 = C_1_1*w28;
4002                    const register double tmp3_1 = C_0_1*w22;
4003                    const register double tmp6_1 = C_1_0*w29;
4004                    const register double tmp11_1 = C_1_0*w28;
4005                    const register double tmp2_1 = C_0_1*w23;
4006                    const register double tmp12_1 = C_0_1*w25;
4007                    const register double tmp17_1 = C_1_0*w24;
4008                    const register double tmp8_1 = C_1_0*w26;
4009                    const register double tmp4_1 = C_0_0*w21;
4010                    const register double tmp10_1 = C_1_1*w29;
4011                    const register double tmp13_1 = C_0_0*w23;
4012                    const register double tmp1_1 = C_0_0*w25;
4013                    const register double tmp14_1 = C_0_0*w22;
4014                    EM_S[INDEX2(1,2,4)]+=tmp0_1;
4015                    EM_S[INDEX2(3,2,4)]+=tmp1_1 + tmp2_1;
4016                    EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1;
4017                    EM_S[INDEX2(3,0,4)]+=tmp0_1;
4018                    EM_S[INDEX2(3,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
4019                    EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1 + tmp7_1;
4020                    EM_S[INDEX2(1,0,4)]+=tmp12_1 + tmp13_1;
4021                    EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1;
4022                } else { /* constant data */
4023                    const register double C_0 = C_p[0];
4024                    const register double C_1 = C_p[1];
4025                    const register double tmp4_1 = C_0*w30;
4026                    const register double tmp1_1 = C_0*w32;
4027                    const register double tmp0_1 = C_0*w34;
4028                    const register double tmp5_1 = C_1*w33;
4029                    const register double tmp2_1 = C_1*w35;
4030                    const register double tmp3_1 = C_0*w31;
4031                    EM_S[INDEX2(1,2,4)]+=tmp0_1;
4032                    EM_S[INDEX2(3,2,4)]+=tmp1_1;
4033                    EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;
4034                    EM_S[INDEX2(3,0,4)]+=tmp0_1;
4035                    EM_S[INDEX2(3,1,4)]+=tmp4_1 + tmp5_1;
4036                    EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp4_1;
4037                    EM_S[INDEX2(1,0,4)]+=tmp1_1;
4038                    EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;
4039                }
4040            }
4041            ///////////////
4042            // process D //
4043            ///////////////
4044            if (!D.isEmpty()) {
4045                add_EM_S=true;
4046                const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
4047                if (D.actsExpanded()) {
4048                    const register double D_0 = D_p[0];
4049                    const register double D_1 = D_p[1];
4050                    const register double tmp0_0 = D_0 + D_1;
4051                    const register double tmp4_1 = D_1*w37;
4052                    const register double tmp2_1 = tmp0_0*w38;
4053                    const register double tmp1_1 = D_0*w37;
4054                    const register double tmp3_1 = D_0*w36;
4055                    const register double tmp0_1 = D_1*w36;
4056                    EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp1_1;
4057                    EM_S[INDEX2(3,1,4)]+=tmp2_1;
4058                    EM_S[INDEX2(1,3,4)]+=tmp2_1;
4059                    EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;
4060                } else { /* constant data */
4061                    const register double D_0 = D_p[0];
4062                    const register double tmp0_1 = D_0*w39;
4063                    const register double tmp1_1 = D_0*w40;
4064                    EM_S[INDEX2(3,3,4)]+=tmp0_1;
4065                    EM_S[INDEX2(3,1,4)]+=tmp1_1;
4066                    EM_S[INDEX2(1,3,4)]+=tmp1_1;
4067                    EM_S[INDEX2(1,1,4)]+=tmp0_1;
4068                }
4069            }
4070            ///////////////
4071            // process X //
4072            ///////////////
4073            if (!X.isEmpty()) {
4074                add_EM_F=true;
4075                const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
4076                if (X.actsExpanded()) {
4077                    const register double X_0_0 = X_p[INDEX2(0,0,2)];
4078                    const register double X_1_0 = X_p[INDEX2(1,0,2)];
4079                    const register double X_0_1 = X_p[INDEX2(0,1,2)];
4080                    const register double X_1_1 = X_p[INDEX2(1,1,2)];
4081                    const register double tmp0_0 = X_1_0 + X_1_1;
4082                    const register double tmp2_1 = X_0_1*w44;
4083                    const register double tmp1_1 = X_0_1*w42;
4084                    const register double tmp7_1 = X_0_0*w44;
4085                    const register double tmp3_1 = X_0_0*w43;
4086                    const register double tmp6_1 = X_0_0*w42;
4087                    const register double tmp8_1 = tmp0_0*w35;
4088                    const register double tmp0_1 = X_0_0*w41;
4089                    const register double tmp9_1 = X_0_1*w43;
4090                    const register double tmp4_1 = tmp0_0*w33;
4091                    const register double tmp5_1 = X_0_1*w41;
4092                    EM_F[0]+=tmp0_1 + tmp1_1;
4093                    EM_F[1]+=tmp2_1 + tmp3_1 + tmp4_1;
4094                    EM_F[2]+=tmp5_1 + tmp6_1;
4095                    EM_F[3]+=tmp7_1 + tmp8_1 + tmp9_1;
4096                } else { /* constant data */
4097                    const register double X_0 = X_p[0];
4098                    const register double X_1 = X_p[1];
4099                    const register double tmp1_1 = X_0*w47;
4100                    const register double tmp0_1 = X_0*w45;
4101                    const register double tmp3_1 = X_1*w48;
4102                    const register double tmp2_1 = X_1*w46;
4103                    EM_F[0]+=tmp0_1;
4104                    EM_F[1]+=tmp1_1 + tmp2_1;
4105                    EM_F[2]+=tmp0_1;
4106                    EM_F[3]+=tmp1_1 + tmp3_1;
4107                }
4108            }
4109            ///////////////
4110            // process Y //
4111            ///////////////
4112            if (!Y.isEmpty()) {
4113                add_EM_F=true;
4114                const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
4115                if (Y.actsExpanded()) {
4116                    const register double Y_0 = Y_p[0];
4117                    const register double Y_1 = Y_p[1];
4118                    const register double tmp0_1 = Y_0*w49;
4119                    const register double tmp1_1 = Y_1*w50;
4120                    const register double tmp3_1 = Y_0*w50;
4121                    const register double tmp2_1 = Y_1*w49;
4122                    EM_F[1]+=tmp0_1 + tmp1_1;
4123                    EM_F[3]+=tmp2_1 + tmp3_1;
4124                } else { /* constant data */
4125                    const register double Y_0 = Y_p[0];
4126                    const register double tmp0_1 = Y_0*w51;
4127                    EM_F[1]+=tmp0_1;
4128                    EM_F[3]+=tmp0_1;
4129                }
4130            }
4131                          /* GENERATOR SNIP_PDEBC_SINGLE_1 BOTTOM */
4132                        }
4133                         // ADD EM_F  and EM_S
4134                    } // end colouring
4135                }
4136    
4137                if (m_faceOffset[2] > -1) {
4138               for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4139    #pragma omp for nowait
4140                      for (index_t k0 = 0; k0 < m_NE0; ++k0) {
4141                          /* GENERATOR SNIP_PDEBC_SINGLE_2 TOP */
4142                          ///////////////
4143                          // process A //
4144                          ///////////////
4145                          if (!A.isEmpty()) {
4146                              add_EM_S=true;
4147                              const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
4148                              if (A.actsExpanded()) {
4149                                  const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
4150                                  const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
4151                                  const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
4152                                  const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
4153                                  const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
4154                                  const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
4155                                  const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
4156                                  const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
4157                                  const register double tmp1_0 = A_11_0 + A_11_1;
4158                                  const register double tmp2_0 = A_01_1 + A_10_1;
4159                                  const register double tmp0_0 = A_00_0 + A_00_1;
4160                                  const register double tmp3_0 = A_01_0 + A_10_0;
4161                                  const register double tmp27_1 = A_11_1*w62;
4162                                  const register double tmp23_1 = A_10_0*w56;
4163                                  const register double tmp15_1 = A_11_1*w60;
4164                                  const register double tmp4_1 = A_10_0*w54;
4165                                  const register double tmp12_1 = tmp0_0*w59;
4166                                  const register double tmp25_1 = A_01_0*w54;
4167                                  const register double tmp21_1 = A_10_1*w56;
4168                                  const register double tmp10_1 = tmp2_0*w53;
4169                                  const register double tmp16_1 = A_10_0*w57;
4170                                  const register double tmp18_1 = A_11_0*w62;
4171                                  const register double tmp11_1 = A_11_0*w60;
4172                                  const register double tmp8_1 = A_01_0*w56;
4173                                  const register double tmp6_1 = A_01_1*w53;
4174                                  const register double tmp20_1 = A_10_0*w53;
4175                                  const register double tmp0_1 = A_10_1*w57;
4176                                  const register double tmp19_1 = A_11_1*w63;
4177                                  const register double tmp5_1 = A_01_0*w53;
4178                                  const register double tmp9_1 = A_11_1*w61;
4179                                  const register double tmp1_1 = tmp0_0*w52;
4180                                  const register double tmp24_1 = A_01_1*w57;
4181                                  const register double tmp7_1 = tmp1_0*w58;
4182                                  const register double tmp17_1 = A_10_1*w54;
4183                                  const register double tmp29_1 = A_01_1*w54;
4184                                  const register double tmp3_1 = tmp1_0*w55;
4185                                  const register double tmp26_1 = A_11_0*w63;
4186                                  const register double tmp22_1 = A_10_1*w53;
4187                                  const register double tmp14_1 = A_11_0*w61;
4188                                  const register double tmp31_1 = tmp2_0*w54;
4189                                  const register double tmp30_1 = tmp3_0*w57;
4190                                  const register double tmp28_1 = A_01_0*w57;
4191                                  const register double tmp2_1 = A_01_1*w56;
4192                                  const register double tmp13_1 = tmp3_0*w56;
4193                                  EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
4194                                  EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1;
4195                                  EM_S[INDEX2(3,2,4)]+=tmp3_1;
4196                                  EM_S[INDEX2(0,0,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;
4197                                  EM_S[INDEX2(3,3,4)]+=tmp14_1 + tmp15_1;
4198                                  EM_S[INDEX2(3,0,4)]+=tmp16_1 + tmp17_1 + tmp7_1;
4199                                  EM_S[INDEX2(3,1,4)]+=tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1;
4200                                  EM_S[INDEX2(2,1,4)]+=tmp22_1 + tmp23_1 + tmp7_1;
4201                                  EM_S[INDEX2(0,2,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1;
4202                                  EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp26_1 + tmp27_1 + tmp4_1;
4203                                  EM_S[INDEX2(1,3,4)]+=tmp18_1 + tmp19_1 + tmp2_1 + tmp5_1;
4204                                  EM_S[INDEX2(2,3,4)]+=tmp3_1;
4205                                  EM_S[INDEX2(2,2,4)]+=tmp11_1 + tmp9_1;
4206                                  EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp20_1 + tmp21_1 + tmp24_1 + tmp25_1 + tmp3_1;
4207                                  EM_S[INDEX2(0,3,4)]+=tmp28_1 + tmp29_1 + tmp7_1;
4208                                  EM_S[INDEX2(1,1,4)]+=tmp12_1 + tmp14_1 + tmp15_1 + tmp30_1 + tmp31_1;
4209                              } else { /* constant data */
4210                                  const register double A_00 = A_p[INDEX2(0,0,2)];
4211                                  const register double A_01 = A_p[INDEX2(0,1,2)];
4212                                  const register double A_10 = A_p[INDEX2(1,0,2)];
4213                                  const register double A_11 = A_p[INDEX2(1,1,2)];
4214                                  const register double tmp0_0 = A_01 + A_10;
4215                                  const register double tmp0_1 = A_10*w66;
4216                                  const register double tmp11_1 = tmp0_0*w66;
4217                                  const register double tmp6_1 = A_00*w69;
4218                                  const register double tmp9_1 = A_10*w65;
4219                                  const register double tmp7_1 = A_11*w70;
4220                                  const register double tmp8_1 = A_11*w71;
4221                                  const register double tmp3_1 = A_11*w67;
4222                                  const register double tmp5_1 = tmp0_0*w65;
4223                                  const register double tmp2_1 = A_00*w64;
4224                                  const register double tmp4_1 = A_11*w68;
4225                                  const register double tmp10_1 = A_01*w66;
4226                                  const register double tmp1_1 = A_01*w65;
4227                                  EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
4228                                  EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1;
4229                                  EM_S[INDEX2(3,2,4)]+=tmp3_1;
4230                                  EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
4231                                  EM_S[INDEX2(3,3,4)]+=tmp7_1;
4232                                  EM_S[INDEX2(3,0,4)]+=tmp0_1 + tmp4_1;
4233                                  EM_S[INDEX2(3,1,4)]+=tmp8_1 + tmp9_1;
4234                                  EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp9_1;
4235                                  EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp8_1;
4236                                  EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp8_1;
4237                                  EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp8_1;
4238                                  EM_S[INDEX2(2,3,4)]+=tmp3_1;
4239                                  EM_S[INDEX2(2,2,4)]+=tmp7_1;
4240                                  EM_S[INDEX2(1,0,4)]+=tmp10_1 + tmp2_1 + tmp3_1 + tmp9_1;
4241                                  EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp4_1;
4242                                  EM_S[INDEX2(1,1,4)]+=tmp11_1 + tmp6_1 + tmp7_1;
4243                              }
4244                          }
4245                          ///////////////
4246                          // process B //
4247                          ///////////////
4248                          if (!B.isEmpty()) {
4249                              add_EM_S=true;
4250                              const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
4251                              if (B.actsExpanded()) {
4252                                  const register double B_0_0 = B_p[INDEX2(0,0,2)];
4253                                  const register double B_1_0 = B_p[INDEX2(1,0,2)];
4254                                  const register double B_0_1 = B_p[INDEX2(0,1,2)];
4255                                  const register double B_1_1 = B_p[INDEX2(1,1,2)];
4256                                  const register double tmp0_0 = B_1_0 + B_1_1;
4257                                  const register double tmp15_1 = B_0_1*w28;
4258                                  const register double tmp14_1 = B_0_0*w29;
4259                                  const register double tmp9_1 = B_1_1*w77;
4260                                  const register double tmp8_1 = B_1_0*w76;
4261                                  const register double tmp16_1 = B_1_0*w74;
4262                                  const register double tmp7_1 = tmp0_0*w75;
4263                                  const register double tmp2_1 = B_0_1*w24;
4264                                  const register double tmp0_1 = tmp0_0*w72;
4265                                  const register double tmp5_1 = B_0_1*w26;
4266                                  const register double tmp17_1 = B_1_1*w73;
4267                                  const register double tmp3_1 = B_1_0*w73;
4268                                  const register double tmp13_1 = B_0_1*w29;
4269                                  const register double tmp11_1 = B_1_0*w77;
4270                                  const register double tmp10_1 = B_1_1*w76;
4271                                  const register double tmp4_1 = B_1_1*w74;
4272                                  const register double tmp6_1 = B_0_0*w24;
4273                                  const register double tmp1_1 = B_0_0*w26;
4274                                  const register double tmp12_1 = B_0_0*w28;
4275                                  EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1;
4276                                  EM_S[INDEX2(0,0,4)]+=tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1;
4277                                  EM_S[INDEX2(3,0,4)]+=tmp7_1;
4278                                  EM_S[INDEX2(3,1,4)]+=tmp8_1 + tmp9_1;
4279                                  EM_S[INDEX2(2,1,4)]+=tmp7_1;
4280                                  EM_S[INDEX2(2,0,4)]+=tmp10_1 + tmp11_1;
4281                                  EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp12_1 + tmp13_1;
4282                                  EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1;
4283                              } else { /* constant data */
4284                                  const register double B_0 = B_p[0];
4285                                  const register double B_1 = B_p[1];
4286                                  const register double tmp4_1 = B_1*w81;
4287                                  const register double tmp0_1 = B_0*w33;
4288                                  const register double tmp3_1 = B_1*w80;
4289                                  const register double tmp2_1 = B_1*w79;
4290                                  const register double tmp5_1 = B_0*w35;
4291                                  const register double tmp1_1 = B_1*w78;
4292                                  EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
4293                                  EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp2_1;
4294                                  EM_S[INDEX2(3,0,4)]+=tmp3_1;
4295                                  EM_S[INDEX2(3,1,4)]+=tmp4_1;
4296                                  EM_S[INDEX2(2,1,4)]+=tmp3_1;
4297                                  EM_S[INDEX2(2,0,4)]+=tmp4_1;
4298                                  EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp5_1;
4299                                  EM_S[INDEX2(1,1,4)]+=tmp2_1 + tmp5_1;
4300                              }
4301                          }
4302                          ///////////////
4303                          // process C //
4304                          ///////////////
4305                          if (!C.isEmpty()) {
4306                              add_EM_S=true;
4307                              const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
4308                              if (C.actsExpanded()) {
4309                                  const register double C_0_0 = C_p[INDEX2(0,0,2)];
4310                                  const register double C_1_0 = C_p[INDEX2(1,0,2)];
4311                                  const register double C_0_1 = C_p[INDEX2(0,1,2)];
4312                                  const register double C_1_1 = C_p[INDEX2(1,1,2)];
4313                                  const register double tmp0_0 = C_1_0 + C_1_1;
4314                                  const register double tmp1_1 = C_0_0*w28;
4315                                  const register double tmp14_1 = C_0_0*w29;
4316                                  const register double tmp9_1 = C_1_0*w77;
4317                                  const register double tmp8_1 = C_1_1*w76;
4318                                  const register double tmp5_1 = C_1_1*w74;
4319                                  const register double tmp3_1 = tmp0_0*w75;
4320                                  const register double tmp13_1 = C_0_1*w24;
4321                                  const register double tmp0_1 = tmp0_0*w72;
4322                                  const register double tmp4_1 = C_1_0*w73;
4323                                  const register double tmp17_1 = C_1_1*w73;
4324                                  const register double tmp2_1 = C_0_1*w29;
4325                                  const register double tmp11_1 = C_1_1*w77;
4326                                  const register double tmp10_1 = C_1_0*w76;
4327                                  const register double tmp16_1 = C_1_0*w74;
4328                                  const register double tmp6_1 = C_0_1*w26;
4329                                  const register double tmp15_1 = C_0_1*w28;
4330                                  const register double tmp7_1 = C_0_0*w24;
4331                                  const register double tmp12_1 = C_0_0*w26;
4332                                  EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1;
4333                                  EM_S[INDEX2(1,2,4)]+=tmp3_1;
4334                                  EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
4335                                  EM_S[INDEX2(0,2,4)]+=tmp8_1 + tmp9_1;
4336                                  EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
4337                                  EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp12_1 + tmp13_1;
4338                                  EM_S[INDEX2(0,3,4)]+=tmp3_1;
4339                                  EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1;
4340                              } else { /* constant data */
4341                                  const register double C_0 = C_p[0];
4342                                  const register double C_1 = C_p[1];
4343                                  const register double tmp5_1 = C_1*w81;
4344                                  const register double tmp2_1 = C_1*w80;
4345                                  const register double tmp4_1 = C_0*w33;
4346                                  const register double tmp1_1 = C_1*w78;
4347                                  const register double tmp0_1 = C_0*w35;
4348                                  const register double tmp3_1 = C_1*w79;
4349                                  EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
4350                                  EM_S[INDEX2(1,2,4)]+=tmp2_1;
4351                                  EM_S[INDEX2(0,0,4)]+=tmp3_1 + tmp4_1;
4352                                  EM_S[INDEX2(0,2,4)]+=tmp5_1;
4353                                  EM_S[INDEX2(1,3,4)]+=tmp5_1;
4354                                  EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;
4355                                  EM_S[INDEX2(0,3,4)]+=tmp2_1;
4356                                  EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp3_1;
4357                              }
4358                          }
4359                          ///////////////
4360                          // process D //
4361                          ///////////////
4362                          if (!D.isEmpty()) {
4363                              add_EM_S=true;
4364                              const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
4365                              if (D.actsExpanded()) {
4366                                  const register double D_0 = D_p[0];
4367                                  const register double D_1 = D_p[1];
4368                                  const register double tmp0_0 = D_0 + D_1;
4369                                  const register double tmp2_1 = D_1*w84;
4370                                  const register double tmp4_1 = D_0*w84;
4371                                  const register double tmp0_1 = tmp0_0*w82;
4372                                  const register double tmp1_1 = D_0*w83;
4373                                  const register double tmp3_1 = D_1*w83;
4374                                  EM_S[INDEX2(0,1,4)]+=tmp0_1;
4375                                  EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp2_1;
4376                                  EM_S[INDEX2(1,0,4)]+=tmp0_1;
4377                                  EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;
4378                              } else { /* constant data */
4379                                  const register double D_0 = D_p[0];
4380                                  const register double tmp1_1 = D_0*w86;
4381                                  const register double tmp0_1 = D_0*w85;
4382                                  EM_S[INDEX2(0,1,4)]+=tmp0_1;
4383                                  EM_S[INDEX2(0,0,4)]+=tmp1_1;
4384                                  EM_S[INDEX2(1,0,4)]+=tmp0_1;
4385                                  EM_S[INDEX2(1,1,4)]+=tmp1_1;
4386                              }
4387                          }
4388                          ///////////////
4389                          // process X //
4390                          ///////////////
4391                          if (!X.isEmpty()) {
4392                              add_EM_F=true;
4393                              const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
4394                              if (X.actsExpanded()) {
4395                                  const register double X_0_0 = X_p[INDEX2(0,0,2)];
4396                                  const register double X_1_0 = X_p[INDEX2(1,0,2)];
4397                                  const register double X_0_1 = X_p[INDEX2(0,1,2)];
4398                                  const register double X_1_1 = X_p[INDEX2(1,1,2)];
4399                                  const register double tmp0_0 = X_0_0 + X_0_1;
4400                                  const register double tmp4_1 = X_1_0*w88;
4401                                  const register double tmp6_1 = X_1_1*w90;
4402                                  const register double tmp9_1 = X_1_1*w89;
4403                                  const register double tmp0_1 = X_1_1*w88;
4404                                  const register double tmp8_1 = X_1_0*w90;
4405                                  const register double tmp3_1 = tmp0_0*w35;
4406                                  const register double tmp5_1 = X_1_1*w87;
4407                                  const register double tmp7_1 = X_1_0*w89;
4408                                  const register double tmp2_1 = tmp0_0*w33;
4409                                  const register double tmp1_1 = X_1_0*w87;
4410                                  EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;
4411                                  EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;
4412                                  EM_F[2]+=tmp6_1 + tmp7_1;
4413                                  EM_F[3]+=tmp8_1 + tmp9_1;
4414                              } else { /* constant data */
4415                                  const register double X_0 = X_p[0];
4416                                  const register double X_1 = X_p[1];
4417                                  const register double tmp0_1 = X_0*w46;
4418                                  const register double tmp3_1 = X_1*w92;
4419                                  const register double tmp2_1 = X_0*w48;
4420                                  const register double tmp1_1 = X_1*w91;
4421                                  EM_F[0]+=tmp0_1 + tmp1_1;
4422                                  EM_F[1]+=tmp1_1 + tmp2_1;
4423                                  EM_F[2]+=tmp3_1;
4424                                  EM_F[3]+=tmp3_1;
4425                              }
4426                          }
4427                          ///////////////
4428                          // process Y //
4429                          ///////////////
4430                          if (!Y.isEmpty()) {
4431                              add_EM_F=true;
4432                              const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
4433                              if (Y.actsExpanded()) {
4434                                  const register double Y_0 = Y_p[0];
4435                                  const register double Y_1 = Y_p[1];
4436                                  const register double tmp0_1 = Y_0*w93;
4437                                  const register double tmp2_1 = Y_1*w93;
4438                                  const register double tmp1_1 = Y_1*w94;
4439                                  const register double tmp3_1 = Y_0*w94;
4440                                  EM_F[0]+=tmp0_1 + tmp1_1;
4441                                  EM_F[1]+=tmp2_1 + tmp3_1;
4442                              } else { /* constant data */
4443                                  const register double Y_0 = Y_p[0];
4444                                  const register double tmp0_1 = Y_0*w95;
4445                                  EM_F[0]+=tmp0_1;
4446                                  EM_F[1]+=tmp0_1;
4447                              }
4448                          }
4449                          /* GENERATOR SNIP_PDEBC_SINGLE_2 BOTTOM */
4450                      }
4451                      // ADD EM_F  and EM_S
4452                   } // end colouring
4453                }
4454    
4455                if (m_faceOffset[3] > -1) {
4456                 for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4457    #pragma omp for nowait
4458                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
4459                           /* GENERATOR SNIP_PDEBC_SINGLE_3 TOP */
4460                           ///////////////
4461                           // process A //
4462                           ///////////////
4463                           if (!A.isEmpty()) {
4464                               add_EM_S=true;
4465                               const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
4466                               if (A.actsExpanded()) {
4467                                   const register double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
4468                                   const register double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
4469                                   const register double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
4470                                   const register double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
4471                                   const register double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
4472                                   const register double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
4473                                   const register double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
4474                                   const register double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
4475                                   const register double tmp0_0 = A_11_0 + A_11_1;
4476                                   const register double tmp3_0 = A_01_1 + A_10_1;
4477                                   const register double tmp1_0 = A_00_0 + A_00_1;
4478                                   const register double tmp2_0 = A_01_0 + A_10_0;
4479                                   const register double tmp24_1 = A_11_1*w62;
4480                                   const register double tmp27_1 = A_10_0*w56;
4481                                   const register double tmp11_1 = A_11_1*w60;
4482                                   const register double tmp31_1 = A_10_0*w54;
4483                                   const register double tmp16_1 = A_01_1*w57;
4484                                   const register double tmp3_1 = A_10_1*w56;
4485                                   const register double tmp12_1 = tmp2_0*w53;
4486                                   const register double tmp4_1 = A_10_0*w57;
4487                                   const register double tmp21_1 = A_11_0*w62;
4488                                   const register double tmp10_1 = A_11_0*w60;
4489                                   const register double tmp6_1 = A_01_0*w56;
4490                                   const register double tmp5_1 = tmp1_0*w52;
4491                                   const register double tmp14_1 = tmp1_0*w59;
4492                                   const register double tmp8_1 = A_01_1*w53;
4493                                   const register double tmp1_1 = A_10_0*w53;
4494                                   const register double tmp2_1 = tmp0_0*w58;
4495                                   const register double tmp30_1 = A_10_1*w57;
4496                                   const register double tmp20_1 = A_11_1*w63;
4497                                   const register double tmp22_1 = A_01_0*w53;
4498                                   const register double tmp9_1 = A_11_1*w61;
4499                                   const register double tmp17_1 = A_01_0*w54;
4500                                   const register double tmp7_1 = A_10_1*w54;
4501                                   const register double tmp19_1 = A_01_1*w54;
4502                                   const register double tmp25_1 = A_11_0*w63;
4503                                   const register double tmp26_1 = A_10_1*w53;
4504                                   const register double tmp13_1 = A_11_0*w61;
4505                                   const register double tmp23_1 = A_01_1*w56;
4506                                   const register double tmp29_1 = tmp2_0*w54;
4507                                   const register double tmp28_1 = tmp3_0*w57;
4508                                   const register double tmp18_1 = A_01_0*w57;
4509                                   const register double tmp0_1 = tmp0_0*w55;
4510                                   const register double tmp15_1 = tmp3_0*w56;
4511                                   EM_S[INDEX2(0,1,4)]+=tmp0_1;
4512                                   EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1 + tmp3_1;
4513                                   EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;
4514                                   EM_S[INDEX2(0,0,4)]+=tmp10_1 + tmp9_1;
4515                                   EM_S[INDEX2(3,3,4)]+=tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
4516                                   EM_S[INDEX2(3,0,4)]+=tmp16_1 + tmp17_1 + tmp2_1;
4517                                   EM_S[INDEX2(3,1,4)]+=tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1;
4518                                   EM_S[INDEX2(2,1,4)]+=tmp22_1 + tmp23_1 + tmp2_1;
4519                                   EM_S[INDEX2(0,2,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1;
4520                                   EM_S[INDEX2(2,0,4)]+=tmp24_1 + tmp25_1 + tmp6_1 + tmp8_1;
4521                                   EM_S[INDEX2(1,3,4)]+=tmp20_1 + tmp21_1 + tmp4_1 + tmp7_1;
4522                                   EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp18_1 + tmp19_1 + tmp26_1 + tmp27_1 + tmp5_1;
4523                                   EM_S[INDEX2(2,2,4)]+=tmp10_1 + tmp14_1 + tmp28_1 + tmp29_1 + tmp9_1;
4524                                   EM_S[INDEX2(1,0,4)]+=tmp0_1;
4525                                   EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp30_1 + tmp31_1;
4526                                   EM_S[INDEX2(1,1,4)]+=tmp11_1 + tmp13_1;
4527                               } else { /* constant data */
4528                                   const register double A_00 = A_p[INDEX2(0,0,2)];
4529                                   const register double A_01 = A_p[INDEX2(0,1,2)];
4530                                   const register double A_10 = A_p[INDEX2(1,0,2)];
4531                                   const register double A_11 = A_p[INDEX2(1,1,2)];
4532                                   const register double tmp0_0 = A_01 + A_10;
4533                                   const register double tmp3_1 = A_10*w66;
4534                                   const register double tmp11_1 = tmp0_0*w66;
4535                                   const register double tmp6_1 = A_11*w70;
4536                                   const register double tmp1_1 = A_10*w65;
4537                                   const register double tmp8_1 = A_00*w69;
4538                                   const register double tmp10_1 = A_11*w71;
4539                                   const register double tmp0_1 = A_11*w67;
4540                                   const register double tmp7_1 = tmp0_0*w65;
4541                                   const register double tmp5_1 = A_00*w64;
4542                                   const register double tmp2_1 = A_11*w68;
4543                                   const register double tmp9_1 = A_01*w66;
4544                                   const register double tmp4_1 = A_01*w65;
4545                                   EM_S[INDEX2(0,1,4)]+=tmp0_1;
4546                                   EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
4547                                   EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp3_1 + tmp4_1 + tmp5_1;
4548                                   EM_S[INDEX2(0,0,4)]+=tmp6_1;
4549                                   EM_S[INDEX2(3,3,4)]+=tmp6_1 + tmp7_1 + tmp8_1;
4550                                   EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp9_1;
4551                                   EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp9_1;
4552                                   EM_S[INDEX2(2,1,4)]+=tmp2_1 + tmp4_1;
4553                                   EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp1_1;
4554                                   EM_S[INDEX2(2,0,4)]+=tmp10_1 + tmp4_1;
4555                                   EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp3_1;
4556                                   EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp5_1 + tmp9_1;
4557                                   EM_S[INDEX2(2,2,4)]+=tmp11_1 + tmp6_1 + tmp8_1;
4558                                   EM_S[INDEX2(1,0,4)]+=tmp0_1;
4559                                   EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp3_1;
4560                                   EM_S[INDEX2(1,1,4)]+=tmp6_1;
4561                               }
4562                           }
4563                           ///////////////
4564                           // process B //
4565                           ///////////////
4566                           if (!B.isEmpty()) {
4567                               add_EM_S=true;
4568                               const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
4569                               if (B.actsExpanded()) {
4570                                   const register double B_0_0 = B_p[INDEX2(0,0,2)];
4571                                   const register double B_1_0 = B_p[INDEX2(1,0,2)];
4572                                   const register double B_0_1 = B_p[INDEX2(0,1,2)];
4573                                   const register double B_1_1 = B_p[INDEX2(1,1,2)];
4574                                   const register double tmp0_0 = B_1_0 + B_1_1;
4575                                   const register double tmp8_1 = B_1_0*w73;
4576                                   const register double tmp4_1 = B_0_0*w29;
4577                                   const register double tmp7_1 = B_1_1*w77;
4578                                   const register double tmp6_1 = B_1_0*w76;
4579                                   const register double tmp11_1 = B_1_0*w74;
4580                                   const register double tmp3_1 = tmp0_0*w75;
4581                                   const register double tmp13_1 = B_0_1*w24;
4582                                   const register double tmp0_1 = tmp0_0*w72;
4583                                   const register double tmp14_1 = B_0_1*w26;
4584                                   const register double tmp10_1 = B_1_1*w73;
4585                                   const register double tmp5_1 = B_0_1*w28;
4586                                   const register double tmp1_1 = B_0_1*w29;
4587                                   const register double tmp15_1 = B_1_0*w77;
4588                                   const register double tmp17_1 = B_1_1*w76;
4589                                   const register double tmp9_1 = B_1_1*w74;
4590                                   const register double tmp16_1 = B_0_0*w24;
4591                                   const register double tmp12_1 = B_0_0*w26;
4592                                   const register double tmp2_1 = B_0_0*w28;
4593                                   EM_S[INDEX2(1,2,4)]+=tmp0_1;
4594                                   EM_S[INDEX2(3,2,4)]+=tmp1_1 + tmp2_1 + tmp3_1;
4595                                   EM_S[INDEX2(3,3,4)]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;
4596                                   EM_S[INDEX2(0,2,4)]+=tmp8_1 + tmp9_1;
4597                                   EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;
4598                                   EM_S[INDEX2(2,3,4)]+=tmp12_1 + tmp13_1 + tmp3_1;
4599                                   EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1;
4600                                   EM_S[INDEX2(0,3,4)]+=tmp0_1;
4601                               } else { /* constant data */
4602                                   const register double B_0 = B_p[0];
4603                                   const register double B_1 = B_p[1];
4604                                   const register double tmp3_1 = B_1*w81;
4605                                   const register double tmp5_1 = B_0*w33;
4606                                   const register double tmp2_1 = B_1*w80;
4607                                   const register double tmp4_1 = B_1*w79;
4608                                   const register double tmp1_1 = B_0*w35;
4609                                   const register double tmp0_1 = B_1*w78;
4610                                   EM_S[INDEX2(1,2,4)]+=tmp0_1;
4611                                   EM_S[INDEX2(3,2,4)]+=tmp1_1 + tmp2_1;
4612                                   EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp3_1;
4613                                   EM_S[INDEX2(0,2,4)]+=tmp4_1;
4614                                   EM_S[INDEX2(1,3,4)]+=tmp4_1;
4615                                   EM_S[INDEX2(2,3,4)]+=tmp2_1 + tmp5_1;
4616                                   EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp5_1;
4617                                   EM_S[INDEX2(0,3,4)]+=tmp0_1;
4618                               }
4619                           }
4620                           ///////////////
4621                           // process C //
4622                           ///////////////
4623                           if (!C.isEmpty()) {
4624                               add_EM_S=true;
4625                               const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
4626                               if (C.actsExpanded()) {
4627                                   const register double C_0_0 = C_p[INDEX2(0,0,2)];
4628                                   const register double C_1_0 = C_p[INDEX2(1,0,2)];
4629                                   const register double C_0_1 = C_p[INDEX2(0,1,2)];
4630                                   const register double C_1_1 = C_p[INDEX2(1,1,2)];
4631                                   const register double tmp0_0 = C_1_0 + C_1_1;
4632                                   const register double tmp13_1 = C_0_0*w28;
4633                                   const register double tmp3_1 = C_0_0*w29;
4634                                   const register double tmp15_1 = C_1_0*w77;
4635                                   const register double tmp17_1 = C_1_1*w76;
4636                                   const register double tmp11_1 = C_1_1*w74;
4637                                   const register double tmp0_1 = tmp0_0*w75;
4638                                   const register double tmp2_1 = C_0_1*w24;
4639                                   const register double tmp7_1 = tmp0_0*w72;
4640                                   const register double tmp4_1 = C_0_1*w28;
4641                                   const register double tmp8_1 = C_1_1*w73;
4642                                   const register double tmp12_1 = C_0_1*w29;
4643                                   const register double tmp6_1 = C_1_1*w77;
4644                                   const register double tmp5_1 = C_1_0*w76;
4645                                   const register double tmp9_1 = C_1_0*w74;
4646                                   const register double tmp14_1 = C_0_1*w26;
4647                                   const register double tmp10_1 = C_1_0*w73;
4648                                   const register double tmp16_1 = C_0_0*w24;
4649                                   const register double tmp1_1 = C_0_0*w26;
4650                                   EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1;
4651                                   EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1;
4652                                   EM_S[INDEX2(3,0,4)]+=tmp7_1;
4653                                   EM_S[INDEX2(3,1,4)]+=tmp8_1 + tmp9_1;
4654                                   EM_S[INDEX2(2,1,4)]+=tmp7_1;
4655                                   EM_S[INDEX2(2,0,4)]+=tmp10_1 + tmp11_1;
4656                                   EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp12_1 + tmp13_1;
4657                                   EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1;
4658                               } else { /* constant data */
4659                                   const register double C_0 = C_p[0];
4660                                   const register double C_1 = C_p[1];
4661                                   const register double tmp3_1 = C_1*w81;
4662                                   const register double tmp0_1 = C_1*w80;
4663                                   const register double tmp4_1 = C_1*w78;
4664                                   const register double tmp1_1 = C_0*w33;
4665                                   const register double tmp2_1 = C_0*w35;
4666                                   const register double tmp5_1 = C_1*w79;
4667                                   EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1;
4668                                   EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;
4669                                   EM_S[INDEX2(3,0,4)]+=tmp4_1;
4670                                   EM_S[INDEX2(3,1,4)]+=tmp5_1;
4671                                   EM_S[INDEX2(2,1,4)]+=tmp4_1;
4672                                   EM_S[INDEX2(2,0,4)]+=tmp5_1;
4673                                   EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
4674                                   EM_S[INDEX2(2,2,4)]+=tmp1_1 + tmp3_1;
4675                               }
4676                           }
4677                           ///////////////
4678                           // process D //
4679                           ///////////////
4680                           if (!D.isEmpty()) {
4681                               add_EM_S=true;
4682                               const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
4683                               if (D.actsExpanded()) {
4684                                   const register double D_0 = D_p[0];
4685                                   const register double D_1 = D_p[1];
4686                                   const register double tmp0_0 = D_0 + D_1;
4687                                   const register double tmp4_1 = D_1*w84;
4688                                   const register double tmp2_1 = D_0*w84;
4689                                   const register double tmp0_1 = tmp0_0*w82;
4690                                   const register double tmp3_1 = D_0*w83;
4691                                   const register double tmp1_1 = D_1*w83;
4692                                   EM_S[INDEX2(3,2,4)]+=tmp0_1;
4693                                   EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;
4694                                   EM_S[INDEX2(2,3,4)]+=tmp0_1;
4695                                   EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;
4696                               } else { /* constant data */
4697                                   const register double D_0 = D_p[0];
4698                                   const register double tmp1_1 = D_0*w86;
4699                                   const register double tmp0_1 = D_0*w85;
4700                                   EM_S[INDEX2(3,2,4)]+=tmp0_1;
4701                                   EM_S[INDEX2(3,3,4)]+=tmp1_1;
4702                                   EM_S[INDEX2(2,3,4)]+=tmp0_1;
4703                                   EM_S[INDEX2(2,2,4)]+=tmp1_1;
4704                               }
4705                           }
4706                           ///////////////
4707                           // process X //
4708                           ///////////////
4709                           if (!X.isEmpty()) {
4710                               add_EM_F=true;
4711                               const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
4712                               if (X.actsExpanded()) {
4713                                   const register double X_0_0 = X_p[INDEX2(0,0,2)];
4714                                   const register double X_1_0 = X_p[INDEX2(1,0,2)];
4715                                   const register double X_0_1 = X_p[INDEX2(0,1,2)];
4716                                   const register double X_1_1 = X_p[INDEX2(1,1,2)];
4717                                   const register double tmp0_0 = X_0_0 + X_0_1;
4718                                   const register double tmp2_1 = X_1_0*w88;
4719                                   const register double tmp4_1 = X_1_1*w90;
4720                                   const register double tmp9_1 = X_1_1*w89;
4721                                   const register double tmp0_1 = X_1_1*w88;
4722                                   const register double tmp8_1 = X_1_0*w90;
4723                                   const register double tmp7_1 = tmp0_0*w35;
4724                                   const register double tmp3_1 = X_1_1*w87;
4725                                   const register double tmp6_1 = X_1_0*w89;
4726                                   const register double tmp5_1 = tmp0_0*w33;
4727                                   const register double tmp1_1 = X_1_0*w87;
4728                                   EM_F[0]+=tmp0_1 + tmp1_1;
4729                                   EM_F[1]+=tmp2_1 + tmp3_1;
4730                                   EM_F[2]+=tmp4_1 + tmp5_1 + tmp6_1;
4731                                   EM_F[3]+=tmp7_1 + tmp8_1 + tmp9_1;
4732                               } else { /* constant data */
4733                                   const register double X_0 = X_p[0];
4734                                   const register double X_1 = X_p[1];
4735                                   const register double tmp1_1 = X_0*w46;
4736                                   const register double tmp2_1 = X_1*w92;
4737                                   const register double tmp3_1 = X_0*w48;
4738                                   const register double tmp0_1 = X_1*w91;
4739                                   EM_F[0]+=tmp0_1;
4740                                   EM_F[1]+=tmp0_1;
4741                                   EM_F[2]+=tmp1_1 + tmp2_1;
4742                                   EM_F[3]+=tmp2_1 + tmp3_1;
4743                               }
4744                           }
4745                           ///////////////
4746                           // process Y //
4747                           ///////////////
4748                           if (!Y.isEmpty()) {
4749                               add_EM_F=true;
4750                               const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
4751                               if (Y.actsExpanded()) {
4752                                   const register double Y_0 = Y_p[0];
4753                                   const register double Y_1 = Y_p[1];
4754                                   const register double tmp0_1 = Y_0*w93;
4755                                   const register double tmp2_1 = Y_1*w93;
4756                                   const register double tmp1_1 = Y_1*w94;
4757                                   const register double tmp3_1 = Y_0*w94;
4758                                   EM_F[2]+=tmp0_1 + tmp1_1;
4759                                   EM_F[3]+=tmp2_1 + tmp3_1;
4760                               } else { /* constant data */
4761                                   const register double Y_0 = Y_p[0];
4762                                   const register double tmp0_1 = Y_0*w95;
4763                                   EM_F[2]+=tmp0_1;
4764                                   EM_F[3]+=tmp0_1;
4765                               }
4766                           }
4767                          /* GENERATOR SNIP_PDEBC_SINGLE_3 BOTTOM */
4768                    }
4769                    // ADD EM_F  and EM_S
4770                    } // end colouring
4771                }
4772            } // end of parallel section
4773            
4774  }  }
4775    
4776  //protected  //protected
4777  void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
4778        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
4779            const escript::Data& C, const escript::Data& D,
4780            const escript::Data& X, const escript::Data& Y) const
4781  {  {
4782        const double h0 = m_l0/m_gNE0;
4783        const double h1 = m_l1/m_gNE1;
4784        dim_t numEq, numComp;
4785        if (!mat)
4786            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
4787        else {
4788            numEq=mat->logical_row_block_size;
4789            numComp=mat->logical_col_block_size;
4790        }
4791        /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_PRE TOP */
4792        const double w0 = -0.25*h1/pow(h0,2);
4793        const double w1 = -0.5/h0;
4794        const double w10 = 0.25*h1;
4795        const double w11 = -0.5*h1/h0;
4796        const double w12 = -1.0000000000000000000;
4797        const double w13 = 0.5*h1/h0;
4798        const double w14 = 1.0000000000000000000;
4799        const double w15 = 0.5*h1;
4800        const double w16 = -1/h0;
4801        const double w17 = 0.5/h1;
4802        const double w18 = -0.5/h1;
4803        const double w19 = 0.25*h0/pow(h1,2);
4804        const double w2 = 0.5/h0;
4805        const double w20 = -0.25*h0/pow(h1,2);
4806        const double w21 = 1.0/h0;
4807        const double w22 = -0.25*h0/h1;
4808        const double w23 = 0.25*h0/h1;
4809        const double w24 = 0.25*h0;
4810        const double w25 = -0.5*h0/h1;
4811        const double w26 = 0.5*h0/h1;
4812        const double w27 = 0.5*h0;
4813        const double w3 = 0.25*h1/pow(h0,2);
4814        const double w4 = 1.0/h1;
4815        const double w5 = -1/h1;
4816        const double w6 = 0.25*h1/h0;
4817        const double w7 = -0.25*h1/h0;
4818        const double w8 = -0.50000000000000000000;
4819        const double w9 = 0.50000000000000000000;
4820      /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_PRE BOTTOM */
4821      #pragma omp parallel
4822            {
4823                if (m_faceOffset[0] > -1) {
4824                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
4825    #pragma omp for nowait
4826                        for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4827                  /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_0 TOP */
4828            ///////////////
4829            // process A //
4830            ///////////////
4831            if (!A.isEmpty()) {
4832                add_EM_S=true;
4833                const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
4834                const register double A_00 = A_p[INDEX2(0,0,2)];
4835                const register double A_01 = A_p[INDEX2(0,1,2)];
4836                const register double A_10 = A_p[INDEX2(1,0,2)];
4837                const register double A_11 = A_p[INDEX2(1,1,2)];
4838                const register double tmp0_0 = A_01 + A_10;
4839                const register double tmp6_1 = A_01*w1;
4840                const register double tmp5_1 = tmp0_0*w2;
4841                const register double tmp7_1 = A_10*w2;
4842                const register double tmp0_1 = A_00*w0;
4843                const register double tmp9_1 = tmp0_0*w1;
4844                const register double tmp8_1 = A_11*w5;
4845                const register double tmp2_1 = A_01*w2;
4846                const register double tmp4_1 = A_11*w4;
4847                const register double tmp1_1 = A_10*w1;
4848                const register double tmp3_1 = A_00*w3;
4849                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
4850                EM_S[INDEX2(1,2,4)]+=tmp0_1 + tmp2_1;
4851                EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp2_1;
4852                EM_S[INDEX2(0,0,4)]+=tmp3_1 + tmp4_1 + tmp5_1;
4853                EM_S[INDEX2(3,3,4)]+=tmp3_1;
4854                EM_S[INDEX2(3,0,4)]+=tmp0_1 + tmp6_1;
4855                EM_S[INDEX2(3,1,4)]+=tmp3_1;
4856                EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp7_1;
4857                EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1 + tmp7_1 + tmp8_1;
4858                EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp2_1 + tmp3_1 + tmp8_1;
4859                EM_S[INDEX2(1,3,4)]+=tmp3_1;
4860                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp7_1;
4861                EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1 + tmp9_1;
4862                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp6_1;
4863                EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;
4864                EM_S[INDEX2(1,1,4)]+=tmp3_1;
4865            }
4866            ///////////////
4867            // process B //
4868            ///////////////
4869            if (!B.isEmpty()) {
4870                add_EM_S=true;
4871                const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
4872                const register double B_0 = B_p[0];
4873                const register double B_1 = B_p[1];
4874                const register double tmp2_1 = B_0*w7;
4875                const register double tmp1_1 = B_1*w8;
4876                const register double tmp0_1 = B_0*w6;
4877                const register double tmp3_1 = B_1*w9;
4878                EM_S[INDEX2(1,2,4)]+=tmp0_1;
4879                EM_S[INDEX2(3,2,4)]+=tmp0_1;
4880                EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp2_1;
4881                EM_S[INDEX2(3,0,4)]+=tmp0_1;
4882                EM_S[INDEX2(0,2,4)]+=tmp1_1 + tmp2_1;
4883                EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1;
4884                EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;
4885                EM_S[INDEX2(1,0,4)]+=tmp0_1;
4886            }
4887            ///////////////
4888            // process C //
4889            ///////////////
4890            if (!C.isEmpty()) {
4891                add_EM_S=true;
4892                const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
4893                const register double C_0 = C_p[0];
4894                const register double C_1 = C_p[1];
4895                const register double tmp1_1 = C_1*w8;
4896                const register double tmp3_1 = C_1*w9;
4897                const register double tmp0_1 = C_0*w6;
4898                const register double tmp2_1 = C_0*w7;
4899                EM_S[INDEX2(0,1,4)]+=tmp0_1;
4900                EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp2_1;
4901                EM_S[INDEX2(2,1,4)]+=tmp0_1;
4902                EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;
4903                EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp2_1;
4904                EM_S[INDEX2(2,3,4)]+=tmp0_1;
4905                EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;
4906                EM_S[INDEX2(0,3,4)]+=tmp0_1;
4907            }
4908            ///////////////
4909            // process D //
4910            ///////////////
4911            if (!D.isEmpty()) {
4912                add_EM_S=true;
4913                const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
4914                const register double D_0 = D_p[0];
4915                const register double tmp0_1 = D_0*w10;
4916                EM_S[INDEX2(0,0,4)]+=tmp0_1;
4917                EM_S[INDEX2(0,2,4)]+=tmp0_1;
4918                EM_S[INDEX2(2,0,4)]+=tmp0_1;
4919                EM_S[INDEX2(2,2,4)]+=tmp0_1;
4920            }
4921            ///////////////
4922            // process X //
4923            ///////////////
4924            if (!X.isEmpty()) {
4925                add_EM_F=true;
4926                const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
4927                const register double X_0 = X_p[0];
4928                const register double X_1 = X_p[1];
4929                const register double tmp0_1 = X_0*w11;
4930                const register double tmp2_1 = X_0*w13;
4931                const register double tmp1_1 = X_1*w12;
4932                const register double tmp3_1 = X_1*w14;
4933                EM_F[0]+=tmp0_1 + tmp1_1;
4934                EM_F[1]+=tmp2_1;
4935                EM_F[2]+=tmp0_1 + tmp3_1;
4936                EM_F[3]+=tmp2_1;
4937            }
4938            ///////////////
4939            // process Y //
4940            ///////////////
4941            if (!Y.isEmpty()) {
4942                add_EM_F=true;
4943                const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
4944                const register double Y_0 = Y_p[0];
4945                const register double tmp0_1 = Y_0*w15;
4946                EM_F[0]+=tmp0_1;
4947                EM_F[2]+=tmp0_1;
4948            }
4949                          /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_0 BOTTOM */
4950                        }
4951                        // ADD EM_F  and EM_S
4952            } // end colouring
4953                }
4954    
4955                if (m_faceOffset[1] > -1) {
4956                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
4957    #pragma omp for nowait
4958                        for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4959                  /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_1 TOP */
4960            ///////////////
4961            // process A //
4962            ///////////////
4963            if (!A.isEmpty()) {
4964                add_EM_S=true;
4965                const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
4966                const register double A_00 = A_p[INDEX2(0,0,2)];
4967                const register double A_01 = A_p[INDEX2(0,1,2)];
4968                const register double A_10 = A_p[INDEX2(1,0,2)];
4969                const register double A_11 = A_p[INDEX2(1,1,2)];
4970                const register double tmp0_0 = A_01 + A_10;
4971                const register double tmp7_1 = A_01*w1;
4972                const register double tmp6_1 = tmp0_0*w2;
4973                const register double tmp2_1 = A_10*w2;
4974                const register double tmp0_1 = A_00*w0;
4975                const register double tmp9_1 = tmp0_0*w1;
4976                const register double tmp8_1 = A_11*w5;
4977                const register double tmp1_1 = A_01*w2;
4978                const register double tmp5_1 = A_11*w4;
4979                const register double tmp3_1 = A_10*w1;
4980                const register double tmp4_1 = A_00*w3;
4981                EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
4982                EM_S[INDEX2(1,2,4)]+=tmp0_1 + tmp2_1;
4983                EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp3_1;
4984                EM_S[INDEX2(0,0,4)]+=tmp4_1;
4985                EM_S[INDEX2(3,3,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
4986                EM_S[INDEX2(3,0,4)]+=tmp0_1 + tmp3_1;
4987                EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp4_1 + tmp7_1 + tmp8_1;
4988                EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;
4989                EM_S[INDEX2(0,2,4)]+=tmp4_1;
4990                EM_S[INDEX2(2,0,4)]+=tmp4_1;
4991                EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp3_1 + tmp4_1 + tmp8_1;
4992                EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp7_1;
4993                EM_S[INDEX2(2,2,4)]+=tmp4_1;
4994                EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1;
4995                EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp7_1;
4996                EM_S[INDEX2(1,1,4)]+=tmp4_1 + tmp5_1 + tmp9_1;
4997            }
4998            ///////////////
4999            // process B //
5000            ///////////////
5001            if (!B.isEmpty()) {
5002                add_EM_S=true;
5003                const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
5004                const register double B_0 = B_p[0];
5005                const register double B_1 = B_p[1];
5006                const register double tmp0_1 = B_0*w7;
5007                const register double tmp3_1 = B_1*w8;
5008                const register double tmp1_1 = B_1*w9;
5009                const register double tmp2_1 = B_0*w6;
5010                EM_S[INDEX2(0,1,4)]+=tmp0_1;
5011                EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;
5012                EM_S[INDEX2(3,1,4)]+=tmp1_1 + tmp2_1;
5013                EM_S[INDEX2(2,1,4)]+=tmp0_1;
5014                EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1;
5015                EM_S[INDEX2(2,3,4)]+=tmp0_1;
5016                EM_S[INDEX2(0,3,4)]+=tmp0_1;
5017                EM_S[INDEX2(1,1,4)]+=tmp2_1 + tmp3_1;
5018            }
5019            ///////////////
5020            // process C //
5021            ///////////////
5022            if (!C.isEmpty()) {
5023                add_EM_S=true;
5024                const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
5025                const register double C_0 = C_p[0];
5026                const register double C_1 = C_p[1];
5027                const register double tmp3_1 = C_1*w8;
5028                const register double tmp2_1 = C_0*w6;
5029                const register double tmp1_1 = C_1*w9;
5030                const register double tmp0_1 = C_0*w7;
5031                EM_S[INDEX2(1,2,4)]+=tmp0_1;
5032                EM_S[INDEX2(3,2,4)]+=tmp0_1;
5033                EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;
5034                EM_S[INDEX2(3,0,4)]+=tmp0_1;
5035                EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;
5036                EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;
5037                EM_S[INDEX2(1,0,4)]+=tmp0_1;
5038                EM_S[INDEX2(1,1,4)]+=tmp2_1 + tmp3_1;
5039            }
5040            ///////////////
5041            // process D //
5042            ///////////////
5043            if (!D.isEmpty()) {
5044                add_EM_S=true;
5045                const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
5046                const register double D_0 = D_p[0];
5047                const register double tmp0_1 = D_0*w10;
5048                EM_S[INDEX2(3,3,4)]+=tmp0_1;
5049                EM_S[INDEX2(3,1,4)]+=tmp0_1;
5050                EM_S[INDEX2(1,3,4)]+=tmp0_1;
5051                EM_S[INDEX2(1,1,4)]+=tmp0_1;
5052            }
5053            ///////////////
5054            // process X //
5055            ///////////////
5056            if (!X.isEmpty()) {
5057                add_EM_F=true;
5058                const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
5059                const register double X_0 = X_p[0];
5060                const register double X_1 = X_p[1];
5061                const register double tmp0_1 = X_0*w11;
5062                const register double tmp1_1 = X_0*w13;
5063                const register double tmp2_1 = X_1*w12;
5064                const register double tmp3_1 = X_1*w14;
5065                EM_F[0]+=tmp0_1;
5066                EM_F[1]+=tmp1_1 + tmp2_1;
5067                EM_F[2]+=tmp0_1;
5068                EM_F[3]+=tmp1_1 + tmp3_1;
5069            }
5070            ///////////////
5071            // process Y //
5072            ///////////////
5073            if (!Y.isEmpty()) {
5074                add_EM_F=true;
5075                const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
5076                const register double Y_0 = Y_p[0];
5077                const register double tmp0_1 = Y_0*w15;
5078                EM_F[1]+=tmp0_1;
5079                EM_F[3]+=tmp0_1;
5080            }
5081                          /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_1 BOTTOM */
5082                        }
5083                         // ADD EM_F  and EM_S
5084                    } // end colouring
5085                }
5086    
5087                if (m_faceOffset[2] > -1) {
5088               for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
5089    #pragma omp for nowait
5090                      for (index_t k0 = 0; k0 < m_NE0; ++k0) {
5091                          /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_2 TOP */
5092                          ///////////////
5093                          // process A //
5094                          ///////////////
5095                          if (!A.isEmpty()) {
5096                              add_EM_S=true;
5097                              const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
5098                              const register double A_00 = A_p[INDEX2(0,0,2)];
5099                              const register double A_01 = A_p[INDEX2(0,1,2)];
5100                              const register double A_10 = A_p[INDEX2(1,0,2)];
5101                              const register double A_11 = A_p[INDEX2(1,1,2)];
5102                              const register double tmp0_0 = A_01 + A_10;
5103                              const register double tmp7_1 = A_10*w17;
5104                              const register double tmp3_1 = A_01*w17;
5105                              const register double tmp5_1 = A_00*w21;
5106                              const register double tmp1_1 = A_11*w19;
5107                              const register double tmp2_1 = A_00*w16;
5108                              const register double tmp9_1 = tmp0_0*w18;
5109                              const register double tmp6_1 = tmp0_0*w17;
5110                              const register double tmp8_1 = A_01*w18;
5111                              const register double tmp4_1 = A_11*w20;
5112                              const register double tmp0_1 = A_10*w18;
5113                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
5114                              EM_S[INDEX2(1,2,4)]+=tmp3_1 + tmp4_1;
5115                              EM_S[INDEX2(3,2,4)]+=tmp1_1;
5116                              EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp5_1 + tmp6_1;
5117                              EM_S[INDEX2(3,3,4)]+=tmp1_1;
5118                              EM_S[INDEX2(3,0,4)]+=tmp0_1 + tmp4_1;
5119                              EM_S[INDEX2(3,1,4)]+=tmp4_1 + tmp7_1;
5120                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;
5121                              EM_S[INDEX2(0,2,4)]+=tmp4_1 + tmp8_1;
5122                              EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp4_1;
5123                              EM_S[INDEX2(1,3,4)]+=tmp3_1 + tmp4_1;
5124                              EM_S[INDEX2(2,3,4)]+=tmp1_1;
5125                              EM_S[INDEX2(2,2,4)]+=tmp1_1;
5126                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1 + tmp7_1 + tmp8_1;
5127                              EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp8_1;
5128                              EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp5_1 + tmp9_1;
5129                          }
5130                          ///////////////
5131                          // process B //
5132                          ///////////////
5133                          if (!B.isEmpty()) {
5134                              add_EM_S=true;
5135                              const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
5136                              const register double B_0 = B_p[0];
5137                              const register double B_1 = B_p[1];
5138                              const register double tmp1_1 = B_0*w8;
5139                              const register double tmp2_1 = B_1*w23;
5140                              const register double tmp0_1 = B_1*w22;
5141                              const register double tmp3_1 = B_0*w9;
5142                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
5143                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
5144                              EM_S[INDEX2(3,0,4)]+=tmp2_1;
5145                              EM_S[INDEX2(3,1,4)]+=tmp2_1;
5146                              EM_S[INDEX2(2,1,4)]+=tmp2_1;
5147                              EM_S[INDEX2(2,0,4)]+=tmp2_1;
5148                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp3_1;
5149                              EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp3_1;
5150                          }
5151                          ///////////////
5152                          // process C //
5153                          ///////////////
5154                          if (!C.isEmpty()) {
5155                              add_EM_S=true;
5156                              const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
5157                              const register double C_0 = C_p[0];
5158                              const register double C_1 = C_p[1];
5159                              const register double tmp2_1 = C_1*w23;
5160                              const register double tmp0_1 = C_1*w22;
5161                              const register double tmp1_1 = C_0*w9;
5162                              const register double tmp3_1 = C_0*w8;
5163                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
5164                              EM_S[INDEX2(1,2,4)]+=tmp2_1;
5165                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp3_1;
5166                              EM_S[INDEX2(0,2,4)]+=tmp2_1;
5167                              EM_S[INDEX2(1,3,4)]+=tmp2_1;
5168                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp3_1;
5169                              EM_S[INDEX2(0,3,4)]+=tmp2_1;
5170                              EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;
5171                          }
5172                          ///////////////
5173                          // process D //
5174                          ///////////////
5175                          if (!D.isEmpty()) {
5176                              add_EM_S=true;
5177                              const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
5178                              const register double D_0 = D_p[0];
5179                              const register double tmp0_1 = D_0*w24;
5180                              EM_S[INDEX2(0,1,4)]+=tmp0_1;
5181                              EM_S[INDEX2(0,0,4)]+=tmp0_1;
5182                              EM_S[INDEX2(1,0,4)]+=tmp0_1;
5183                              EM_S[INDEX2(1,1,4)]+=tmp0_1;
5184                          }
5185                          ///////////////
5186                          // process X //
5187                          ///////////////
5188                          if (!X.isEmpty()) {
5189                              add_EM_F=true;
5190                              const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
5191                              const register double X_0 = X_p[0];
5192                              const register double X_1 = X_p[1];
5193                              const register double tmp1_1 = X_0*w12;
5194                              const register double tmp2_1 = X_0*w14;
5195                              const register double tmp0_1 = X_1*w25;
5196                              const register double tmp3_1 = X_1*w26;
5197                              EM_F[0]+=tmp0_1 + tmp1_1;
5198                              EM_F[1]+=tmp0_1 + tmp2_1;
5199                              EM_F[2]+=tmp3_1;
5200                              EM_F[3]+=tmp3_1;
5201                          }
5202                          ///////////////
5203                          // process Y //
5204                          ///////////////
5205                          if (!Y.isEmpty()) {
5206                              add_EM_F=true;
5207                              const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
5208                              const register double Y_0 = Y_p[0];
5209                              const register double tmp0_1 = Y_0*w27;
5210                              EM_F[0]+=tmp0_1;
5211                              EM_F[1]+=tmp0_1;
5212                          }
5213                          /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_2 BOTTOM */
5214                      }
5215                      // ADD EM_F  and EM_S
5216                   } // end colouring
5217                }
5218    
5219                if (m_faceOffset[3] > -1) {
5220                 for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
5221    #pragma omp for nowait
5222                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
5223                           /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_3 TOP */
5224                           ///////////////
5225                           // process A //
5226                           ///////////////
5227                           if (!A.isEmpty()) {
5228                               add_EM_S=true;
5229                               const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
5230                               const register double A_00 = A_p[INDEX2(0,0,2)];
5231                               const register double A_01 = A_p[INDEX2(0,1,2)];
5232                               const register double A_10 = A_p[INDEX2(1,0,2)];
5233                               const register double A_11 = A_p[INDEX2(1,1,2)];
5234                               const register double tmp0_0 = A_01 + A_10;
5235                               const register double tmp2_1 = A_10*w17;
5236                               const register double tmp5_1 = A_01*w17;
5237                               const register double tmp6_1 = A_00*w21;
5238                               const register double tmp0_1 = A_11*w19;
5239                               const register double tmp4_1 = A_00*w16;
5240                               const register double tmp9_1 = tmp0_0*w18;
5241                               const register double tmp7_1 = tmp0_0*w17;
5242                               const register double tmp8_1 = A_01*w18;
5243                               const register double tmp1_1 = A_11*w20;
5244                               const register double tmp3_1 = A_10*w18;
5245                               EM_S[INDEX2(0,1,4)]+=tmp0_1;
5246                               EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
5247                               EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp3_1 + tmp4_1 + tmp5_1;
5248                               EM_S[INDEX2(0,0,4)]+=tmp0_1;
5249                               EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
5250                               EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp8_1;
5251                               EM_S[INDEX2(3,1,4)]+=tmp1_1 + tmp8_1;
5252                               EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp5_1;
5253                               EM_S[INDEX2(0,2,4)]+=tmp1_1 + tmp2_1;
5254                               EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp5_1;
5255                               EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp3_1;
5256                               EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1 + tmp4_1 + tmp8_1;
5257                               EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp9_1;
5258                               EM_S[INDEX2(1,0,4)]+=tmp0_1;
5259                               EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp3_1;
5260                               EM_S[INDEX2(1,1,4)]+=tmp0_1;
5261                           }
5262                           ///////////////
5263                           // process B //
5264                           ///////////////
5265                           if (!B.isEmpty()) {
5266                               add_EM_S=true;
5267                               const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
5268                               const register double B_0 = B_p[0];
5269                               const register double B_1 = B_p[1];
5270                               const register double tmp3_1 = B_0*w8;
5271                               const register double tmp2_1 = B_1*w23;
5272                               const register double tmp0_1 = B_1*w22;
5273                               const register double tmp1_1 = B_0*w9;
5274                               EM_S[INDEX2(1,2,4)]+=tmp0_1;
5275                               EM_S[INDEX2(3,2,4)]+=tmp1_1 + tmp2_1;
5276                               EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;
5277                               EM_S[INDEX2(0,2,4)]+=tmp0_1;
5278                               EM_S[INDEX2(1,3,4)]+=tmp0_1;
5279                               EM_S[INDEX2(2,3,4)]+=tmp2_1 + tmp3_1;
5280                               EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;
5281                               EM_S[INDEX2(0,3,4)]+=tmp0_1;
5282                           }
5283                           ///////////////
5284                           // process C //
5285                           ///////////////
5286                           if (!C.isEmpty()) {
5287                               add_EM_S=true;
5288                               const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
5289                               const register double C_0 = C_p[0];
5290                               const register double C_1 = C_p[1];
5291                               const register double tmp1_1 = C_1*w23;
5292                               const register double tmp3_1 = C_1*w22;
5293                               const register double tmp0_1 = C_0*w8;
5294                               const register double tmp2_1 = C_0*w9;
5295                               EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1;
5296                               EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;
5297                               EM_S[INDEX2(3,0,4)]+=tmp3_1;
5298                               EM_S[INDEX2(3,1,4)]+=tmp3_1;
5299                               EM_S[INDEX2(2,1,4)]+=tmp3_1;
5300                               EM_S[INDEX2(2,0,4)]+=tmp3_1;
5301                               EM_S[INDEX2(2,3,4)]+=tmp1_1 + tmp2_1;
5302                               EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp1_1;
5303                           }
5304                           ///////////////
5305                           // process D //
5306                           ///////////////
5307                           if (!D.isEmpty()) {
5308                               add_EM_S=true;
5309                               const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
5310                               const register double D_0 = D_p[0];
5311                               const register double tmp0_1 = D_0*w24;
5312                               EM_S[INDEX2(3,2,4)]+=tmp0_1;
5313                               EM_S[INDEX2(3,3,4)]+=tmp0_1;
5314                               EM_S[INDEX2(2,3,4)]+=tmp0_1;
5315                               EM_S[INDEX2(2,2,4)]+=tmp0_1;
5316                           }
5317                           ///////////////
5318                           // process X //
5319                           ///////////////
5320                           if (!X.isEmpty()) {
5321                               add_EM_F=true;
5322                               const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
5323                               const register double X_0 = X_p[0];
5324                               const register double X_1 = X_p[1];
5325                               const register double tmp2_1 = X_0*w12;
5326                               const register double tmp1_1 = X_1*w26;
5327                               const register double tmp0_1 = X_1*w25;
5328                               const register double tmp3_1 = X_0*w14;
5329                               EM_F[0]+=tmp0_1;
5330                               EM_F[1]+=tmp0_1;
5331                               EM_F[2]+=tmp1_1 + tmp2_1;
5332                               EM_F[3]+=tmp1_1 + tmp3_1;
5333                           }
5334                           ///////////////
5335                           // process Y //
5336                           ///////////////
5337                           if (!Y.isEmpty()) {
5338                               add_EM_F=true;
5339                               const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
5340                               const register double Y_0 = Y_p[0];
5341                               const register double tmp0_1 = Y_0*w27;
5342                               EM_F[2]+=tmp0_1;
5343                               EM_F[3]+=tmp0_1;
5344                           }
5345                          /* GENERATOR SNIP_PDEBC_SINGLE_REDUCED_3 BOTTOM */
5346                    }
5347                    // ADD EM_F  and EM_S
5348                    } // end colouring
5349                }
5350            } // end of parallel section
5351    
5352  }  }
5353    
5354  //protected  //protected
5355  void Rectangle::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
5356        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
5357            const escript::Data& C, const escript::Data& D,
5358            const escript::Data& X, const escript::Data& Y) const
5359  {  {
5360        const double h0 = m_l0/m_gNE0;
5361        const double h1 = m_l1/m_gNE1;
5362        dim_t numEq, numComp;
5363        if (!mat) {
5364            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
5365        } else {
5366            numEq=mat->logical_row_block_size;
5367            numComp=mat->logical_col_block_size;
5368        }
5369        /* GENERATOR SNIP_PDEBC_SYSTEM_PRE TOP */
5370        const double w0 = -0.31100423396407310779*h1/pow(h0,2);
5371        const double w1 = -0.39433756729740644113/h0;
5372        const double w10 = 0.083333333333333333333*h1/pow(h0,2);
5373        const double w11 = -0.5/h1;
5374        const double w12 = -0.33333333333333333333*h1/pow(h0,2);
5375        const double w13 = -0.5/h0;
5376        const double w14 = -0.16666666666666666667*h1/pow(h0,2);
5377        const double w15 = 0.5/h0;
5378        const double w16 = 0.33333333333333333333*h1/pow(h0,2);
5379        const double w17 = 1.0/h1;
5380        const double w18 = 0.16666666666666666667*h1/pow(h0,2);
5381        const double w19 = -1.0/h1;
5382        const double w2 = -0.022329099369260225539*h1/pow(h0,2);
5383        const double w20 = 0.083333333333333333333*h1/h0;
5384        const double w21 = 0.022329099369260225539*h1/h0;
5385        const double w22 = 0.31100423396407310779*h1/h0;
5386        const double w23 = -0.31100423396407310779*h1/h0;
5387        const double w24 = -0.39433756729740643311;
5388        const double w25 = -0.022329099369260225539*h1/h0;
5389        const double w26 = -0.10566243270259356689;
5390        const double w27 = -0.083333333333333333333*h1/h0;
5391        const double w28 = 0.39433756729740643311;
5392        const double w29 = 0.10566243270259356689;
5393        const double w3 = -0.10566243270259355887/h0;
5394        const double w30 = 0.16666666666666666667*h1/h0;
5395        const double w31 = 0.33333333333333333333*h1/h0;
5396        const double w32 = -0.33333333333333333333*h1/h0;
5397        const double w33 = -0.50000000000000000000;
5398        const double w34 = -0.16666666666666666667*h1/h0;
5399        const double w35 = 0.50000000000000000000;
5400        const double w36 = 0.31100423396407310779*h1;
5401        const double w37 = 0.022329099369260225539*h1;
5402        const double w38 = 0.083333333333333333333*h1;
5403        const double w39 = 0.33333333333333333333*h1;
5404        const double w4 = -0.083333333333333333333*h1/pow(h0,2);
5405        const double w40 = 0.16666666666666666667*h1;
5406        const double w41 = -0.39433756729740644113*h1/h0;
5407        const double w42 = -0.10566243270259355887*h1/h0;
5408        const double w43 = 0.39433756729740644113*h1/h0;
5409        const double w44 = 0.10566243270259355887*h1/h0;
5410        const double w45 = -0.5*h1/h0;
5411        const double w46 = -1.0000000000000000000;
5412        const double w47 = 0.5*h1/h0;
5413        const double w48 = 1.0000000000000000000;
5414        const double w49 = 0.39433756729740644113*h1;
5415        const double w5 = 0.39433756729740644113/h0;
5416        const double w50 = 0.10566243270259355887*h1;
5417        const double w51 = 0.5*h1;
5418        const double w52 = -0.5/h0;
5419        const double w53 = 0.10566243270259355887/h1;
5420        const double w54 = -0.39433756729740644113/h1;
5421        const double w55 = 0.083333333333333333333*h0/pow(h1,2);
5422        const double w56 = 0.39433756729740644113/h1;
5423        const double w57 = -0.10566243270259355887/h1;
5424        const double w58 = -0.083333333333333333333*h0/pow(h1,2);
5425        const double w59 = 0.5/h0;
5426        const double w6 = 0.10566243270259355887/h0;
5427        const double w60 = 0.31100423396407310779*h0/pow(h1,2);
5428        const double w61 = 0.022329099369260225539*h0/pow(h1,2);
5429        const double w62 = -0.022329099369260225539*h0/pow(h1,2);
5430        const double w63 = -0.31100423396407310779*h0/pow(h1,2);
5431        const double w64 = -1.0/h0;
5432        const double w65 = 0.5/h1;
5433        const double w66 = -0.5/h1;
5434        const double w67 = 0.16666666666666666667*h0/pow(h1,2);
5435        const double w68 = -0.16666666666666666667*h0/pow(h1,2);
5436        const double w69 = 1.0/h0;
5437        const double w7 = 0.31100423396407310779*h1/pow(h0,2);
5438        const double w70 = 0.33333333333333333333*h0/pow(h1,2);
5439        const double w71 = -0.33333333333333333333*h0/pow(h1,2);
5440        const double w72 = -0.083333333333333333333*h0/h1;
5441        const double w73 = -0.31100423396407310779*h0/h1;
5442        const double w74 = -0.022329099369260225539*h0/h1;
5443        const double w75 = 0.083333333333333333333*h0/h1;
5444        const double w76 = 0.022329099369260225539*h0/h1;
5445        const double w77 = 0.31100423396407310779*h0/h1;
5446        const double w78 = -0.16666666666666666667*h0/h1;
5447        const double w79 = -0.33333333333333333333*h0/h1;
5448        const double w8 = 0.5/h1;
5449        const double w80 = 0.16666666666666666667*h0/h1;
5450        const double w81 = 0.33333333333333333333*h0/h1;
5451        const double w82 = 0.083333333333333333333*h0;
5452        const double w83 = 0.31100423396407310779*h0;
5453        const double w84 = 0.022329099369260225539*h0;
5454        const double w85 = 0.16666666666666666667*h0;
5455        const double w86 = 0.33333333333333333333*h0;
5456        const double w87 = -0.39433756729740644113*h0/h1;
5457        const double w88 = -0.10566243270259355887*h0/h1;
5458        const double w89 = 0.39433756729740644113*h0/h1;
5459        const double w9 = 0.022329099369260225539*h1/pow(h0,2);
5460        const double w90 = 0.10566243270259355887*h0/h1;
5461        const double w91 = -0.5*h0/h1;
5462        const double w92 = 0.5*h0/h1;
5463        const double w93 = 0.39433756729740644113*h0;
5464        const double w94 = 0.10566243270259355887*h0;
5465        const double w95 = 0.5*h0;
5466      /* GENERATOR SNIP_PDEBC_SYSTEM_PRE BOTTOM */
5467      #pragma omp parallel
5468            {
5469                if (m_faceOffset[0] > -1) {
5470                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
5471    #pragma omp for nowait
5472                        for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
5473                  /* GENERATOR SNIP_PDEBC_SYSTEM_0 TOP */
5474            ///////////////
5475            // process A //
5476            ///////////////
5477            if (!A.isEmpty()) {
5478                add_EM_S=true;
5479                const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
5480                if (A.actsExpanded()) {
5481                    for (index_t k=0; k<numEq; k++) {
5482                        for (index_t m=0; m<numComp; m++) {
5483                            const register double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];
5484                            const register double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];
5485                            const register double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];
5486                            const register double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];
5487                            const register double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];
5488                            const register double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];
5489                            const register double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];
5490                            const register double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];
5491                            const register double tmp1_0 = A_11_0 + A_11_1;
5492                            const register double tmp2_0 = A_01_1 + A_10_1;
5493                            const register double tmp0_0 = A_00_0 + A_00_1;
5494                            const register double tmp3_0 = A_01_0 + A_10_0;
5495                            const register double tmp7_1 = A_00_1*w0;
5496                            const register double tmp16_1 = A_00_0*w9;
5497                            const register double tmp23_1 = A_01_0*w1;
5498                            const register double tmp4_1 = A_01_1*w6;
5499                            const register double tmp13_1 = tmp2_0*w6;
5500                            const register double tmp30_1 = A_10_0*w3;
5501                            const register double tmp2_1 = A_10_0*w1;
5502                            const register double tmp14_1 = A_00_0*w7;
5503                            const register double tmp5_1 = tmp0_0*w4;
5504                            const register double tmp15_1 = tmp3_0*w5;
5505                            const register double tmp22_1 = A_10_0*w5;
5506                            const register double tmp10_1 = A_00_0*w2;
5507                            const register double tmp20_1 = tmp0_0*w10;
5508                            const register double tmp29_1 = tmp2_0*w1;
5509                            const register double tmp0_1 = A_10_1*w3;
5510                            const register double tmp31_1 = A_10_1*w1;
5511                            const register double tmp11_1 = tmp1_0*w8;
5512                            const register double tmp21_1 = A_10_1*w6;
5513                            const register double tmp9_1 = A_01_1*w5;
5514                            const register double tmp12_1 = A_00_1*w9;
5515                            const register double tmp27_1 = A_10_1*w5;
5516                            const register double tmp6_1 = A_01_0*w5;
5517                            const register double tmp8_1 = A_01_0*w6;
5518                            const register double tmp1_1 = A_00_0*w0;
5519                            const register double tmp3_1 = A_00_1*w2;
5520                            const register double tmp19_1 = A_01_1*w1;
5521                            const register double tmp25_1 = A_01_1*w3;
5522                            const register double tmp18_1 = A_01_0*w3;
5523                            const register double tmp24_1 = tmp1_0*w11;
5524                            const register double tmp26_1 = A_10_0*w6;
5525                            const register double tmp17_1 = A_00_1*w7;
5526                            const register double tmp28_1 = tmp3_0*w3;
5527                            EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
5528                            EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
5529                            EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp7_1 + tmp8_1 + tmp9_1;
5530                            EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
5531                            EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp16_1 + tmp17_1;
5532                            EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp18_1 + tmp19_1 + tmp5_1;
5533                            EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp20_1;
5534                            EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp21_1 + tmp22_1 + tmp5_1;
5535                            EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp20_1 + tmp23_1 + tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1;
5536                            EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp20_1 + tmp24_1 + tmp2_1 + tmp8_1 + tmp9_1;
5537                            EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp20_1;
5538                            EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp10_1 + tmp26_1 + tmp27_1 + tmp7_1;
5539                            EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp11_1 + tmp16_1 + tmp17_1 + tmp28_1 + tmp29_1;
5540                            EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp23_1 + tmp25_1 + tmp3_1;
5541                            EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp30_1 + tmp31_1 + tmp5_1;
5542                            EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp12_1 + tmp14_1;
5543                        }
5544                    }
5545                } else { /* constant data */
5546                    for (index_t k=0; k<numEq; k++) {
5547                              for (index_t m=0; m<numComp; m++) {
5548                                  const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
5549                            const register double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
5550                            const register double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
5551                            const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
5552                            const register double tmp0_0 = A_01 + A_10;
5553                            const register double tmp5_1 = A_11*w17;
5554                            const register double tmp8_1 = A_00*w18;
5555                            const register double tmp0_1 = A_10*w13;
5556                            const register double tmp11_1 = tmp0_0*w13;
5557                            const register double tmp3_1 = A_01*w15;
5558                            const register double tmp10_1 = A_11*w19;
5559                            const register double tmp4_1 = A_00*w16;
5560                            const register double tmp2_1 = A_00*w14;
5561                            const register double tmp6_1 = tmp0_0*w15;
5562                            const register double tmp1_1 = A_00*w12;
5563                            const register double tmp7_1 = A_01*w13;
5564                            const register double tmp9_1 = A_10*w15;
5565                            EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
5566                            EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
5567                            EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
5568                            EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
5569                            EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp4_1;
5570                            EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp7_1;
5571                            EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp8_1;
5572                            EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp2_1 + tmp9_1;
5573                            EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp10_1 + tmp7_1 + tmp8_1 + tmp9_1;
5574                            EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp3_1 + tmp8_1;
5575                            EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp8_1;
5576                            EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp1_1 + tmp9_1;
5577                            EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp11_1 + tmp4_1 + tmp5_1;
5578                            EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp7_1;
5579                            EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
5580                            EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp4_1;
5581                        }
5582                    }
5583                }
5584            }
5585            ///////////////
5586            // process B //
5587            ///////////////
5588            if (!B.isEmpty()) {
5589                add_EM_S=true;
5590                const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
5591                if (B.actsExpanded()) {
5592                    for (index_t k=0; k<numEq; k++) {
5593                        for (index_t m=0; m<numComp; m++) {
5594                            const register double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];
5595                            const register double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];
5596                            const register double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];
5597                            const register double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];
5598                            const register double tmp0_0 = B_0_0 + B_0_1;
5599                            const register double tmp7_1 = tmp0_0*w27;
5600                            const register double tmp5_1 = B_1_0*w24;
5601                            const register double tmp8_1 = B_1_0*w26;
5602                            const register double tmp11_1 = B_1_0*w28;
5603                            const register double tmp10_1 = B_1_1*w29;
5604                            const register double tmp15_1 = B_0_1*w23;
5605                            const register double tmp0_1 = tmp0_0*w20;
5606                            const register double tmp4_1 = B_0_1*w25;
5607                            const register double tmp9_1 = B_1_1*w24;
5608                            const register double tmp16_1 = B_0_0*w22;
5609                            const register double tmp3_1 = B_1_1*w26;
5610                            const register double tmp13_1 = B_1_1*w28;
5611                            const register double tmp12_1 = B_1_0*w29;
5612                            const register double tmp1_1 = B_0_1*w22;
5613                            const register double tmp17_1 = B_0_1*w21;
5614                            const register double tmp6_1 = B_0_0*w23;
5615                            const register double tmp2_1 = B_0_0*w21;
5616                            const register double tmp14_1 = B_0_0*w25;
5617                            EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0_1;
5618                            EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
5619                            EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1;
5620                            EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0_1;
5621                            EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
5622                            EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp7_1;
5623                            EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;
5624                            EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1;
5625                        }
5626                    }
5627                } else { /* constant data */
5628                    for (index_t k=0; k<numEq; k++) {
5629                        for (index_t m=0; m<numComp; m++) {
5630                            const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
5631                            const register double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
5632                            const register double tmp5_1 = B_1*w35;
5633                            const register double tmp1_1 = B_0*w31;
5634                            const register double tmp3_1 = B_1*w33;
5635                            const register double tmp0_1 = B_0*w30;
5636                            const register double tmp2_1 = B_0*w32;
5637                            const register double tmp4_1 = B_0*w34;
5638                            EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0_1;
5639                            EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp1_1;
5640                            EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
5641                            EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0_1;
5642                            EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
5643                            EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp4_1 + tmp5_1;
5644                            EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp2_1 + tmp5_1;
5645                            EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1;
5646                        }
5647                    }
5648                }
5649            }
5650            ///////////////
5651            // process C //
5652            ///////////////
5653            if (!C.isEmpty()) {
5654                add_EM_S=true;
5655                const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
5656                if (C.actsExpanded()) {
5657                    for (index_t k=0; k<numEq; k++) {
5658                        for (index_t m=0; m<numComp; m++) {
5659                            const register double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];
5660                            const register double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];
5661                            const register double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];
5662                            const register double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];
5663                            const register double tmp0_0 = C_0_0 + C_0_1;
5664                            const register double tmp14_1 = C_1_0*w29;
5665                            const register double tmp9_1 = tmp0_0*w27;
5666                            const register double tmp11_1 = C_1_1*w24;
5667                            const register double tmp2_1 = C_1_1*w26;
5668                            const register double tmp15_1 = C_1_1*w28;
5669                            const register double tmp12_1 = C_0_1*w22;
5670                            const register double tmp0_1 = C_0_0*w22;
5671                            const register double tmp6_1 = tmp0_0*w20;
5672                            const register double tmp17_1 = C_0_1*w23;
5673                            const register double tmp3_1 = C_0_1*w25;
5674                            const register double tmp16_1 = C_0_0*w25;
5675                            const register double tmp4_1 = C_1_0*w24;
5676                            const register double tmp10_1 = C_1_0*w26;
5677                            const register double tmp8_1 = C_1_0*w28;
5678                            const register double tmp7_1 = C_1_1*w29;
5679                            const register double tmp5_1 = C_0_0*w23;
5680                            const register double tmp1_1 = C_0_1*w21;
5681                            const register double tmp13_1 = C_0_0*w21;
5682                            EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
5683                            EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
5684                            EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp6_1;
5685                            EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
5686                            EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp9_1;
5687                            EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp12_1 + tmp13_1;
5688                            EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1;
5689                            EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp6_1;
5690                        }
5691                    }
5692                } else { /* constant data */
5693                    for (index_t k=0; k<numEq; k++) {
5694                        for (index_t m=0; m<numComp; m++) {
5695                            const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
5696                            const register double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
5697                            const register double tmp3_1 = C_0*w30;
5698                            const register double tmp1_1 = C_0*w32;
5699                            const register double tmp5_1 = C_0*w34;
5700                            const register double tmp2_1 = C_1*w33;
5701                            const register double tmp4_1 = C_1*w35;
5702                            const register double tmp0_1 = C_0*w31;
5703                            EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
5704                            EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
5705                            EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp3_1;
5706                            EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1;
5707                            EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1 + tmp5_1;
5708                            EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
5709                            EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp1_1 + tmp4_1;
5710                            EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp3_1;
5711                        }
5712                    }
5713                }
5714            }
5715            ///////////////
5716            // process D //
5717            ///////////////
5718            if (!D.isEmpty()) {
5719                add_EM_S=true;
5720                const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
5721                if (D.actsExpanded()) {
5722                    for (index_t k=0; k<numEq; k++) {
5723                        for (index_t m=0; m<numComp; m++) {
5724                            const register double D_0 = D_p[INDEX3(k, m, 0, numEq, numComp)];
5725                            const register double D_1 = D_p[INDEX3(k, m, 1, numEq, numComp)];
5726                            const register double tmp0_0 = D_0 + D_1;
5727                            const register double tmp1_1 = D_1*w37;
5728                            const register double tmp4_1 = D_0*w37;
5729                            const register double tmp0_1 = D_0*w36;
5730                            const register double tmp3_1 = D_1*w36;
5731                            const register double tmp2_1 = tmp0_0*w38;
5732                            EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
5733                            EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1;
5734                            EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1;
5735                            EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
5736                        }
5737                     }
5738                } else { /* constant data */
5739                    for (index_t k=0; k<numEq; k++) {
5740                        for (index_t m=0; m<numComp; m++) {
5741                            const register double D_0 = D_p[INDEX2(k, m, numEq)];
5742                            const register double tmp0_1 = D_0*w39;
5743                            const register double tmp1_1 = D_0*w40;
5744                            EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
5745                            EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp1_1;
5746                            EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1;
5747                            EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;
5748                        }
5749                    }
5750                }
5751            }
5752            ///////////////
5753            // process X //
5754            ///////////////
5755            if (!X.isEmpty()) {
5756                add_EM_F=true;
5757                const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
5758                if (X.actsExpanded()) {
5759                    for (index_t k=0; k<numEq; k++) {
5760                        const register double X_0_0 = X_p[INDEX3(k, 0, 0, numEq, 2)];
5761                        const register double X_1_0 = X_p[INDEX3(k, 1, 0, numEq, 2)];
5762                        const register double X_0_1 = X_p[INDEX3(k, 0, 1, numEq, 2)];
5763                        const register double X_1_1 = X_p[INDEX3(k, 1, 1, numEq, 2)];
5764                        const register double tmp0_0 = X_1_0 + X_1_1;
5765                        const register double tmp3_1 = X_0_1*w44;
5766                        const register double tmp1_1 = X_0_1*w42;
5767                        const register double tmp8_1 = X_0_0*w44;
5768                        const register double tmp4_1 = X_0_0*w43;
5769                        const register double tmp7_1 = X_0_0*w42;
5770                        const register double tmp5_1 = tmp0_0*w35;
5771                        const register double tmp0_1 = X_0_0*w41;
5772                        const register double tmp9_1 = X_0_1*w43;
5773                        const register double tmp2_1 = tmp0_0*w33;
5774                        const register double tmp6_1 = X_0_1*w41;
5775                        EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;
5776                        EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1;
5777                        EM_F[INDEX2(k,2,numEq)]+=tmp5_1 + tmp6_1 + tmp7_1;
5778                        EM_F[INDEX2(k,3,numEq)]+=tmp8_1 + tmp9_1;
5779                    }
5780                } else { /* constant data */
5781                    for (index_t k=0; k<numEq; k++) {
5782                        const register double X_0 = X_p[INDEX2(k, 0, numEq)];
5783                        const register double X_1 = X_p[INDEX2(k, 1, numEq)];
5784                        const register double tmp2_1 = X_0*w47;
5785                        const register double tmp1_1 = X_0*w45;
5786                        const register double tmp3_1 = X_1*w48;
5787                        const register double tmp0_1 = X_1*w46;
5788                        EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
5789                        EM_F[INDEX2(k,1,numEq)]+=tmp2_1;
5790                        EM_F[INDEX2(k,2,numEq)]+=tmp1_1 + tmp3_1;
5791                        EM_F[INDEX2(k,3,numEq)]+=tmp2_1;
5792                    }
5793                }
5794            }
5795            ///////////////
5796            // process Y //
5797            ///////////////
5798            if (!Y.isEmpty()) {
5799                add_EM_F=true;
5800                const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
5801                if (Y.actsExpanded()) {
5802                    for (index_t k=0; k<numEq; k++) {
5803                        const register double Y_0 = Y_p[INDEX2(k, 0, numEq)];
5804                        const register double Y_1 = Y_p[INDEX2(k, 1, numEq)];
5805                        const register double tmp0_1 = Y_0*w49;
5806                        const register double tmp1_1 = Y_1*w50;
5807                        const register double tmp3_1 = Y_0*w50;
5808                        const register double tmp2_1 = Y_1*w49;
5809                        EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
5810                        EM_F[INDEX2(k,2,numEq)]+=tmp2_1 + tmp3_1;
5811                    }
5812                } else { /* constant data */
5813                    for (index_t k=0; k<numEq; k++) {
5814                        const register double Y_0 = Y_p[k];
5815                        const register double tmp0_1 = Y_0*w51;
5816                        EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
5817                        EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
5818                    }
5819                }
5820            }
5821                          /* GENERATOR SNIP_PDEBC_SYSTEM_0 BOTTOM */
5822                        }
5823                        // ADD EM_F  and EM_S
5824            } // end colouring
5825                }
5826    
5827                if (m_faceOffset[1] > -1) {
5828                for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
5829    #pragma omp for nowait
5830                        for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
5831                  /* GENERATOR SNIP_PDEBC_SYSTEM_1 TOP */
5832            ///////////////
5833            // process A //
5834            ///////////////
5835            if (!A.isEmpty()) {
5836                add_EM_S=true;
5837                const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
5838                if (A.actsExpanded()) {
5839                    for (index_t k=0; k<numEq; k++) {
5840                        for (index_t m=0; m<numComp; m++) {
5841                            const register double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];
5842                            const register double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];
5843                            const register double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];
5844                            const register double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];
5845                            const register double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];
5846                            const register double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];
5847                            const register double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];
5848                            const register double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];
5849                            const register double tmp1_0 = A_11_0 + A_11_1;
5850                            const register double tmp3_0 = A_01_1 + A_10_1;
5851                            const register double tmp0_0 = A_00_0 + A_00_1;
5852                            const register double tmp2_0 = A_01_0 + A_10_0;
5853                            const register double tmp8_1 = A_00_1*w0;
5854                            const register double tmp14_1 = A_00_0*w9;
5855                            const register double tmp29_1 = A_01_0*w1;
5856                            const register double tmp1_1 = A_01_1*w6;
5857                            const register double tmp15_1 = tmp2_0*w6;
5858                            const register double tmp7_1 = A_10_0*w3;
5859                            const register double tmp19_1 = A_10_0*w1;
5860                            const register double tmp12_1 = A_00_0*w7;
5861                            const register double tmp5_1 = tmp0_0*w4;
5862                            const register double tmp17_1 = tmp3_0*w5;
5863                            const register double tmp25_1 = A_10_0*w5;
5864                            const register double tmp10_1 = A_00_0*w2;
5865                            const register double tmp23_1 = tmp0_0*w10;
5866                            const register double tmp31_1 = tmp2_0*w1;
5867                            const register double tmp18_1 = A_10_1*w3;
5868                            const register double tmp9_1 = A_10_1*w1;
5869                            const register double tmp13_1 = tmp1_0*w8;
5870                            const register double tmp24_1 = A_10_1*w6;
5871                            const register double tmp27_1 = A_01_1*w5;
5872                            const register double tmp11_1 = A_00_1*w9;
5873                            const register double tmp6_1 = A_10_1*w5;
5874                            const register double tmp2_1 = A_01_0*w5;
5875                            const register double tmp26_1 = A_01_0*w6;
5876                            const register double tmp0_1 = A_00_0*w0;
5877                            const register double tmp3_1 = A_00_1*w2;
5878                            const register double tmp20_1 = A_01_1*w1;
5879                            const register double tmp28_1 = A_01_1*w3;
5880                            const register double tmp22_1 = A_01_0*w3;
5881                            const register double tmp21_1 = tmp1_0*w11;
5882                            const register double tmp4_1 = A_10_0*w6;
5883                            const register double tmp16_1 = A_00_1*w7;
5884                            const register double tmp30_1 = tmp3_0*w3;
5885                            EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
5886                            EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
5887                            EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp7_1 + tmp8_1 + tmp9_1;
5888                            EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1;
5889                            EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1;
5890                            EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp18_1 + tmp19_1 + tmp5_1;
5891                            EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;
5892                            EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp26_1 + tmp27_1 + tmp5_1;
5893                            EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp23_1;
5894                            EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp23_1;
5895                            EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1 + tmp21_1 + tmp23_1 + tmp2_1 + tmp7_1 + tmp9_1;
5896                            EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp10_1 + tmp20_1 + tmp22_1 + tmp8_1;
5897                            EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp14_1 + tmp16_1;
5898                            EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp24_1 + tmp25_1 + tmp3_1;
5899                            EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp28_1 + tmp29_1 + tmp5_1;
5900                            EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp13_1 + tmp30_1 + tmp31_1;
5901                        }
5902                    }
5903                } else { /* constant data */
5904                    for (index_t k=0; k<numEq; k++) {
5905                              for (index_t m=0; m<numComp; m++) {
5906                                  const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
5907                            const register double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
5908                            const register double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
5909                            const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
5910                            const register double tmp0_0 = A_01 + A_10;
5911                            const register double tmp6_1 = A_11*w17;
5912                            const register double tmp9_1 = A_00*w18;
5913                            const register double tmp4_1 = A_10*w13;
5914                            const register double tmp11_1 = tmp0_0*w13;
5915                            const register double tmp0_1 = A_01*w15;
5916                            const register double tmp10_1 = A_11*w19;
5917                            const register double tmp5_1 = A_00*w16;
5918                            const register double tmp2_1 = A_00*w14;
5919                            const register double tmp7_1 = tmp0_0*w15;
5920                            const register double tmp1_1 = A_00*w12;
5921                            const register double tmp8_1 = A_01*w13;
5922                            const register double tmp3_1 = A_10*w15;
5923                            EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
5924                            EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
5925                            EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp1_1 + tmp4_1;
5926                            EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1;
5927                            EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp7_1;
5928                            EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;
5929                            EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp3_1 + tmp8_1 + tmp9_1;
5930                            EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
5931                            EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp9_1;
5932                            EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp9_1;
5933                            EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp4_1 + tmp9_1;
5934                            EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp1_1 + tmp8_1;
5935                            EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp5_1;
5936                            EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
5937                            EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp2_1 + tmp8_1;
5938                            EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp11_1 + tmp5_1 + tmp6_1;
5939                        }
5940                    }
5941                }
5942            }
5943            ///////////////
5944            // process B //
5945            ///////////////
5946            if (!B.isEmpty()) {
5947                add_EM_S=true;
5948                const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
5949                if (B.actsExpanded()) {
5950                    for (index_t k=0; k<numEq; k++) {
5951                        for (index_t m=0; m<numComp; m++) {
5952                            const register double B_0_0 = B_p[INDEX4(k,0,m,0, numEq,2,numComp)];
5953                            const register double B_1_0 = B_p[INDEX4(k,1,m,0, numEq,2,numComp)];
5954                            const register double B_0_1 = B_p[INDEX4(k,0,m,1, numEq,2,numComp)];
5955                            const register double B_1_1 = B_p[INDEX4(k,1,m,1, numEq,2,numComp)];
5956                            const register double tmp0_0 = B_0_0 + B_0_1;
5957                            const register double tmp9_1 = tmp0_0*w27;
5958                            const register double tmp17_1 = B_1_0*w24;
5959                            const register double tmp10_1 = B_1_0*w26;
5960                            const register double tmp8_1 = B_1_0*w28;
5961                            const register double tmp6_1 = B_1_1*w29;
5962                            const register double tmp13_1 = B_0_1*w23;
5963                            const register double tmp7_1 = tmp0_0*w20;
5964                            const register double tmp0_1 = B_0_1*w25;
5965                            const register double tmp11_1 = B_1_1*w24;
5966                            const register double tmp14_1 = B_0_0*w22;
5967                            const register double tmp16_1 = B_1_1*w26;
5968                            const register double tmp4_1 = B_1_1*w28;
5969                            const register double tmp5_1 = B_1_0*w29;
5970                            const register double tmp3_1 = B_0_0*w21;
5971                            const register double tmp15_1 = B_0_1*w21;
5972                            const register double tmp1_1 = B_0_0*w23;
5973                            const register double tmp2_1 = B_0_1*w22;
5974                            const register double tmp12_1 = B_0_0*w25;
5975                            EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
5976                            EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
5977                            EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1;
5978                            EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp9_1;
5979                            EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp7_1;
5980                            EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp12_1 + tmp13_1;
5981                            EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp9_1;
5982                            EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1;
5983                        }
5984                    }
5985                } else { /* constant data */
5986                    for (index_t k=0; k<numEq; k++) {
5987                        for (index_t m=0; m<numComp; m++) {
5988                            const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
5989                            const register double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
5990                            const register double tmp1_1 = B_1*w35;
5991                            const register double tmp2_1 = B_0*w31;
5992                            const register double tmp5_1 = B_1*w33;
5993                            const register double tmp3_1 = B_0*w30;
5994                            const register double tmp0_1 = B_0*w32;
5995                            const register double tmp4_1 = B_0*w34;
5996                            EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
5997                            EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
5998                            EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
5999                            EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1;
6000                            EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp3_1 + tmp5_1;
6001                            EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
6002                            EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1;
6003                            EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp2_1 + tmp5_1;
6004                        }
6005                    }
6006                }
6007            }
6008            ///////////////
6009            // process C //
6010            ///////////////
6011            if (!C.isEmpty()) {
6012                add_EM_S=true;
6013                const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
6014                if (C.actsExpanded()) {
6015                    for (index_t k=0; k<numEq; k++) {
6016                        for (index_t m=0; m<numComp; m++) {
6017                            const register double C_0_0 = C_p[INDEX4(k,m,0, 0, numEq,numComp,2)];
6018                            const register double C_1_0 = C_p[INDEX4(k,m,1, 0, numEq,numComp,2)];
6019                            const register double C_0_1 = C_p[INDEX4(k,m,0, 1, numEq,numComp,2)];
6020                            const register double C_1_1 = C_p[INDEX4(k,m,1, 1, numEq,numComp,2)];
6021                            const register double tmp0_0 = C_0_0 + C_0_1;
6022                            const register double tmp7_1 = tmp0_0*w20;
6023                            const register double tmp15_1 = C_0_1*w21;
6024                            const register double tmp0_1 = tmp0_0*w27;
6025                            const register double tmp9_1 = C_1_1*w24;
6026                            const register double tmp16_1 = C_1_1*w26;
6027                            const register double tmp5_1 = C_1_1*w28;
6028                            const register double tmp3_1 = C_0_1*w22;
6029                            const register double tmp6_1 = C_1_0*w29;
6030                            const register double tmp11_1 = C_1_0*w28;
6031                            const register double tmp2_1 = C_0_1*w23;
6032                            const register double tmp12_1 = C_0_1*w25;
6033                            const register double tmp17_1 = C_1_0*w24;
6034                            const register double tmp8_1 = C_1_0*w26;
6035                            const register double tmp4_1 = C_0_0*w21;
6036                            const register double tmp10_1 = C_1_1*w29;
6037                            const register double tmp13_1 = C_0_0*w23;
6038                            const register double tmp1_1 = C_0_0*w25;
6039                            const register double tmp14_1 = C_0_0*w22;
6040                            EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0_1;
6041                            EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
6042                            EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1;
6043                            EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0_1;
6044                            EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
6045                            EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp7_1;
6046                            EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp12_1 + tmp13_1;
6047                            EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1;
6048                        }
6049                    }
6050                } else { /* constant data */
6051                    for (index_t k=0; k<numEq; k++) {
6052                        for (index_t m=0; m<numComp; m++) {
6053                            const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
6054                            const register double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
6055                            const register double tmp4_1 = C_0*w30;
6056                            const register double tmp1_1 = C_0*w32;
6057                            const register double tmp0_1 = C_0*w34;
6058                            const register double tmp5_1 = C_1*w33;
6059                            const register double tmp2_1 = C_1*w35;
6060                            const register double tmp3_1 = C_0*w31;
6061                            EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0_1;
6062                            EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp1_1;
6063                            EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
6064                            EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0_1;
6065                            EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1;
6066                            EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp4_1;
6067                            EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1;
6068                            EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;
6069                        }
6070                    }
6071                }
6072            }
6073            ///////////////
6074            // process D //
6075            ///////////////
6076            if (!D.isEmpty()) {
6077                add_EM_S=true;
6078                const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
6079                if (D.actsExpanded()) {
6080                    for (index_t k=0; k<numEq; k++) {
6081                        for (index_t m=0; m<numComp; m++) {
6082                            const register double D_0 = D_p[INDEX3(k, m, 0, numEq, numComp)];
6083                            const register double D_1 = D_p[INDEX3(k, m, 1, numEq, numComp)];
6084                            const register double tmp0_0 = D_0 + D_1;
6085                            const register double tmp4_1 = D_1*w37;
6086                            const register double tmp2_1 = tmp0_0*w38;
6087                            const register double tmp1_1 = D_0*w37;
6088                            const register double tmp3_1 = D_0*w36;
6089                            const register double tmp0_1 = D_1*w36;
6090                            EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
6091                            EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1;
6092                            EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1;
6093                            EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_1;
6094                        }
6095                     }
6096                } else { /* constant data */
6097                    for (index_t k=0; k<numEq; k++) {
6098                        for (index_t m=0; m<numComp; m++) {
6099                            const register double D_0 = D_p[INDEX2(k, m, numEq)];
6100                            const register double tmp0_1 = D_0*w39;
6101                            const register double tmp1_1 = D_0*w40;
6102                            EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
6103                            EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp1_1;
6104                            EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1;
6105                            EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;
6106                        }
6107                    }
6108                }
6109            }
6110            ///////////////
6111            // process X //
6112            ///////////////
6113            if (!X.isEmpty()) {
6114                add_EM_F=true;
6115                const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
6116                if (X.actsExpanded()) {
6117                    for (index_t k=0; k<numEq; k++) {
6118                        const register double X_0_0 = X_p[INDEX3(k, 0, 0, numEq, 2)];
6119                        const register double X_1_0 = X_p[INDEX3(k, 1, 0, numEq, 2)];
6120                        const register double X_0_1 = X_p[INDEX3(k, 0, 1, numEq, 2)];
6121                        const register double X_1_1 = X_p[INDEX3(k, 1, 1, numEq, 2)];
6122                        const register double tmp0_0 = X_1_0 + X_1_1;
6123                        const register double tmp2_1 = X_0_1*w44;
6124                        const register double tmp1_1 = X_0_1*w42;
6125                        const register double tmp7_1 = X_0_0*w44;
6126                        const register double tmp3_1 = X_0_0*w43;
6127                        const register double tmp6_1 = X_0_0*w42;
6128                        const register double tmp8_1 = tmp0_0*w35;
6129                        const register double tmp0_1 = X_0_0*w41;
6130                        const register double tmp9_1 = X_0_1*w43;
6131                        const register double tmp4_1 = tmp0_0*w33;
6132                        const register double tmp5_1 = X_0_1*w41;
6133                        EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
6134                        EM_F[INDEX2(k,1,numEq)]+=tmp2_1 + tmp3_1 + tmp4_1;
6135                        EM_F[INDEX2(k,2,numEq)]+=tmp5_1 + tmp6_1;
6136                        EM_F[INDEX2(k,3,numEq)]+=tmp7_1 + tmp8_1 + tmp9_1;
6137                    }
6138                } else { /* constant data */
6139                    for (index_t k=0; k<numEq; k++) {
6140                        const register double X_0 = X_p[INDEX2(k, 0, numEq)];
6141                        const register double X_1 = X_p[INDEX2(k, 1, numEq)];
6142                        const register double tmp1_1 = X_0*w47;
6143                        const register double tmp0_1 = X_0*w45;
6144                        const register double tmp3_1 = X_1*w48;
6145                        const register double tmp2_1 = X_1*w46;
6146                        EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
6147                        EM_F[INDEX2(k,1,numEq)]+=tmp1_1 + tmp2_1;
6148                        EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
6149                        EM_F[INDEX2(k,3,numEq)]+=tmp1_1 + tmp3_1;
6150                    }
6151                }
6152            }
6153            ///////////////
6154            // process Y //
6155            ///////////////
6156            if (!Y.isEmpty()) {
6157                add_EM_F=true;
6158                const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
6159                if (Y.actsExpanded()) {
6160                    for (index_t k=0; k<numEq; k++) {
6161                        const register double Y_0 = Y_p[INDEX2(k, 0, numEq)];
6162                        const register double Y_1 = Y_p[INDEX2(k, 1, numEq)];
6163                        const register double tmp0_1 = Y_0*w49;
6164                        const register double tmp1_1 = Y_1*w50;
6165                        const register double tmp3_1 = Y_0*w50;
6166                        const register double tmp2_1 = Y_1*w49;
6167                        EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp1_1;
6168                        EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;
6169                    }
6170                } else { /* constant data */
6171                    for (index_t k=0; k<numEq; k++) {
6172                        const register double Y_0 = Y_p[k];
6173                        const register double tmp0_1 = Y_0*w51;
6174                        EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
6175