/[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 4375 by caltinay, Mon Apr 22 05:35:52 2013 UTC revision 4378 by caltinay, Wed Apr 24 06:21:42 2013 UTC
# Line 1808  void Rectangle::assemblePDESingle(Paso_S Line 1808  void Rectangle::assemblePDESingle(Paso_S
1808          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
1809          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
1810  {  {
1811      /* GENERATOR SNIP_PDE_SINGLE_PRE TOP */      const double SQRT3 = 1.73205080756887719318;
1812      const double w0 = -0.1555021169820365539*m_dx[1]/m_dx[0];      const double w1 = 1.0/24.0;
1813      const double w1 = 0.041666666666666666667;      const double w5 = -SQRT3/24 + 1.0/12;
1814      const double w2 = -0.15550211698203655390;      const double w2 = -SQRT3/24 - 1.0/12;
1815      const double w3 = 0.041666666666666666667*m_dx[0]/m_dx[1];      const double w19 = -m_dx[0]/12;
1816      const double w4 = -0.01116454968463011277*m_dx[1]/m_dx[0];      const double w11 = w19*(SQRT3 + 3)/12;
1817      const double w5 = 0.011164549684630112770;      const double w14 = w19*(-SQRT3 + 3)/12;
1818      const double w6 = -0.041666666666666666667*m_dx[1]/m_dx[0];      const double w16 = w19*(5*SQRT3 + 9)/12;
1819      const double w7 = 0.1555021169820365539*m_dx[0]/m_dx[1];      const double w17 = w19*(-5*SQRT3 + 9)/12;
1820      const double w8 = 0.01116454968463011277*m_dx[0]/m_dx[1];      const double w27 = w19*(-SQRT3 - 3)/2;
1821      const double w9 = -1./4;      const double w28 = w19*(SQRT3 - 3)/2;
1822      const double w10 = -0.16666666666666666667*m_dx[0]/m_dx[1];      const double w18 = -m_dx[1]/12;
1823      const double w11 = -0.032861463941450536761*m_dx[1];      const double w12 = w18*(5*SQRT3 + 9)/12;
1824      const double w12 = -0.032861463941450536761*m_dx[0];      const double w13 = w18*(-5*SQRT3 + 9)/12;
1825      const double w13 = -0.12264065304058601714*m_dx[1];      const double w10 = w18*(SQRT3 + 3)/12;
1826      const double w14 = -0.0023593469594139828636*m_dx[1];      const double w15 = w18*(-SQRT3 + 3)/12;
1827      const double w15 = -0.008805202725216129906*m_dx[0];      const double w25 = w18*(-SQRT3 - 3)/2;
1828      const double w16 = -0.008805202725216129906*m_dx[1];      const double w26 = w18*(SQRT3 - 3)/2;
1829      const double w17 = -0.12264065304058601714*m_dx[0];      const double w22 = m_dx[0]*m_dx[1]/144;
1830      const double w18 = -0.0023593469594139828636*m_dx[0];      const double w20 = w22*(SQRT3 + 2);
1831      const double w19 = -0.16666666666666666667*m_dx[1];      const double w21 = w22*(-SQRT3 + 2);
1832      const double w20 = -0.083333333333333333333*m_dx[0];      const double w23 = w22*(4*SQRT3 + 7);
1833      const double w21 = 0.025917019497006092316*m_dx[0]*m_dx[1];      const double w24 = w22*(-4*SQRT3 + 7);
1834      const double w22 = 0.0018607582807716854616*m_dx[0]*m_dx[1];      const double w3 = m_dx[0]/(24*m_dx[1]);
1835      const double w23 = 0.0069444444444444444444*m_dx[0]*m_dx[1];      const double w7 = w3*(SQRT3 + 2);
1836      const double w24 = 0.09672363354357992482*m_dx[0]*m_dx[1];      const double w8 = w3*(-SQRT3 + 2);
1837      const double w25 = 0.000498588678642297402*m_dx[0]*m_dx[1];      const double w6 = -m_dx[1]/(24*m_dx[0]);
1838      const double w26 = 0.19716878364870322056*m_dx[1];      const double w0 = w6*(SQRT3 + 2);
1839      const double w27 = 0.052831216351296779436*m_dx[1];      const double w4 = w6*(-SQRT3 + 2);
     const double w28 = 0.19716878364870322056*m_dx[0];  
     const double w29 = 0.052831216351296779436*m_dx[0];  
     /* GENERATOR SNIP_PDE_SINGLE_PRE BOTTOM */  
1840    
1841      rhs.requireWrite();      rhs.requireWrite();
1842  #pragma omp parallel  #pragma omp parallel
# Line 1853  void Rectangle::assemblePDESingle(Paso_S Line 1850  void Rectangle::assemblePDESingle(Paso_S
1850                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
1851                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
1852                      const index_t e = k0 + m_NE[0]*k1;                      const index_t e = k0 + m_NE[0]*k1;
                     /* GENERATOR SNIP_PDE_SINGLE TOP */  
1853                      ///////////////                      ///////////////
1854                      // process A //                      // process A //
1855                      ///////////////                      ///////////////
# Line 1954  void Rectangle::assemblePDESingle(Paso_S Line 1950  void Rectangle::assemblePDESingle(Paso_S
1950                              const double A_01 = A_p[INDEX2(0,1,2)];                              const double A_01 = A_p[INDEX2(0,1,2)];
1951                              const double A_10 = A_p[INDEX2(1,0,2)];                              const double A_10 = A_p[INDEX2(1,0,2)];
1952                              const double A_11 = A_p[INDEX2(1,1,2)];                              const double A_11 = A_p[INDEX2(1,1,2)];
1953                              const double tmp0 = w9*(-A_01 - A_10);                              const double tmp0 = 6*w1*(A_01 - A_10);
1954                              const double tmp1 = w9*(-A_01 + A_10);                              const double tmp1 = 6*w1*(A_01 + A_10);
1955                              const double tmp2 = w9*(A_01 + A_10);                              const double tmp2 = 6*w1*(-A_01 - A_10);
1956                              const double tmp3 = w9*(A_01 - A_10);                              const double tmp3 = 6*w1*(-A_01 + A_10);
1957                              EM_S[INDEX2(0,0,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;                              EM_S[INDEX2(0,0,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp1;
1958                              EM_S[INDEX2(0,1,4)]+=8*A_00*w6 + 6*A_01*w1 + A_10*w9 + 4*A_11*w3;                              EM_S[INDEX2(0,1,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp0;
1959                              EM_S[INDEX2(0,2,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;                              EM_S[INDEX2(0,2,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp3;
1960                              EM_S[INDEX2(0,3,4)]+=4*A_00*w6 + A_11*w10 + tmp2;                              EM_S[INDEX2(0,3,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp2;
1961                              EM_S[INDEX2(1,0,4)]+=8*A_00*w6 - A_11*w10 + tmp3;                              EM_S[INDEX2(1,0,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp3;
1962                              EM_S[INDEX2(1,1,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;                              EM_S[INDEX2(1,1,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp2;
1963                              EM_S[INDEX2(1,2,4)]+=4*A_00*w6 + A_11*w10 + tmp0;                              EM_S[INDEX2(1,2,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp1;
1964                              EM_S[INDEX2(1,3,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;                              EM_S[INDEX2(1,3,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp0;
1965                              EM_S[INDEX2(2,0,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;                              EM_S[INDEX2(2,0,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp0;
1966                              EM_S[INDEX2(2,1,4)]+=4*A_00*w6 + A_11*w10 + tmp0;                              EM_S[INDEX2(2,1,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp1;
1967                              EM_S[INDEX2(2,2,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;                              EM_S[INDEX2(2,2,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp2;
1968                              EM_S[INDEX2(2,3,4)]+=8*A_00*w6 - A_11*w10 + tmp3;                              EM_S[INDEX2(2,3,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp3;
1969                              EM_S[INDEX2(3,0,4)]+=4*A_00*w6 + A_11*w10 + tmp2;                              EM_S[INDEX2(3,0,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp2;
1970                              EM_S[INDEX2(3,1,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;                              EM_S[INDEX2(3,1,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp3;
1971                              EM_S[INDEX2(3,2,4)]+=8*A_00*w6 - A_11*w10 + tmp1;                              EM_S[INDEX2(3,2,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp0;
1972                              EM_S[INDEX2(3,3,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;                              EM_S[INDEX2(3,3,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp1;
1973                          }                          }
1974                      }                      }
1975                      ///////////////                      ///////////////
# Line 1991  void Rectangle::assemblePDESingle(Paso_S Line 1987  void Rectangle::assemblePDESingle(Paso_S
1987                              const double B_1_2 = B_p[INDEX2(1,2,2)];                              const double B_1_2 = B_p[INDEX2(1,2,2)];
1988                              const double B_0_3 = B_p[INDEX2(0,3,2)];                              const double B_0_3 = B_p[INDEX2(0,3,2)];
1989                              const double B_1_3 = B_p[INDEX2(1,3,2)];                              const double B_1_3 = B_p[INDEX2(1,3,2)];
1990                              const double tmp0 = w15*(B_1_2 + B_1_3);                              const double tmp0 = w11*(B_1_0 + B_1_1);
1991                              const double tmp1 = w12*(B_1_0 + B_1_1);                              const double tmp1 = w14*(B_1_2 + B_1_3);
1992                              const double tmp2 = w15*(B_1_0 + B_1_1);                              const double tmp2 = w15*(-B_0_1 - B_0_3);
1993                              const double tmp3 = w16*(-B_0_1 - B_0_3);                              const double tmp3 = w10*(-B_0_0 - B_0_2);
1994                              const double tmp4 = w11*(-B_0_0 - B_0_2);                              const double tmp4 = w11*(B_1_2 + B_1_3);
1995                              const double tmp5 = w12*(B_1_2 + B_1_3);                              const double tmp5 = w14*(B_1_0 + B_1_1);
1996                              const double tmp6 = w15*(-B_1_0 - B_1_1);                              const double tmp6 = w11*(-B_1_2 - B_1_3);
1997                              const double tmp7 = w12*(-B_1_2 - B_1_3);                              const double tmp7 = w14*(-B_1_0 - B_1_1);
1998                              const double tmp8 = w15*(-B_1_2 - B_1_3);                              const double tmp8 = w11*(-B_1_0 - B_1_1);
1999                              const double tmp9 = w12*(-B_1_0 - B_1_1);                              const double tmp9 = w14*(-B_1_2 - B_1_3);
2000                              const double tmp10 = w11*(-B_0_1 - B_0_3);                              const double tmp10 = w10*(-B_0_1 - B_0_3);
2001                              const double tmp11 = w16*(-B_0_0 - B_0_2);                              const double tmp11 = w15*(-B_0_0 - B_0_2);
2002                              const double tmp12 = w16*(B_0_0 + B_0_2);                              const double tmp12 = w15*(B_0_0 + B_0_2);
2003                              const double tmp13 = w11*(B_0_1 + B_0_3);                              const double tmp13 = w10*(B_0_1 + B_0_3);
2004                              const double tmp14 = w11*(B_0_0 + B_0_2);                              const double tmp14 = w10*(B_0_0 + B_0_2);
2005                              const double tmp15 = w16*(B_0_1 + B_0_3);                              const double tmp15 = w15*(B_0_1 + B_0_3);
2006                              EM_S[INDEX2(0,0,4)]+=B_0_0*w13 + B_0_1*w11 + B_0_2*w16 + B_0_3*w14 + B_1_0*w17 + B_1_1*w15 + B_1_2*w12 + B_1_3*w18;                              EM_S[INDEX2(0,0,4)]+=B_0_0*w12 + B_0_1*w10 + B_0_2*w15 + B_0_3*w13 + B_1_0*w16 + B_1_1*w14 + B_1_2*w11 + B_1_3*w17;
2007                              EM_S[INDEX2(0,1,4)]+=B_0_0*w11 + B_0_1*w13 + B_0_2*w14 + B_0_3*w16 + tmp0 + tmp1;                              EM_S[INDEX2(0,1,4)]+=B_0_0*w10 + B_0_1*w12 + B_0_2*w13 + B_0_3*w15 + tmp0 + tmp1;
2008                              EM_S[INDEX2(0,2,4)]+=B_1_0*w12 + B_1_1*w18 + B_1_2*w17 + B_1_3*w15 + tmp14 + tmp15;                              EM_S[INDEX2(0,2,4)]+=B_1_0*w11 + B_1_1*w17 + B_1_2*w16 + B_1_3*w14 + tmp14 + tmp15;
2009                              EM_S[INDEX2(0,3,4)]+=tmp12 + tmp13 + tmp2 + tmp5;                              EM_S[INDEX2(0,3,4)]+=tmp12 + tmp13 + tmp4 + tmp5;
2010                              EM_S[INDEX2(1,0,4)]+=-B_0_0*w13 - B_0_1*w11 - B_0_2*w16 - B_0_3*w14 + tmp0 + tmp1;                              EM_S[INDEX2(1,0,4)]+=-B_0_0*w12 - B_0_1*w10 - B_0_2*w15 - B_0_3*w13 + tmp0 + tmp1;
2011                              EM_S[INDEX2(1,1,4)]+=-B_0_0*w11 - B_0_1*w13 - B_0_2*w14 - B_0_3*w16 + B_1_0*w15 + B_1_1*w17 + B_1_2*w18 + B_1_3*w12;                              EM_S[INDEX2(1,1,4)]+=-B_0_0*w10 - B_0_1*w12 - B_0_2*w13 - B_0_3*w15 + B_1_0*w14 + B_1_1*w16 + B_1_2*w17 + B_1_3*w11;
2012                              EM_S[INDEX2(1,2,4)]+=tmp2 + tmp3 + tmp4 + tmp5;                              EM_S[INDEX2(1,2,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2013                              EM_S[INDEX2(1,3,4)]+=B_1_0*w18 + B_1_1*w12 + B_1_2*w15 + B_1_3*w17 + tmp10 + tmp11;                              EM_S[INDEX2(1,3,4)]+=B_1_0*w17 + B_1_1*w11 + B_1_2*w14 + B_1_3*w16 + tmp10 + tmp11;
2014                              EM_S[INDEX2(2,0,4)]+=-B_1_0*w17 - B_1_1*w15 - B_1_2*w12 - B_1_3*w18 + tmp14 + tmp15;                              EM_S[INDEX2(2,0,4)]+=-B_1_0*w16 - B_1_1*w14 - B_1_2*w11 - B_1_3*w17 + tmp14 + tmp15;
2015                              EM_S[INDEX2(2,1,4)]+=tmp12 + tmp13 + tmp8 + tmp9;                              EM_S[INDEX2(2,1,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2016                              EM_S[INDEX2(2,2,4)]+=B_0_0*w16 + B_0_1*w14 + B_0_2*w13 + B_0_3*w11 - B_1_0*w12 - B_1_1*w18 - B_1_2*w17 - B_1_3*w15;                              EM_S[INDEX2(2,2,4)]+=B_0_0*w15 + B_0_1*w13 + B_0_2*w12 + B_0_3*w10 - B_1_0*w11 - B_1_1*w17 - B_1_2*w16 - B_1_3*w14;
2017                              EM_S[INDEX2(2,3,4)]+=B_0_0*w14 + B_0_1*w16 + B_0_2*w11 + B_0_3*w13 + tmp6 + tmp7;                              EM_S[INDEX2(2,3,4)]+=B_0_0*w13 + B_0_1*w15 + B_0_2*w10 + B_0_3*w12 + tmp6 + tmp7;
2018                              EM_S[INDEX2(3,0,4)]+=tmp3 + tmp4 + tmp8 + tmp9;                              EM_S[INDEX2(3,0,4)]+=tmp2 + tmp3 + tmp8 + tmp9;
2019                              EM_S[INDEX2(3,1,4)]+=-B_1_0*w15 - B_1_1*w17 - B_1_2*w18 - B_1_3*w12 + tmp10 + tmp11;                              EM_S[INDEX2(3,1,4)]+=-B_1_0*w14 - B_1_1*w16 - B_1_2*w17 - B_1_3*w11 + tmp10 + tmp11;
2020                              EM_S[INDEX2(3,2,4)]+=-B_0_0*w16 - B_0_1*w14 - B_0_2*w13 - B_0_3*w11 + tmp6 + tmp7;                              EM_S[INDEX2(3,2,4)]+=-B_0_0*w15 - B_0_1*w13 - B_0_2*w12 - B_0_3*w10 + tmp6 + tmp7;
2021                              EM_S[INDEX2(3,3,4)]+=-B_0_0*w14 - B_0_1*w16 - B_0_2*w11 - B_0_3*w13 - B_1_0*w18 - B_1_1*w12 - B_1_2*w15 - B_1_3*w17;                              EM_S[INDEX2(3,3,4)]+=-B_0_0*w13 - B_0_1*w15 - B_0_2*w10 - B_0_3*w12 - B_1_0*w17 - B_1_1*w11 - B_1_2*w14 - B_1_3*w16;
2022                          } else { // constant data                          } else { // constant data
2023                              const double B_0 = B_p[0];                              const double B_0 = B_p[0];
2024                              const double B_1 = B_p[1];                              const double B_1 = B_p[1];
2025                              EM_S[INDEX2(0,0,4)]+=B_0*w19 + 2*B_1*w20;                              EM_S[INDEX2(0,0,4)]+= 2*B_0*w18 + 2*B_1*w19;
2026                              EM_S[INDEX2(0,1,4)]+=B_0*w19 + B_1*w20;                              EM_S[INDEX2(0,1,4)]+= 2*B_0*w18 +   B_1*w19;
2027                              EM_S[INDEX2(0,2,4)]+=B_0*w19/2 + 2*B_1*w20;                              EM_S[INDEX2(0,2,4)]+=   B_0*w18 + 2*B_1*w19;
2028                              EM_S[INDEX2(0,3,4)]+=B_0*w19/2 + B_1*w20;                              EM_S[INDEX2(0,3,4)]+=   B_0*w18 +   B_1*w19;
2029                              EM_S[INDEX2(1,0,4)]+=-B_0*w19 + B_1*w20;                              EM_S[INDEX2(1,0,4)]+=-2*B_0*w18 +   B_1*w19;
2030                              EM_S[INDEX2(1,1,4)]+=-B_0*w19 + 2*B_1*w20;                              EM_S[INDEX2(1,1,4)]+=-2*B_0*w18 + 2*B_1*w19;
2031                              EM_S[INDEX2(1,2,4)]+=-B_0*w19/2 + B_1*w20;                              EM_S[INDEX2(1,2,4)]+=  -B_0*w18 +   B_1*w19;
2032                              EM_S[INDEX2(1,3,4)]+=-B_0*w19/2 + 2*B_1*w20;                              EM_S[INDEX2(1,3,4)]+=  -B_0*w18 + 2*B_1*w19;
2033                              EM_S[INDEX2(2,0,4)]+=B_0*w19/2 - 2*B_1*w20;                              EM_S[INDEX2(2,0,4)]+=   B_0*w18 - 2*B_1*w19;
2034                              EM_S[INDEX2(2,1,4)]+=B_0*w19/2 - B_1*w20;                              EM_S[INDEX2(2,1,4)]+=   B_0*w18 -   B_1*w19;
2035                              EM_S[INDEX2(2,2,4)]+=B_0*w19 - 2*B_1*w20;                              EM_S[INDEX2(2,2,4)]+= 2*B_0*w18 - 2*B_1*w19;
2036                              EM_S[INDEX2(2,3,4)]+=B_0*w19 - B_1*w20;                              EM_S[INDEX2(2,3,4)]+= 2*B_0*w18 -   B_1*w19;
2037                              EM_S[INDEX2(3,0,4)]+=-B_0*w19/2 - B_1*w20;                              EM_S[INDEX2(3,0,4)]+=  -B_0*w18 -   B_1*w19;
2038                              EM_S[INDEX2(3,1,4)]+=-B_0*w19/2 - 2*B_1*w20;                              EM_S[INDEX2(3,1,4)]+=  -B_0*w18 - 2*B_1*w19;
2039                              EM_S[INDEX2(3,2,4)]+=-B_0*w19 - B_1*w20;                              EM_S[INDEX2(3,2,4)]+=-2*B_0*w18 -   B_1*w19;
2040                              EM_S[INDEX2(3,3,4)]+=-B_0*w19 - 2*B_1*w20;                              EM_S[INDEX2(3,3,4)]+=-2*B_0*w18 - 2*B_1*w19;
2041                          }                          }
2042                      }                      }
2043                      ///////////////                      ///////////////
# Line 2059  void Rectangle::assemblePDESingle(Paso_S Line 2055  void Rectangle::assemblePDESingle(Paso_S
2055                              const double C_1_2 = C_p[INDEX2(1,2,2)];                              const double C_1_2 = C_p[INDEX2(1,2,2)];
2056                              const double C_0_3 = C_p[INDEX2(0,3,2)];                              const double C_0_3 = C_p[INDEX2(0,3,2)];
2057                              const double C_1_3 = C_p[INDEX2(1,3,2)];                              const double C_1_3 = C_p[INDEX2(1,3,2)];
2058                              const double tmp0 = w15*(C_1_2 + C_1_3);                              const double tmp0 = w11*(C_1_0 + C_1_1);
2059                              const double tmp1 = w12*(C_1_0 + C_1_1);                              const double tmp1 = w14*(C_1_2 + C_1_3);
2060                              const double tmp2 = w15*(-C_1_2 - C_1_3);                              const double tmp2 = w15*(C_0_0 + C_0_2);
2061                              const double tmp3 = w16*(C_0_0 + C_0_2);                              const double tmp3 = w10*(C_0_1 + C_0_3);
2062                              const double tmp4 = w11*(C_0_1 + C_0_3);                              const double tmp4 = w11*(-C_1_0 - C_1_1);
2063                              const double tmp5 = w12*(-C_1_0 - C_1_1);                              const double tmp5 = w14*(-C_1_2 - C_1_3);
2064                              const double tmp6 = w15*(-C_1_0 - C_1_1);                              const double tmp6 = w11*(-C_1_2 - C_1_3);
2065                              const double tmp7 = w12*(-C_1_2 - C_1_3);                              const double tmp7 = w14*(-C_1_0 - C_1_1);
2066                              const double tmp8 = w15*(C_1_0 + C_1_1);                              const double tmp8 = w11*(C_1_2 + C_1_3);
2067                              const double tmp9 = w12*(C_1_2 + C_1_3);                              const double tmp9 = w14*(C_1_0 + C_1_1);
2068                              const double tmp10 = w11*(-C_0_1 - C_0_3);                              const double tmp10 = w10*(-C_0_1 - C_0_3);
2069                              const double tmp11 = w16*(-C_0_0 - C_0_2);                              const double tmp11 = w15*(-C_0_0 - C_0_2);
2070                              const double tmp12 = w16*(-C_0_1 - C_0_3);                              const double tmp12 = w15*(-C_0_1 - C_0_3);
2071                              const double tmp13 = w11*(-C_0_0 - C_0_2);                              const double tmp13 = w10*(-C_0_0 - C_0_2);
2072                              const double tmp14 = w11*(C_0_0 + C_0_2);                              const double tmp14 = w10*(C_0_0 + C_0_2);
2073                              const double tmp15 = w16*(C_0_1 + C_0_3);                              const double tmp15 = w15*(C_0_1 + C_0_3);
2074                              EM_S[INDEX2(0,0,4)]+=C_0_0*w13 + C_0_1*w11 + C_0_2*w16 + C_0_3*w14 + C_1_0*w17 + C_1_1*w15 + C_1_2*w12 + C_1_3*w18;                              EM_S[INDEX2(0,0,4)]+=C_0_0*w12 + C_0_1*w10 + C_0_2*w15 + C_0_3*w13 + C_1_0*w16 + C_1_1*w14 + C_1_2*w11 + C_1_3*w17;
2075                              EM_S[INDEX2(0,1,4)]+=-C_0_0*w13 - C_0_1*w11 - C_0_2*w16 - C_0_3*w14 + tmp0 + tmp1;                              EM_S[INDEX2(0,1,4)]+=-C_0_0*w12 - C_0_1*w10 - C_0_2*w15 - C_0_3*w13 + tmp0 + tmp1;
2076                              EM_S[INDEX2(0,2,4)]+=-C_1_0*w17 - C_1_1*w15 - C_1_2*w12 - C_1_3*w18 + tmp14 + tmp15;                              EM_S[INDEX2(0,2,4)]+=-C_1_0*w16 - C_1_1*w14 - C_1_2*w11 - C_1_3*w17 + tmp14 + tmp15;
2077                              EM_S[INDEX2(0,3,4)]+=tmp12 + tmp13 + tmp2 + tmp5;                              EM_S[INDEX2(0,3,4)]+=tmp12 + tmp13 + tmp4 + tmp5;
2078                              EM_S[INDEX2(1,0,4)]+=C_0_0*w11 + C_0_1*w13 + C_0_2*w14 + C_0_3*w16 + tmp0 + tmp1;                              EM_S[INDEX2(1,0,4)]+=C_0_0*w10 + C_0_1*w12 + C_0_2*w13 + C_0_3*w15 + tmp0 + tmp1;
2079                              EM_S[INDEX2(1,1,4)]+=-C_0_0*w11 - C_0_1*w13 - C_0_2*w14 - C_0_3*w16 + C_1_0*w15 + C_1_1*w17 + C_1_2*w18 + C_1_3*w12;                              EM_S[INDEX2(1,1,4)]+=-C_0_0*w10 - C_0_1*w12 - C_0_2*w13 - C_0_3*w15 + C_1_0*w14 + C_1_1*w16 + C_1_2*w17 + C_1_3*w11;
2080                              EM_S[INDEX2(1,2,4)]+=tmp2 + tmp3 + tmp4 + tmp5;                              EM_S[INDEX2(1,2,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2081                              EM_S[INDEX2(1,3,4)]+=-C_1_0*w15 - C_1_1*w17 - C_1_2*w18 - C_1_3*w12 + tmp10 + tmp11;                              EM_S[INDEX2(1,3,4)]+=-C_1_0*w14 - C_1_1*w16 - C_1_2*w17 - C_1_3*w11 + tmp10 + tmp11;
2082                              EM_S[INDEX2(2,0,4)]+=C_1_0*w12 + C_1_1*w18 + C_1_2*w17 + C_1_3*w15 + tmp14 + tmp15;                              EM_S[INDEX2(2,0,4)]+=C_1_0*w11 + C_1_1*w17 + C_1_2*w16 + C_1_3*w14 + tmp14 + tmp15;
2083                              EM_S[INDEX2(2,1,4)]+=tmp12 + tmp13 + tmp8 + tmp9;                              EM_S[INDEX2(2,1,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2084                              EM_S[INDEX2(2,2,4)]+=C_0_0*w16 + C_0_1*w14 + C_0_2*w13 + C_0_3*w11 - C_1_0*w12 - C_1_1*w18 - C_1_2*w17 - C_1_3*w15;                              EM_S[INDEX2(2,2,4)]+=C_0_0*w15 + C_0_1*w13 + C_0_2*w12 + C_0_3*w10 - C_1_0*w11 - C_1_1*w17 - C_1_2*w16 - C_1_3*w14;
2085                              EM_S[INDEX2(2,3,4)]+=-C_0_0*w16 - C_0_1*w14 - C_0_2*w13 - C_0_3*w11 + tmp6 + tmp7;                              EM_S[INDEX2(2,3,4)]+=-C_0_0*w15 - C_0_1*w13 - C_0_2*w12 - C_0_3*w10 + tmp6 + tmp7;
2086                              EM_S[INDEX2(3,0,4)]+=tmp3 + tmp4 + tmp8 + tmp9;                              EM_S[INDEX2(3,0,4)]+=tmp2 + tmp3 + tmp8 + tmp9;
2087                              EM_S[INDEX2(3,1,4)]+=C_1_0*w18 + C_1_1*w12 + C_1_2*w15 + C_1_3*w17 + tmp10 + tmp11;                              EM_S[INDEX2(3,1,4)]+=C_1_0*w17 + C_1_1*w11 + C_1_2*w14 + C_1_3*w16 + tmp10 + tmp11;
2088                              EM_S[INDEX2(3,2,4)]+=C_0_0*w14 + C_0_1*w16 + C_0_2*w11 + C_0_3*w13 + tmp6 + tmp7;                              EM_S[INDEX2(3,2,4)]+=C_0_0*w13 + C_0_1*w15 + C_0_2*w10 + C_0_3*w12 + tmp6 + tmp7;
2089                              EM_S[INDEX2(3,3,4)]+=-C_0_0*w14 - C_0_1*w16 - C_0_2*w11 - C_0_3*w13 - C_1_0*w18 - C_1_1*w12 - C_1_2*w15 - C_1_3*w17;                              EM_S[INDEX2(3,3,4)]+=-C_0_0*w13 - C_0_1*w15 - C_0_2*w10 - C_0_3*w12 - C_1_0*w17 - C_1_1*w11 - C_1_2*w14 - C_1_3*w16;
2090                          } else { // constant data                          } else { // constant data
2091                              const double C_0 = C_p[0];                              const double C_0 = C_p[0];
2092                              const double C_1 = C_p[1];                              const double C_1 = C_p[1];
2093                              EM_S[INDEX2(0,0,4)]+=C_0*w19 + 2*C_1*w20;                              EM_S[INDEX2(0,0,4)]+= 2*C_0*w18 + 2*C_1*w19;
2094                              EM_S[INDEX2(0,1,4)]+=-C_0*w19 + C_1*w20;                              EM_S[INDEX2(0,1,4)]+=-2*C_0*w18 +   C_1*w19;
2095                              EM_S[INDEX2(0,2,4)]+=C_0*w19/2 - 2*C_1*w20;                              EM_S[INDEX2(0,2,4)]+=   C_0*w18 - 2*C_1*w19;
2096                              EM_S[INDEX2(0,3,4)]+=-C_0*w19/2 - C_1*w20;                              EM_S[INDEX2(0,3,4)]+=  -C_0*w18 -   C_1*w19;
2097                              EM_S[INDEX2(1,0,4)]+=C_0*w19 + C_1*w20;                              EM_S[INDEX2(1,0,4)]+= 2*C_0*w18 +   C_1*w19;
2098                              EM_S[INDEX2(1,1,4)]+=-C_0*w19 + 2*C_1*w20;                              EM_S[INDEX2(1,1,4)]+=-2*C_0*w18 + 2*C_1*w19;
2099                              EM_S[INDEX2(1,2,4)]+=C_0*w19/2 - C_1*w20;                              EM_S[INDEX2(1,2,4)]+=   C_0*w18 -   C_1*w19;
2100                              EM_S[INDEX2(1,3,4)]+=-C_0*w19/2 - 2*C_1*w20;                              EM_S[INDEX2(1,3,4)]+=  -C_0*w18 - 2*C_1*w19;
2101                              EM_S[INDEX2(2,0,4)]+=C_0*w19/2 + 2*C_1*w20;                              EM_S[INDEX2(2,0,4)]+=   C_0*w18 + 2*C_1*w19;
2102                              EM_S[INDEX2(2,1,4)]+=-C_0*w19/2 + C_1*w20;                              EM_S[INDEX2(2,1,4)]+=  -C_0*w18 +   C_1*w19;
2103                              EM_S[INDEX2(2,2,4)]+=C_0*w19 - 2*C_1*w20;                              EM_S[INDEX2(2,2,4)]+= 2*C_0*w18 - 2*C_1*w19;
2104                              EM_S[INDEX2(2,3,4)]+=-C_0*w19 - C_1*w20;                              EM_S[INDEX2(2,3,4)]+=-2*C_0*w18 -   C_1*w19;
2105                              EM_S[INDEX2(3,0,4)]+=C_0*w19/2 + C_1*w20;                              EM_S[INDEX2(3,0,4)]+=   C_0*w18 +   C_1*w19;
2106                              EM_S[INDEX2(3,1,4)]+=-C_0*w19/2 + 2*C_1*w20;                              EM_S[INDEX2(3,1,4)]+=  -C_0*w18 + 2*C_1*w19;
2107                              EM_S[INDEX2(3,2,4)]+=C_0*w19 - C_1*w20;                              EM_S[INDEX2(3,2,4)]+= 2*C_0*w18 -   C_1*w19;
2108                              EM_S[INDEX2(3,3,4)]+=-C_0*w19 - 2*C_1*w20;                              EM_S[INDEX2(3,3,4)]+=-2*C_0*w18 - 2*C_1*w19;
2109                          }                          }
2110                      }                      }
2111                      ///////////////                      ///////////////
# Line 2123  void Rectangle::assemblePDESingle(Paso_S Line 2119  void Rectangle::assemblePDESingle(Paso_S
2119                              const double D_1 = D_p[1];                              const double D_1 = D_p[1];
2120                              const double D_2 = D_p[2];                              const double D_2 = D_p[2];
2121                              const double D_3 = D_p[3];                              const double D_3 = D_p[3];
2122                              const double tmp0 = w21*(D_0 + D_1);                              const double tmp0 = w21*(D_2 + D_3);
2123                              const double tmp1 = w22*(D_2 + D_3);                              const double tmp1 = w20*(D_0 + D_1);
2124                              const double tmp2 = w23*(D_0 + D_1 + D_2 + D_3);                              const double tmp2 = w22*(D_0 + D_1 + D_2 + D_3);
2125                              const double tmp3 = w21*(D_2 + D_3);                              const double tmp3 = w21*(D_0 + D_1);
2126                              const double tmp4 = w22*(D_0 + D_1);                              const double tmp4 = w20*(D_2 + D_3);
2127                              const double tmp5 = w23*(D_1 + D_2);                              const double tmp5 = w22*(D_1 + D_2);
2128                              const double tmp6 = w21*(D_1 + D_3);                              const double tmp6 = w21*(D_0 + D_2);
2129                              const double tmp7 = w22*(D_0 + D_2);                              const double tmp7 = w20*(D_1 + D_3);
2130                              const double tmp8 = w21*(D_0 + D_2);                              const double tmp8 = w21*(D_1 + D_3);
2131                              const double tmp9 = w22*(D_1 + D_3);                              const double tmp9 = w20*(D_0 + D_2);
2132                              const double tmp10 = w23*(D_0 + D_3);                              const double tmp10 = w22*(D_0 + D_3);
2133                              EM_S[INDEX2(0,0,4)]+=D_0*w24 + D_3*w25 + tmp5;                              EM_S[INDEX2(0,0,4)]+=D_0*w23 + D_3*w24 + tmp5;
2134                              EM_S[INDEX2(0,1,4)]+=tmp0 + tmp1;                              EM_S[INDEX2(0,1,4)]+=tmp0 + tmp1;
2135                              EM_S[INDEX2(0,2,4)]+=tmp8 + tmp9;                              EM_S[INDEX2(0,2,4)]+=tmp8 + tmp9;
2136                              EM_S[INDEX2(0,3,4)]+=tmp2;                              EM_S[INDEX2(0,3,4)]+=tmp2;
2137                              EM_S[INDEX2(1,0,4)]+=tmp0 + tmp1;                              EM_S[INDEX2(1,0,4)]+=tmp0 + tmp1;
2138                              EM_S[INDEX2(1,1,4)]+=D_1*w24 + D_2*w25 + tmp10;                              EM_S[INDEX2(1,1,4)]+=D_1*w23 + D_2*w24 + tmp10;
2139                              EM_S[INDEX2(1,2,4)]+=tmp2;                              EM_S[INDEX2(1,2,4)]+=tmp2;
2140                              EM_S[INDEX2(1,3,4)]+=tmp6 + tmp7;                              EM_S[INDEX2(1,3,4)]+=tmp6 + tmp7;
2141                              EM_S[INDEX2(2,0,4)]+=tmp8 + tmp9;                              EM_S[INDEX2(2,0,4)]+=tmp8 + tmp9;
2142                              EM_S[INDEX2(2,1,4)]+=tmp2;                              EM_S[INDEX2(2,1,4)]+=tmp2;
2143                              EM_S[INDEX2(2,2,4)]+=D_1*w25 + D_2*w24 + tmp10;                              EM_S[INDEX2(2,2,4)]+=D_1*w24 + D_2*w23 + tmp10;
2144                              EM_S[INDEX2(2,3,4)]+=tmp3 + tmp4;                              EM_S[INDEX2(2,3,4)]+=tmp3 + tmp4;
2145                              EM_S[INDEX2(3,0,4)]+=tmp2;                              EM_S[INDEX2(3,0,4)]+=tmp2;
2146                              EM_S[INDEX2(3,1,4)]+=tmp6 + tmp7;                              EM_S[INDEX2(3,1,4)]+=tmp6 + tmp7;
2147                              EM_S[INDEX2(3,2,4)]+=tmp3 + tmp4;                              EM_S[INDEX2(3,2,4)]+=tmp3 + tmp4;
2148                              EM_S[INDEX2(3,3,4)]+=D_0*w25 + D_3*w24 + tmp5;                              EM_S[INDEX2(3,3,4)]+=D_0*w24 + D_3*w23 + tmp5;
2149                          } else { // constant data                          } else { // constant data
2150                              const double D_0 = D_p[0];                              const double D_0 = D_p[0];
2151                              EM_S[INDEX2(0,0,4)]+=16*D_0*w23;                              EM_S[INDEX2(0,0,4)]+=16*D_0*w22;
2152                              EM_S[INDEX2(0,1,4)]+=8*D_0*w23;                              EM_S[INDEX2(0,1,4)]+=8*D_0*w22;
2153                              EM_S[INDEX2(0,2,4)]+=8*D_0*w23;                              EM_S[INDEX2(0,2,4)]+=8*D_0*w22;
2154                              EM_S[INDEX2(0,3,4)]+=4*D_0*w23;                              EM_S[INDEX2(0,3,4)]+=4*D_0*w22;
2155                              EM_S[INDEX2(1,0,4)]+=8*D_0*w23;                              EM_S[INDEX2(1,0,4)]+=8*D_0*w22;
2156                              EM_S[INDEX2(1,1,4)]+=16*D_0*w23;                              EM_S[INDEX2(1,1,4)]+=16*D_0*w22;
2157                              EM_S[INDEX2(1,2,4)]+=4*D_0*w23;                              EM_S[INDEX2(1,2,4)]+=4*D_0*w22;
2158                              EM_S[INDEX2(1,3,4)]+=8*D_0*w23;                              EM_S[INDEX2(1,3,4)]+=8*D_0*w22;
2159                              EM_S[INDEX2(2,0,4)]+=8*D_0*w23;                              EM_S[INDEX2(2,0,4)]+=8*D_0*w22;
2160                              EM_S[INDEX2(2,1,4)]+=4*D_0*w23;                              EM_S[INDEX2(2,1,4)]+=4*D_0*w22;
2161                              EM_S[INDEX2(2,2,4)]+=16*D_0*w23;                              EM_S[INDEX2(2,2,4)]+=16*D_0*w22;
2162                              EM_S[INDEX2(2,3,4)]+=8*D_0*w23;                              EM_S[INDEX2(2,3,4)]+=8*D_0*w22;
2163                              EM_S[INDEX2(3,0,4)]+=4*D_0*w23;                              EM_S[INDEX2(3,0,4)]+=4*D_0*w22;
2164                              EM_S[INDEX2(3,1,4)]+=8*D_0*w23;                              EM_S[INDEX2(3,1,4)]+=8*D_0*w22;
2165                              EM_S[INDEX2(3,2,4)]+=8*D_0*w23;                              EM_S[INDEX2(3,2,4)]+=8*D_0*w22;
2166                              EM_S[INDEX2(3,3,4)]+=16*D_0*w23;                              EM_S[INDEX2(3,3,4)]+=16*D_0*w22;
2167                          }                          }
2168                      }                      }
2169                      ///////////////                      ///////////////
# Line 2185  void Rectangle::assemblePDESingle(Paso_S Line 2181  void Rectangle::assemblePDESingle(Paso_S
2181                              const double X_1_2 = X_p[INDEX2(1,2,2)];                              const double X_1_2 = X_p[INDEX2(1,2,2)];
2182                              const double X_0_3 = X_p[INDEX2(0,3,2)];                              const double X_0_3 = X_p[INDEX2(0,3,2)];
2183                              const double X_1_3 = X_p[INDEX2(1,3,2)];                              const double X_1_3 = X_p[INDEX2(1,3,2)];
2184                              const double tmp0 = 6*w15*(X_1_1 + X_1_3);                              const double tmp0 = 6*w15*(X_0_2 + X_0_3);
2185                              const double tmp1 = 6*w16*(X_0_2 + X_0_3);                              const double tmp1 = 6*w10*(X_0_0 + X_0_1);
2186                              const double tmp2 = 6*w11*(X_0_0 + X_0_1);                              const double tmp2 = 6*w11*(X_1_0 + X_1_2);
2187                              const double tmp3 = 6*w12*(X_1_0 + X_1_2);                              const double tmp3 = 6*w14*(X_1_1 + X_1_3);
2188                              const double tmp4 = 6*w15*(X_1_0 + X_1_2);                              const double tmp4 = 6*w11*(X_1_1 + X_1_3);
2189                              const double tmp5 = w27*(X_0_2 + X_0_3);                              const double tmp5 = w25*(X_0_0 + X_0_1);
2190                              const double tmp6 = w26*(X_0_0 + X_0_1);                              const double tmp6 = w26*(X_0_2 + X_0_3);
2191                              const double tmp7 = 6*w12*(X_1_1 + X_1_3);                              const double tmp7 = 6*w14*(X_1_0 + X_1_2);
2192                              const double tmp8 = w29*(X_1_1 + X_1_3);                              const double tmp8 = w27*(X_1_0 + X_1_2);
2193                              const double tmp9 = w27*(-X_0_0 - X_0_1);                              const double tmp9 = w28*(X_1_1 + X_1_3);
2194                              const double tmp10 = w28*(X_1_0 + X_1_2);                              const double tmp10 = w25*(-X_0_2 - X_0_3);
2195                              const double tmp11 = w26*(-X_0_2 - X_0_3);                              const double tmp11 = w26*(-X_0_0 - X_0_1);
2196                              const double tmp12 = w29*(X_1_0 + X_1_2);                              const double tmp12 = w27*(X_1_1 + X_1_3);
2197                              const double tmp13 = w27*(X_0_0 + X_0_1);                              const double tmp13 = w28*(X_1_0 + X_1_2);
2198                              const double tmp14 = w28*(X_1_1 + X_1_3);                              const double tmp14 = w25*(X_0_2 + X_0_3);
2199                              const double tmp15 = w26*(X_0_2 + X_0_3);                              const double tmp15 = w26*(X_0_0 + X_0_1);
2200                              EM_F[0]+=tmp0 + tmp1 + tmp2 + tmp3;                              EM_F[0]+=tmp0 + tmp1 + tmp2 + tmp3;
2201                              EM_F[1]+=tmp4 + tmp5 + tmp6 + tmp7;                              EM_F[1]+=tmp4 + tmp5 + tmp6 + tmp7;
2202                              EM_F[2]+=tmp10 + tmp11 + tmp8 + tmp9;                              EM_F[2]+=tmp10 + tmp11 + tmp8 + tmp9;
# Line 2208  void Rectangle::assemblePDESingle(Paso_S Line 2204  void Rectangle::assemblePDESingle(Paso_S
2204                          } else { // constant data                          } else { // constant data
2205                              const double X_0 = X_p[0];                              const double X_0 = X_p[0];
2206                              const double X_1 = X_p[1];                              const double X_1 = X_p[1];
2207                              EM_F[0]+=3*X_0*w19 + 6*X_1*w20;                              EM_F[0]+= 6*X_0*w18 + 6*X_1*w19;
2208                              EM_F[1]+=-3*X_0*w19 + 6*X_1*w20;                              EM_F[1]+=-6*X_0*w18 + 6*X_1*w19;
2209                              EM_F[2]+=3*X_0*w19 - 6*X_1*w20;                              EM_F[2]+= 6*X_0*w18 - 6*X_1*w19;
2210                              EM_F[3]+=-3*X_0*w19 - 6*X_1*w20;                              EM_F[3]+=-6*X_0*w18 - 6*X_1*w19;
2211                          }                          }
2212                      }                      }
2213                      ///////////////                      ///////////////
# Line 2225  void Rectangle::assemblePDESingle(Paso_S Line 2221  void Rectangle::assemblePDESingle(Paso_S
2221                              const double Y_1 = Y_p[1];                              const double Y_1 = Y_p[1];
2222                              const double Y_2 = Y_p[2];                              const double Y_2 = Y_p[2];
2223                              const double Y_3 = Y_p[3];                              const double Y_3 = Y_p[3];
2224                              const double tmp0 = 6*w23*(Y_1 + Y_2);                              const double tmp0 = 6*w22*(Y_1 + Y_2);
2225                              const double tmp1 = 6*w23*(Y_0 + Y_3);                              const double tmp1 = 6*w22*(Y_0 + Y_3);
2226                              EM_F[0]+=6*Y_0*w21 + 6*Y_3*w22 + tmp0;                              EM_F[0]+=6*Y_0*w20 + 6*Y_3*w21 + tmp0;
2227                              EM_F[1]+=6*Y_1*w21 + 6*Y_2*w22 + tmp1;                              EM_F[1]+=6*Y_1*w20 + 6*Y_2*w21 + tmp1;
2228                              EM_F[2]+=6*Y_1*w22 + 6*Y_2*w21 + tmp1;                              EM_F[2]+=6*Y_1*w21 + 6*Y_2*w20 + tmp1;
2229                              EM_F[3]+=6*Y_0*w22 + 6*Y_3*w21 + tmp0;                              EM_F[3]+=6*Y_0*w21 + 6*Y_3*w20 + tmp0;
2230                          } else { // constant data                          } else { // constant data
2231                              const double Y_0 = Y_p[0];                              EM_F[0]+=36*Y_p[0]*w22;
2232                              EM_F[0]+=36*Y_0*w23;                              EM_F[1]+=36*Y_p[0]*w22;
2233                              EM_F[1]+=36*Y_0*w23;                              EM_F[2]+=36*Y_p[0]*w22;
2234                              EM_F[2]+=36*Y_0*w23;                              EM_F[3]+=36*Y_p[0]*w22;
                             EM_F[3]+=36*Y_0*w23;  
2235                          }                          }
2236                      }                      }
                     /* GENERATOR SNIP_PDE_SINGLE BOTTOM */  
2237    
2238                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2239                      const index_t firstNode=m_NN[0]*k1+k0;                      const index_t firstNode=m_NN[0]*k1+k0;
# Line 2256  void Rectangle::assemblePDESingleReduced Line 2250  void Rectangle::assemblePDESingleReduced
2250          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2251          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
2252  {  {
2253      const double w0 = .25;      const double w0 = 1./4;
2254      const double w1 = .125*m_dx[0];      const double w1 = m_dx[0]/8;
2255      const double w2 = .125*m_dx[1];      const double w2 = m_dx[1]/8;
2256      const double w3 = .0625*m_dx[0]*m_dx[1];      const double w3 = m_dx[0]*m_dx[1]/16;
2257      const double w4 = .25*m_dx[0]/m_dx[1];      const double w4 = m_dx[0]/(4*m_dx[1]);
2258      const double w5 = .25*m_dx[1]/m_dx[0];      const double w5 = m_dx[1]/(4*m_dx[0]);
2259    
2260      rhs.requireWrite();      rhs.requireWrite();
2261  #pragma omp parallel  #pragma omp parallel
# Line 2386  void Rectangle::assemblePDESingleReduced Line 2380  void Rectangle::assemblePDESingleReduced
2380                      if (!X.isEmpty()) {                      if (!X.isEmpty()) {
2381                          add_EM_F=true;                          add_EM_F=true;
2382                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2383                          const double tmp0 = 4*X_p[0]*w2;                          const double wX0 = 4*X_p[0]*w2;
2384                          const double tmp1 = 4*X_p[1]*w1;                          const double wX1 = 4*X_p[1]*w1;
2385                          EM_F[0]+=-tmp0 - tmp1;                          EM_F[0]+=-wX0 - wX1;
2386                          EM_F[1]+=-tmp1 + tmp0;                          EM_F[1]+=-wX1 + wX0;
2387                          EM_F[2]+=-tmp0 + tmp1;                          EM_F[2]+=-wX0 + wX1;
2388                          EM_F[3]+= tmp0 + tmp1;                          EM_F[3]+= wX0 + wX1;
2389                      }                      }
2390                      ///////////////                      ///////////////
2391                      // process Y //                      // process Y //
# Line 2427  void Rectangle::assemblePDESystem(Paso_S Line 2421  void Rectangle::assemblePDESystem(Paso_S
2421          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
2422          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
2423      }      }
2424        const double SQRT3 = 1.73205080756887719318;
2425      /* GENERATOR SNIP_PDE_SYSTEM_PRE TOP */      const double w1 = 1.0/24;
2426      const double w0 = -0.1555021169820365539*m_dx[1]/m_dx[0];      const double w5 = -SQRT3/24 + 1.0/12;
2427      const double w1 = 0.041666666666666666667;      const double w2 = -SQRT3/24 - 1.0/12;
2428      const double w2 = -0.15550211698203655390;      const double w19 = -m_dx[0]/12;
2429      const double w3 = 0.041666666666666666667*m_dx[0]/m_dx[1];      const double w11 = w19*(SQRT3 + 3)/12;
2430      const double w4 = -0.01116454968463011277*m_dx[1]/m_dx[0];      const double w14 = w19*(-SQRT3 + 3)/12;
2431      const double w5 = 0.011164549684630112770;      const double w16 = w19*(5*SQRT3 + 9)/12;
2432      const double w6 = -0.041666666666666666667*m_dx[1]/m_dx[0];      const double w17 = w19*(-5*SQRT3 + 9)/12;
2433      const double w7 = 0.1555021169820365539*m_dx[0]/m_dx[1];      const double w27 = w19*(-SQRT3 - 3)/2;
2434      const double w8 = 0.01116454968463011277*m_dx[0]/m_dx[1];      const double w28 = w19*(SQRT3 - 3)/2;
2435      const double w9 = -0.25;      const double w18 = -m_dx[1]/12;
2436      const double w10 = -0.16666666666666666667*m_dx[0]/m_dx[1];      const double w10 = w18*(SQRT3 + 3)/12;
2437      const double w11 = -0.032861463941450536761*m_dx[1];      const double w15 = w18*(-SQRT3 + 3)/12;
2438      const double w12 = -0.032861463941450536761*m_dx[0];      const double w12 = w18*(5*SQRT3 + 9)/12;
2439      const double w13 = -0.12264065304058601714*m_dx[1];      const double w13 = w18*(-5*SQRT3 + 9)/12;
2440      const double w14 = -0.0023593469594139828636*m_dx[1];      const double w25 = w18*(-SQRT3 - 3)/2;
2441      const double w15 = -0.008805202725216129906*m_dx[0];      const double w26 = w18*(SQRT3 - 3)/2;
2442      const double w16 = -0.008805202725216129906*m_dx[1];      const double w22 = m_dx[0]*m_dx[1]/144;
2443      const double w17 = -0.12264065304058601714*m_dx[0];      const double w20 = w22*(SQRT3 + 2);
2444      const double w18 = -0.0023593469594139828636*m_dx[0];      const double w21 = w22*(-SQRT3 + 2);
2445      const double w19 = -0.16666666666666666667*m_dx[1];      const double w23 = w22*(4*SQRT3 + 7);
2446      const double w20 = -0.083333333333333333333*m_dx[0];      const double w24 = w22*(-4*SQRT3 + 7);
2447      const double w21 = 0.025917019497006092316*m_dx[0]*m_dx[1];      const double w3 = m_dx[0]/(24*m_dx[1]);
2448      const double w22 = 0.0018607582807716854616*m_dx[0]*m_dx[1];      const double w7 = w3*(SQRT3 + 2);
2449      const double w23 = 0.0069444444444444444444*m_dx[0]*m_dx[1];      const double w8 = w3*(-SQRT3 + 2);
2450      const double w24 = 0.09672363354357992482*m_dx[0]*m_dx[1];      const double w6 = -m_dx[1]/(24*m_dx[0]);
2451      const double w25 = 0.000498588678642297402*m_dx[0]*m_dx[1];      const double w0 = w6*(SQRT3 + 2);
2452      const double w26 = 0.19716878364870322056*m_dx[1];      const double w4 = w6*(-SQRT3 + 2);
     const double w27 = 0.052831216351296779436*m_dx[1];  
     const double w28 = 0.19716878364870322056*m_dx[0];  
     const double w29 = 0.052831216351296779436*m_dx[0];  
     /* GENERATOR SNIP_PDE_SYSTEM_PRE BOTTOM */  
2453    
2454      rhs.requireWrite();      rhs.requireWrite();
2455  #pragma omp parallel  #pragma omp parallel
# Line 2473  void Rectangle::assemblePDESystem(Paso_S Line 2463  void Rectangle::assemblePDESystem(Paso_S
2463                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
2464                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
2465                      const index_t e = k0 + m_NE[0]*k1;                      const index_t e = k0 + m_NE[0]*k1;
                     /* GENERATOR SNIP_PDE_SYSTEM TOP */  
2466                      ///////////////                      ///////////////
2467                      // process A //                      // process A //
2468                      ///////////////                      ///////////////
# Line 2580  void Rectangle::assemblePDESystem(Paso_S Line 2569  void Rectangle::assemblePDESystem(Paso_S
2569                                      const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];                                      const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
2570                                      const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];                                      const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
2571                                      const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];                                      const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
2572                                      const double tmp0 = w9*(-A_01 - A_10);                                      const double tmp0 = 6*w1*(A_01 - A_10);
2573                                      const double tmp1 = w9*(-A_01 + A_10);                                      const double tmp1 = 6*w1*(A_01 + A_10);
2574                                      const double tmp2 = w9*(A_01 + A_10);                                      const double tmp2 = 6*w1*(-A_01 - A_10);
2575                                      const double tmp3 = w9*(A_01 - A_10);                                      const double tmp3 = 6*w1*(-A_01 + A_10);
2576                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp1;
2577                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=8*A_00*w6 + 6*A_01*w1 + A_10*w9 + 4*A_11*w3;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp0;
2578                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp3;
2579                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp2;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp2;
2580                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=8*A_00*w6 - A_11*w10 + tmp3;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp3;
2581                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp2;
2582                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp0;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp1;
2583                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp0;
2584                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp0;
2585                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp0;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp1;
2586                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp2;
2587                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=8*A_00*w6 - A_11*w10 + tmp3;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp3;
2588                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp2;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=4*A_00*w6 - 4*A_11*w3 + tmp2;
2589                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-4*A_00*w6 - 8*A_11*w3 + tmp3;
2590                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=8*A_00*w6 - A_11*w10 + tmp1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=8*A_00*w6 + 4*A_11*w3 + tmp0;
2591                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-8*A_00*w6 + 8*A_11*w3 + tmp1;
2592                                  }                                  }
2593                              }                              }
2594                          }                          }
# Line 2621  void Rectangle::assemblePDESystem(Paso_S Line 2610  void Rectangle::assemblePDESystem(Paso_S
2610                                      const double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];                                      const double B_1_2 = B_p[INDEX4(k,1,m,2, numEq,2,numComp)];
2611                                      const double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];                                      const double B_0_3 = B_p[INDEX4(k,0,m,3, numEq,2,numComp)];
2612                                      const double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];                                      const double B_1_3 = B_p[INDEX4(k,1,m,3, numEq,2,numComp)];
2613                                      const double tmp0 = w15*(B_1_2 + B_1_3);                                      const double tmp0 = w11*(B_1_0 + B_1_1);
2614                                      const double tmp1 = w12*(B_1_0 + B_1_1);                                      const double tmp1 = w14*(B_1_2 + B_1_3);
2615                                      const double tmp2 = w15*(B_1_0 + B_1_1);                                      const double tmp2 = w15*(-B_0_1 - B_0_3);
2616                                      const double tmp3 = w16*(-B_0_1 - B_0_3);                                      const double tmp3 = w10*(-B_0_0 - B_0_2);
2617                                      const double tmp4 = w11*(-B_0_0 - B_0_2);                                      const double tmp4 = w11*(B_1_2 + B_1_3);
2618                                      const double tmp5 = w12*(B_1_2 + B_1_3);                                      const double tmp5 = w14*(B_1_0 + B_1_1);
2619                                      const double tmp6 = w15*(-B_1_0 - B_1_1);                                      const double tmp6 = w11*(-B_1_2 - B_1_3);
2620                                      const double tmp7 = w12*(-B_1_2 - B_1_3);                                      const double tmp7 = w14*(-B_1_0 - B_1_1);
2621                                      const double tmp8 = w15*(-B_1_2 - B_1_3);                                      const double tmp8 = w11*(-B_1_0 - B_1_1);
2622                                      const double tmp9 = w12*(-B_1_0 - B_1_1);                                      const double tmp9 = w14*(-B_1_2 - B_1_3);
2623                                      const double tmp10 = w11*(-B_0_1 - B_0_3);                                      const double tmp10 = w10*(-B_0_1 - B_0_3);
2624                                      const double tmp11 = w16*(-B_0_0 - B_0_2);                                      const double tmp11 = w15*(-B_0_0 - B_0_2);
2625                                      const double tmp12 = w16*(B_0_0 + B_0_2);                                      const double tmp12 = w15*(B_0_0 + B_0_2);
2626                                      const double tmp13 = w11*(B_0_1 + B_0_3);                                      const double tmp13 = w10*(B_0_1 + B_0_3);
2627                                      const double tmp14 = w11*(B_0_0 + B_0_2);                                      const double tmp14 = w10*(B_0_0 + B_0_2);
2628                                      const double tmp15 = w16*(B_0_1 + B_0_3);                                      const double tmp15 = w15*(B_0_1 + B_0_3);
2629                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=B_0_0*w13 + B_0_1*w11 + B_0_2*w16 + B_0_3*w14 + B_1_0*w17 + B_1_1*w15 + B_1_2*w12 + B_1_3*w18;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=B_0_0*w12 + B_0_1*w10 + B_0_2*w15 + B_0_3*w13 + B_1_0*w16 + B_1_1*w14 + B_1_2*w11 + B_1_3*w17;
2630                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=B_0_0*w11 + B_0_1*w13 + B_0_2*w14 + B_0_3*w16 + tmp0 + tmp1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=B_0_0*w10 + B_0_1*w12 + B_0_2*w13 + B_0_3*w15 + tmp0 + tmp1;
2631                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=B_1_0*w12 + B_1_1*w18 + B_1_2*w17 + B_1_3*w15 + tmp14 + tmp15;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=B_1_0*w11 + B_1_1*w17 + B_1_2*w16 + B_1_3*w14 + tmp14 + tmp15;
2632                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp12 + tmp13 + tmp2 + tmp5;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp12 + tmp13 + tmp4 + tmp5;
2633                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-B_0_0*w13 - B_0_1*w11 - B_0_2*w16 - B_0_3*w14 + tmp0 + tmp1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-B_0_0*w12 - B_0_1*w10 - B_0_2*w15 - B_0_3*w13 + tmp0 + tmp1;
2634                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-B_0_0*w11 - B_0_1*w13 - B_0_2*w14 - B_0_3*w16 + B_1_0*w15 + B_1_1*w17 + B_1_2*w18 + B_1_3*w12;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-B_0_0*w10 - B_0_1*w12 - B_0_2*w13 - B_0_3*w15 + B_1_0*w14 + B_1_1*w16 + B_1_2*w17 + B_1_3*w11;
2635                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2 + tmp3 + tmp4 + tmp5;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2636                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=B_1_0*w18 + B_1_1*w12 + B_1_2*w15 + B_1_3*w17 + tmp10 + tmp11;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=B_1_0*w17 + B_1_1*w11 + B_1_2*w14 + B_1_3*w16 + tmp10 + tmp11;
2637                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-B_1_0*w17 - B_1_1*w15 - B_1_2*w12 - B_1_3*w18 + tmp14 + tmp15;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-B_1_0*w16 - B_1_1*w14 - B_1_2*w11 - B_1_3*w17 + tmp14 + tmp15;
2638                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp12 + tmp13 + tmp8 + tmp9;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2639                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=B_0_0*w16 + B_0_1*w14 + B_0_2*w13 + B_0_3*w11 - B_1_0*w12 - B_1_1*w18 - B_1_2*w17 - B_1_3*w15;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=B_0_0*w15 + B_0_1*w13 + B_0_2*w12 + B_0_3*w10 - B_1_0*w11 - B_1_1*w17 - B_1_2*w16 - B_1_3*w14;
2640                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=B_0_0*w14 + B_0_1*w16 + B_0_2*w11 + B_0_3*w13 + tmp6 + tmp7;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=B_0_0*w13 + B_0_1*w15 + B_0_2*w10 + B_0_3*w12 + tmp6 + tmp7;
2641                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp3 + tmp4 + tmp8 + tmp9;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2 + tmp3 + tmp8 + tmp9;
2642                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-B_1_0*w15 - B_1_1*w17 - B_1_2*w18 - B_1_3*w12 + tmp10 + tmp11;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-B_1_0*w14 - B_1_1*w16 - B_1_2*w17 - B_1_3*w11 + tmp10 + tmp11;
2643                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-B_0_0*w16 - B_0_1*w14 - B_0_2*w13 - B_0_3*w11 + tmp6 + tmp7;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-B_0_0*w15 - B_0_1*w13 - B_0_2*w12 - B_0_3*w10 + tmp6 + tmp7;
2644                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-B_0_0*w14 - B_0_1*w16 - B_0_2*w11 - B_0_3*w13 - B_1_0*w18 - B_1_1*w12 - B_1_2*w15 - B_1_3*w17;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-B_0_0*w13 - B_0_1*w15 - B_0_2*w10 - B_0_3*w12 - B_1_0*w17 - B_1_1*w11 - B_1_2*w14 - B_1_3*w16;
2645                                  }                                  }
2646                              }                              }
2647                          } else { // constant data                          } else { // constant data
2648                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2649                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2650                                      const double B_0 = B_p[INDEX3(k,0,m,numEq,2)];                                      const double wB0 = B_p[INDEX3(k,0,m,numEq,2)]*w18;
2651                                      const double B_1 = B_p[INDEX3(k,1,m,numEq,2)];                                      const double wB1 = B_p[INDEX3(k,1,m,numEq,2)]*w19;
2652                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=B_0*w19 + 2*B_1*w20;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+= 2*wB0 + 2*wB1;
2653                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=B_0*w19 + B_1*w20;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+= 2*wB0 +   wB1;
2654                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=B_0*w19/2 + 2*B_1*w20;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=   wB0 + 2*wB1;
2655                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=B_0*w19/2 + B_1*w20;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=   wB0 +   wB1;
2656                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-B_0*w19 + B_1*w20;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-2*wB0 +   wB1;
2657                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-B_0*w19 + 2*B_1*w20;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-2*wB0 + 2*wB1;
2658                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=-B_0*w19/2 + B_1*w20;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=  -wB0 +   wB1;
2659                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-B_0*w19/2 + 2*B_1*w20;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=  -wB0 + 2*wB1;
2660                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=B_0*w19/2 - 2*B_1*w20;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=   wB0 - 2*wB1;
2661                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=B_0*w19/2 - B_1*w20;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=   wB0 -   wB1;
2662                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=B_0*w19 - 2*B_1*w20;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+= 2*wB0 - 2*wB1;
2663                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=B_0*w19 - B_1*w20;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+= 2*wB0 -   wB1;
2664                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=-B_0*w19/2 - B_1*w20;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=  -wB0 -   wB1;
2665                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-B_0*w19/2 - 2*B_1*w20;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=  -wB0 - 2*wB1;
2666                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-B_0*w19 - B_1*w20;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-2*wB0 -   wB1;
2667                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-B_0*w19 - 2*B_1*w20;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-2*wB0 - 2*wB1;
2668                                  }                                  }
2669                              }                              }
2670                          }                          }
# Line 2697  void Rectangle::assemblePDESystem(Paso_S Line 2686  void Rectangle::assemblePDESystem(Paso_S
2686                                      const double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];                                      const double C_1_2 = C_p[INDEX4(k,m,1, 2, numEq,numComp,2)];
2687                                      const double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];                                      const double C_0_3 = C_p[INDEX4(k,m,0, 3, numEq,numComp,2)];
2688                                      const double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];                                      const double C_1_3 = C_p[INDEX4(k,m,1, 3, numEq,numComp,2)];
2689                                      const double tmp0 = w15*(C_1_2 + C_1_3);                                      const double tmp0 = w11*(C_1_0 + C_1_1);
2690                                      const double tmp1 = w12*(C_1_0 + C_1_1);                                      const double tmp1 = w14*(C_1_2 + C_1_3);
2691                                      const double tmp2 = w15*(-C_1_2 - C_1_3);                                      const double tmp2 = w15*(C_0_0 + C_0_2);
2692                                      const double tmp3 = w16*(C_0_0 + C_0_2);                                      const double tmp3 = w10*(C_0_1 + C_0_3);
2693                                      const double tmp4 = w11*(C_0_1 + C_0_3);                                      const double tmp4 = w11*(-C_1_0 - C_1_1);
2694                                      const double tmp5 = w12*(-C_1_0 - C_1_1);                                      const double tmp5 = w14*(-C_1_2 - C_1_3);
2695                                      const double tmp6 = w15*(-C_1_0 - C_1_1);                                      const double tmp6 = w11*(-C_1_2 - C_1_3);
2696                                      const double tmp7 = w12*(-C_1_2 - C_1_3);                                      const double tmp7 = w14*(-C_1_0 - C_1_1);
2697                                      const double tmp8 = w15*(C_1_0 + C_1_1);                                      const double tmp8 = w11*(C_1_2 + C_1_3);
2698                                      const double tmp9 = w12*(C_1_2 + C_1_3);                                      const double tmp9 = w14*(C_1_0 + C_1_1);
2699                                      const double tmp10 = w11*(-C_0_1 - C_0_3);                                      const double tmp10 = w10*(-C_0_1 - C_0_3);
2700                                      const double tmp11 = w16*(-C_0_0 - C_0_2);                                      const double tmp11 = w15*(-C_0_0 - C_0_2);
2701                                      const double tmp12 = w16*(-C_0_1 - C_0_3);                                      const double tmp12 = w15*(-C_0_1 - C_0_3);
2702                                      const double tmp13 = w11*(-C_0_0 - C_0_2);                                      const double tmp13 = w10*(-C_0_0 - C_0_2);
2703                                      const double tmp14 = w11*(C_0_0 + C_0_2);                                      const double tmp14 = w10*(C_0_0 + C_0_2);
2704                                      const double tmp15 = w16*(C_0_1 + C_0_3);                                      const double tmp15 = w15*(C_0_1 + C_0_3);
2705                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=C_0_0*w13 + C_0_1*w11 + C_0_2*w16 + C_0_3*w14 + C_1_0*w17 + C_1_1*w15 + C_1_2*w12 + C_1_3*w18;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=C_0_0*w12 + C_0_1*w10 + C_0_2*w15 + C_0_3*w13 + C_1_0*w16 + C_1_1*w14 + C_1_2*w11 + C_1_3*w17;
2706                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-C_0_0*w13 - C_0_1*w11 - C_0_2*w16 - C_0_3*w14 + tmp0 + tmp1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-C_0_0*w12 - C_0_1*w10 - C_0_2*w15 - C_0_3*w13 + tmp0 + tmp1;
2707                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-C_1_0*w17 - C_1_1*w15 - C_1_2*w12 - C_1_3*w18 + tmp14 + tmp15;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-C_1_0*w16 - C_1_1*w14 - C_1_2*w11 - C_1_3*w17 + tmp14 + tmp15;
2708                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp12 + tmp13 + tmp2 + tmp5;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp12 + tmp13 + tmp4 + tmp5;
2709                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=C_0_0*w11 + C_0_1*w13 + C_0_2*w14 + C_0_3*w16 + tmp0 + tmp1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=C_0_0*w10 + C_0_1*w12 + C_0_2*w13 + C_0_3*w15 + tmp0 + tmp1;
2710                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-C_0_0*w11 - C_0_1*w13 - C_0_2*w14 - C_0_3*w16 + C_1_0*w15 + C_1_1*w17 + C_1_2*w18 + C_1_3*w12;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-C_0_0*w10 - C_0_1*w12 - C_0_2*w13 - C_0_3*w15 + C_1_0*w14 + C_1_1*w16 + C_1_2*w17 + C_1_3*w11;
2711                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2 + tmp3 + tmp4 + tmp5;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2712                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-C_1_0*w15 - C_1_1*w17 - C_1_2*w18 - C_1_3*w12 + tmp10 + tmp11;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-C_1_0*w14 - C_1_1*w16 - C_1_2*w17 - C_1_3*w11 + tmp10 + tmp11;
2713                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=C_1_0*w12 + C_1_1*w18 + C_1_2*w17 + C_1_3*w15 + tmp14 + tmp15;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=C_1_0*w11 + C_1_1*w17 + C_1_2*w16 + C_1_3*w14 + tmp14 + tmp15;
2714                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp12 + tmp13 + tmp8 + tmp9;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2715                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=C_0_0*w16 + C_0_1*w14 + C_0_2*w13 + C_0_3*w11 - C_1_0*w12 - C_1_1*w18 - C_1_2*w17 - C_1_3*w15;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=C_0_0*w15 + C_0_1*w13 + C_0_2*w12 + C_0_3*w10 - C_1_0*w11 - C_1_1*w17 - C_1_2*w16 - C_1_3*w14;
2716                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-C_0_0*w16 - C_0_1*w14 - C_0_2*w13 - C_0_3*w11 + tmp6 + tmp7;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-C_0_0*w15 - C_0_1*w13 - C_0_2*w12 - C_0_3*w10 + tmp6 + tmp7;
2717                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp3 + tmp4 + tmp8 + tmp9;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2 + tmp3 + tmp8 + tmp9;
2718                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=C_1_0*w18 + C_1_1*w12 + C_1_2*w15 + C_1_3*w17 + tmp10 + tmp11;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=C_1_0*w17 + C_1_1*w11 + C_1_2*w14 + C_1_3*w16 + tmp10 + tmp11;
2719                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=C_0_0*w14 + C_0_1*w16 + C_0_2*w11 + C_0_3*w13 + tmp6 + tmp7;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=C_0_0*w13 + C_0_1*w15 + C_0_2*w10 + C_0_3*w12 + tmp6 + tmp7;
2720                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-C_0_0*w14 - C_0_1*w16 - C_0_2*w11 - C_0_3*w13 - C_1_0*w18 - C_1_1*w12 - C_1_2*w15 - C_1_3*w17;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-C_0_0*w13 - C_0_1*w15 - C_0_2*w10 - C_0_3*w12 - C_1_0*w17 - C_1_1*w11 - C_1_2*w14 - C_1_3*w16;
2721                                  }                                  }
2722                              }                              }
2723                          } else { // constant data                          } else { // constant data
2724                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2725                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2726                                      const double C_0 = C_p[INDEX3(k,m,0,numEq,numComp)];                                      const double wC0 = C_p[INDEX3(k,m,0,numEq,numComp)]*w18;
2727                                      const double C_1 = C_p[INDEX3(k,m,1,numEq,numComp)];                                      const double wC1 = C_p[INDEX3(k,m,1,numEq,numComp)]*w19;
2728                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=C_0*w19 + 2*C_1*w20;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+= 2*wC0 + 2*wC1;
2729                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-C_0*w19 + C_1*w20;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-2*wC0 +   wC1;
2730                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=C_0*w19/2 - 2*C_1*w20;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=   wC0 - 2*wC1;
2731                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=-C_0*w19/2 - C_1*w20;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=  -wC0 -   wC1;
2732                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=C_0*w19 + C_1*w20;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+= 2*wC0 +   wC1;
2733                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-C_0*w19 + 2*C_1*w20;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-2*wC0 + 2*wC1;
2734                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=C_0*w19/2 - C_1*w20;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=   wC0 -   wC1;
2735                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-C_0*w19/2 - 2*C_1*w20;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=  -wC0 - 2*wC1;
2736                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=C_0*w19/2 + 2*C_1*w20;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=   wC0 + 2*wC1;
2737                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=-C_0*w19/2 + C_1*w20;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=  -wC0 +   wC1;
2738                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=C_0*w19 - 2*C_1*w20;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+= 2*wC0 - 2*wC1;
2739                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-C_0*w19 - C_1*w20;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-2*wC0 -   wC1;
2740                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=C_0*w19/2 + C_1*w20;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=   wC0 +   wC1;
2741                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-C_0*w19/2 + 2*C_1*w20;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=  -wC0 + 2*wC1;
2742                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=C_0*w19 - C_1*w20;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+= 2*wC0 -   wC1;
2743                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-C_0*w19 - 2*C_1*w20;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-2*wC0 - 2*wC1;
2744                                  }                                  }
2745                              }                              }
2746                          }                          }
# Line 2769  void Rectangle::assemblePDESystem(Paso_S Line 2758  void Rectangle::assemblePDESystem(Paso_S
2758                                      const double D_1 = D_p[INDEX3(k,m,1,numEq,numComp)];                                      const double D_1 = D_p[INDEX3(k,m,1,numEq,numComp)];
2759                                      const double D_2 = D_p[INDEX3(k,m,2,numEq,numComp)];                                      const double D_2 = D_p[INDEX3(k,m,2,numEq,numComp)];
2760                                      const double D_3 = D_p[INDEX3(k,m,3,numEq,numComp)];                                      const double D_3 = D_p[INDEX3(k,m,3,numEq,numComp)];
2761                                      const double tmp0 = w21*(D_0 + D_1);                                      const double tmp0 = w21*(D_2 + D_3);
2762                                      const double tmp1 = w22*(D_2 + D_3);                                      const double tmp1 = w20*(D_0 + D_1);
2763                                      const double tmp2 = w23*(D_0 + D_1 + D_2 + D_3);                                      const double tmp2 = w22*(D_0 + D_1 + D_2 + D_3);
2764                                      const double tmp3 = w21*(D_2 + D_3);                                      const double tmp3 = w21*(D_0 + D_1);
2765                                      const double tmp4 = w22*(D_0 + D_1);                                      const double tmp4 = w20*(D_2 + D_3);
2766                                      const double tmp5 = w23*(D_1 + D_2);                                      const double tmp5 = w22*(D_1 + D_2);
2767                                      const double tmp6 = w21*(D_1 + D_3);                                      const double tmp6 = w21*(D_0 + D_2);
2768                                      const double tmp7 = w22*(D_0 + D_2);                                      const double tmp7 = w20*(D_1 + D_3);
2769                                      const double tmp8 = w21*(D_0 + D_2);                                      const double tmp8 = w21*(D_1 + D_3);
2770                                      const double tmp9 = w22*(D_1 + D_3);                                      const double tmp9 = w20*(D_0 + D_2);
2771                                      const double tmp10 = w23*(D_0 + D_3);                                      const double tmp10 = w22*(D_0 + D_3);
2772                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=D_0*w24 + D_3*w25 + tmp5;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=D_0*w23 + D_3*w24 + tmp5;
2773                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0 + tmp1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0 + tmp1;
2774                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp8 + tmp9;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp8 + tmp9;
2775                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp2;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp2;
2776                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0 + tmp1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0 + tmp1;
2777                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=D_1*w24 + D_2*w25 + tmp10;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=D_1*w23 + D_2*w24 + tmp10;
2778                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2;
2779                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp6 + tmp7;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp6 + tmp7;
2780                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp8 + tmp9;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp8 + tmp9;
2781                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp2;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp2;
2782                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=D_1*w25 + D_2*w24 + tmp10;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=D_1*w24 + D_2*w23 + tmp10;
2783                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp3 + tmp4;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp3 + tmp4;
2784                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2;
2785                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp6 + tmp7;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp6 + tmp7;
2786                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3 + tmp4;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3 + tmp4;
2787                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=D_0*w25 + D_3*w24 + tmp5;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=D_0*w24 + D_3*w23 + tmp5;
2788                                  }                                  }
2789                               }                               }
2790                          } else { // constant data                          } else { // constant data
2791                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2792                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2793                                      const double D_0 = D_p[INDEX2(k, m, numEq)];                                      const double D_0 = D_p[INDEX2(k, m, numEq)];
2794                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=16*D_0*w23;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=16*D_0*w22;
2795                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=8*D_0*w23;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=8*D_0*w22;
2796                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=8*D_0*w23;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=8*D_0*w22;
2797                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=4*D_0*w23;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=4*D_0*w22;
2798                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=8*D_0*w23;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=8*D_0*w22;
2799                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=16*D_0*w23;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=16*D_0*w22;
2800                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=4*D_0*w23;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=4*D_0*w22;
2801                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=8*D_0*w23;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=8*D_0*w22;
2802                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=8*D_0*w23;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=8*D_0*w22;
2803                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=4*D_0*w23;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=4*D_0*w22;
2804                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=16*D_0*w23;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=16*D_0*w22;
2805                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=8*D_0*w23;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=8*D_0*w22;
2806                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=4*D_0*w23;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=4*D_0*w22;
2807                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=8*D_0*w23;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=8*D_0*w22;
2808                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=8*D_0*w23;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=8*D_0*w22;
2809                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=16*D_0*w23;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=16*D_0*w22;
2810                                  }                                  }
2811                              }                              }
2812                          }                          }
# Line 2838  void Rectangle::assemblePDESystem(Paso_S Line 2827  void Rectangle::assemblePDESystem(Paso_S
2827                                  const double X_1_2 = X_p[INDEX3(k,1,2,numEq,2)];                                  const double X_1_2 = X_p[INDEX3(k,1,2,numEq,2)];
2828                                  const double X_0_3 = X_p[INDEX3(k,0,3,numEq,2)];                                  const double X_0_3 = X_p[INDEX3(k,0,3,numEq,2)];
2829                                  const double X_1_3 = X_p[INDEX3(k,1,3,numEq,2)];                                  const double X_1_3 = X_p[INDEX3(k,1,3,numEq,2)];
2830                                  const double tmp0 = 6*w15*(X_1_1 + X_1_3);                                  const double tmp0 = 6*w15*(X_0_2 + X_0_3);
2831                                  const double tmp1 = 6*w16*(X_0_2 + X_0_3);                                  const double tmp1 = 6*w10*(X_0_0 + X_0_1);
2832                                  const double tmp2 = 6*w11*(X_0_0 + X_0_1);                                  const double tmp2 = 6*w11*(X_1_0 + X_1_2);
2833                                  const double tmp3 = 6*w12*(X_1_0 + X_1_2);                                  const double tmp3 = 6*w14*(X_1_1 + X_1_3);
2834                                  const double tmp4 = 6*w15*(X_1_0 + X_1_2);                                  const double tmp4 = 6*w11*(X_1_1 + X_1_3);
2835                                  const double tmp5 = w27*(X_0_2 + X_0_3);                                  const double tmp5 = w25*(X_0_0 + X_0_1);
2836                                  const double tmp6 = w26*(X_0_0 + X_0_1);                                  const double tmp6 = w26*(X_0_2 + X_0_3);
2837                                  const double tmp7 = 6*w12*(X_1_1 + X_1_3);                                  const double tmp7 = 6*w14*(X_1_0 + X_1_2);
2838                                  const double tmp8 = w29*(X_1_1 + X_1_3);                                  const double tmp8 = w27*(X_1_0 + X_1_2);
2839                                  const double tmp9 = w27*(-X_0_0 - X_0_1);                                  const double tmp9 = w28*(X_1_1 + X_1_3);
2840                                  const double tmp10 = w28*(X_1_0 + X_1_2);                                  const double tmp10 = w25*(-X_0_2 - X_0_3);
2841                                  const double tmp11 = w26*(-X_0_2 - X_0_3);                                  const double tmp11 = w26*(-X_0_0 - X_0_1);
2842                                  const double tmp12 = w29*(X_1_0 + X_1_2);                                  const double tmp12 = w27*(X_1_1 + X_1_3);
2843                                  const double tmp13 = w27*(X_0_0 + X_0_1);                                  const double tmp13 = w28*(X_1_0 + X_1_2);
2844                                  const double tmp14 = w28*(X_1_1 + X_1_3);                                  const double tmp14 = w25*(X_0_2 + X_0_3);
2845                                  const double tmp15 = w26*(X_0_2 + X_0_3);                                  const double tmp15 = w26*(X_0_0 + X_0_1);
2846                                  EM_F[INDEX2(k,0,numEq)]+=tmp0 + tmp1 + tmp2 + tmp3;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0 + tmp1 + tmp2 + tmp3;
2847                                  EM_F[INDEX2(k,1,numEq)]+=tmp4 + tmp5 + tmp6 + tmp7;                                  EM_F[INDEX2(k,1,numEq)]+=tmp4 + tmp5 + tmp6 + tmp7;
2848                                  EM_F[INDEX2(k,2,numEq)]+=tmp10 + tmp11 + tmp8 + tmp9;                                  EM_F[INDEX2(k,2,numEq)]+=tmp10 + tmp11 + tmp8 + tmp9;
# Line 2861  void Rectangle::assemblePDESystem(Paso_S Line 2850  void Rectangle::assemblePDESystem(Paso_S
2850                              }                              }
2851                          } else { // constant data                          } else { // constant data
2852                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2853                                  const double X_0 = X_p[INDEX2(k, 0, numEq)];                                  const double wX0 = X_p[INDEX2(k, 0, numEq)]*w18;
2854                                  const double X_1 = X_p[INDEX2(k, 1, numEq)];                                  const double wX1 = X_p[INDEX2(k, 1, numEq)]*w19;
2855                                  EM_F[INDEX2(k,0,numEq)]+=3*X_0*w19 + 6*X_1*w20;                                  EM_F[INDEX2(k,0,numEq)]+= 6*wX0 + 6*wX1;
2856                                  EM_F[INDEX2(k,1,numEq)]+=-3*X_0*w19 + 6*X_1*w20;                                  EM_F[INDEX2(k,1,numEq)]+=-6*wX0 + 6*wX1;
2857                                  EM_F[INDEX2(k,2,numEq)]+=3*X_0*w19 - 6*X_1*w20;                                  EM_F[INDEX2(k,2,numEq)]+= 6*wX0 - 6*wX1;
2858                                  EM_F[INDEX2(k,3,numEq)]+=-3*X_0*w19 - 6*X_1*w20;                                  EM_F[INDEX2(k,3,numEq)]+=-6*wX0 - 6*wX1;
2859                              }                              }
2860                          }                          }
2861                      }                      }
# Line 2882  void Rectangle::assemblePDESystem(Paso_S Line 2871  void Rectangle::assemblePDESystem(Paso_S
2871                                  const double Y_1 = Y_p[INDEX2(k, 1, numEq)];                                  const double Y_1 = Y_p[INDEX2(k, 1, numEq)];
2872                                  const double Y_2 = Y_p[INDEX2(k, 2, numEq)];                                  const double Y_2 = Y_p[INDEX2(k, 2, numEq)];
2873                                  const double Y_3 = Y_p[INDEX2(k, 3, numEq)];                                  const double Y_3 = Y_p[INDEX2(k, 3, numEq)];
2874                                  const double tmp0 = 6*w23*(Y_1 + Y_2);                                  const double tmp0 = 6*w22*(Y_1 + Y_2);
2875                                  const double tmp1 = 6*w23*(Y_0 + Y_3);                                  const double tmp1 = 6*w22*(Y_0 + Y_3);
2876                                  EM_F[INDEX2(k,0,numEq)]+=6*Y_0*w21 + 6*Y_3*w22 + tmp0;                                  EM_F[INDEX2(k,0,numEq)]+=6*Y_0*w20 + 6*Y_3*w21 + tmp0;
2877                                  EM_F[INDEX2(k,1,numEq)]+=6*Y_1*w21 + 6*Y_2*w22 + tmp1;                                  EM_F[INDEX2(k,1,numEq)]+=6*Y_1*w20 + 6*Y_2*w21 + tmp1;
2878                                  EM_F[INDEX2(k,2,numEq)]+=6*Y_1*w22 + 6*Y_2*w21 + tmp1;                                  EM_F[INDEX2(k,2,numEq)]+=6*Y_1*w21 + 6*Y_2*w20 + tmp1;
2879                                  EM_F[INDEX2(k,3,numEq)]+=6*Y_0*w22 + 6*Y_3*w21 + tmp0;                                  EM_F[INDEX2(k,3,numEq)]+=6*Y_0*w21 + 6*Y_3*w20 + tmp0;
2880                              }                              }
2881                          } else { // constant data                          } else { // constant data
2882                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2883                                  const double Y_0 = Y_p[k];                                  EM_F[INDEX2(k,0,numEq)]+=36*Y_p[k]*w22;
2884                                  EM_F[INDEX2(k,0,numEq)]+=36*Y_0*w23;                                  EM_F[INDEX2(k,1,numEq)]+=36*Y_p[k]*w22;
2885                                  EM_F[INDEX2(k,1,numEq)]+=36*Y_0*w23;                                  EM_F[INDEX2(k,2,numEq)]+=36*Y_p[k]*w22;
2886                                  EM_F[INDEX2(k,2,numEq)]+=36*Y_0*w23;                                  EM_F[INDEX2(k,3,numEq)]+=36*Y_p[k]*w22;
                                 EM_F[INDEX2(k,3,numEq)]+=36*Y_0*w23;  
2887                              }                              }
2888                          }                          }
2889                      }                      }
                     /* GENERATOR SNIP_PDE_SYSTEM BOTTOM */  
2890    
2891                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2892                      const index_t firstNode=m_NN[0]*k1+k0;                      const index_t firstNode=m_NN[0]*k1+k0;
# Line 2925  void Rectangle::assemblePDESystemReduced Line 2912  void Rectangle::assemblePDESystemReduced
2912          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
2913      }      }
2914    
2915      const double w0 = .25;      const double w0 = 1./4;
2916      const double w1 = .125*m_dx[0];      const double w1 = m_dx[0]/8;
2917      const double w2 = .125*m_dx[1];      const double w2 = m_dx[1]/8;
2918      const double w3 = .0625*m_dx[0]*m_dx[1];      const double w3 = m_dx[0]*m_dx[1]/16;
2919      const double w4 = .25*m_dx[0]/m_dx[1];      const double w4 = m_dx[0]/(4*m_dx[1]);
2920      const double w5 = .25*m_dx[1]/m_dx[0];      const double w5 = m_dx[1]/(4*m_dx[0]);
2921    
2922      rhs.requireWrite();      rhs.requireWrite();
2923  #pragma omp parallel  #pragma omp parallel
# Line 2952  void Rectangle::assemblePDESystemReduced Line 2939  void Rectangle::assemblePDESystemReduced
2939                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2940                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
2941                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
2942                                  const double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];                                  const double Aw00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)]*w5;
2943                                  const double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];                                  const double Aw01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)]*w0;
2944                                  const double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];                                  const double Aw10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)]*w0;
2945                                  const double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];                                  const double Aw11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)]*w4;
2946                                  const double tmp0 = (A_01 + A_10)*w0;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+= Aw00 + Aw01 + Aw10 + Aw11;
2947                                  const double tmp1 = A_00*w5;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-Aw00 - Aw01 + Aw10 + Aw11;
2948                                  const double tmp2 = A_01*w0;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+= Aw00 + Aw01 - Aw10 - Aw11;
2949                                  const double tmp3 = A_10*w0;                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=-Aw00 - Aw01 - Aw10 - Aw11;
2950                                  const double tmp4 = A_11*w4;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-Aw00 + Aw01 - Aw10 + Aw11;
2951                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+= tmp4 + tmp0 + tmp1;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+= Aw00 - Aw01 - Aw10 + Aw11;
2952                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+= tmp4 - tmp1 + tmp3 - tmp2;                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=-Aw00 + Aw01 + Aw10 - Aw11;
2953                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+= tmp2 - tmp3 - tmp4 + tmp1;                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+= Aw00 - Aw01 + Aw10 - Aw11;
2954                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=-tmp1 - tmp4 - tmp0;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+= Aw00 - Aw01 + Aw10 - Aw11;
2955                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+= tmp4 - tmp1 + tmp2 - tmp3;                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=-Aw00 + Aw01 + Aw10 - Aw11;
2956                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+= tmp4 + tmp1 - tmp0;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+= Aw00 - Aw01 - Aw10 + Aw11;
2957                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=-tmp1 + tmp0 - tmp4;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-Aw00 + Aw01 - Aw10 + Aw11;
2958                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=-Aw00 - Aw01 - Aw10 - Aw11;
2959                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+= Aw00 + Aw01 - Aw10 - Aw11;
2960                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=-tmp1 + tmp0 - tmp4;                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-Aw00 - Aw01 + Aw10 + Aw11;
2961                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+= tmp4 + tmp1 - tmp0;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+= Aw00 + Aw01 + Aw10 + Aw11;
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+= tmp4 - tmp1 + tmp2 - tmp3;  
                                 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=-tmp1 - tmp4 - tmp0;  
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+= tmp2 - tmp3 - tmp4 + tmp1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+= tmp4 - tmp1 + tmp3 - tmp2;  
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+= tmp4 + tmp0 + tmp1;  
2962                              }                              }
2963                          }                          }
2964                      }                      }
# Line 2988  void Rectangle::assemblePDESystemReduced Line 2970  void Rectangle::assemblePDESystemReduced
2970                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2971                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
2972                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
2973                                  const double tmp0 = B_p[INDEX3(k,0,m, numEq, 2)]*w2;                                  const double wB0 = B_p[INDEX3(k,0,m, numEq, 2)]*w2;
2974                                  const double tmp1 = B_p[INDEX3(k,1,m, numEq, 2)]*w1;                                  const double wB1 = B_p[INDEX3(k,1,m, numEq, 2)]*w1;
2975                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-tmp0 - tmp1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-wB0 - wB1;
2976                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-tmp1 + tmp0;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-wB0 - wB1;
2977                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-tmp0 + tmp1;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-wB0 - wB1;
2978                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+= tmp0 + tmp1;                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=-wB0 - wB1;
2979                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-tmp0 - tmp1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+= wB0 - wB1;
2980                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-tmp1 + tmp0;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+= wB0 - wB1;
2981                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=-tmp0 + tmp1;                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+= wB0 - wB1;
2982                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+= tmp0 + tmp1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+= wB0 - wB1;
2983                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-tmp0 - tmp1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-wB0 + wB1;
2984                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=-tmp1 + tmp0;                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=-wB0 + wB1;
2985                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=-tmp0 + tmp1;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=-wB0 + wB1;
2986                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+= tmp0 + tmp1;                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-wB0 + wB1;
2987                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=-tmp0 - tmp1;                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+= wB0 + wB1;
2988                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-tmp1 + tmp0;                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+= wB0 + wB1;
2989                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-tmp0 + tmp1;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+= wB0 + wB1;
2990                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+= tmp0 + tmp1;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+= wB0 + wB1;
2991                              }                              }
2992                          }                          }
2993                      }                      }
# Line 3017  void Rectangle::assemblePDESystemReduced Line 2999  void Rectangle::assemblePDESystemReduced
2999                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
3000                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3001                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3002                                  const double tmp0 = C_p[INDEX3(k, m, 0, numEq, numComp)]*w2;                                  const double wC0 = C_p[INDEX3(k, m, 0, numEq, numComp)]*w2;
3003                                  const double tmp1 = C_p[INDEX3(k, m, 1, numEq, numComp)]*w1;                                  const double wC1 = C_p[INDEX3(k, m, 1, numEq, numComp)]*w1;
3004                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-tmp1 - tmp0;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-wC0 - wC1;
3005                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-tmp1 - tmp0;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-wC0 - wC1;
3006                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-tmp1 - tmp0;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-wC0 - wC1;
3007                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=-tmp1 - tmp0;                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=-wC0 - wC1;
3008                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+= tmp0 - tmp1;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+= wC0 - wC1;
3009                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+= tmp0 - tmp1;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+= wC0 - wC1;
3010                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+= tmp0 - tmp1;                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+= wC0 - wC1;
3011                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+= tmp0 - tmp1;                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+= wC0 - wC1;
3012                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+= tmp1 - tmp0;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-wC0 + wC1;
3013                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+= tmp1 - tmp0;                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=-wC0 + wC1;
3014                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+= tmp1 - tmp0;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=-wC0 + wC1;
3015                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+= tmp1 - tmp0;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-wC0 + wC1;
3016                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+= tmp0 + tmp1;                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+= wC0 + wC1;
3017                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+= tmp0 + tmp1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+= wC0 + wC1;
3018                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+= tmp0 + tmp1;                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+= wC0 + wC1;
3019                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+= tmp0 + tmp1;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+= wC0 + wC1;
3020                              }                              }
3021                          }                          }
3022                      }                      }
# Line 3046  void Rectangle::assemblePDESystemReduced Line 3028  void Rectangle::assemblePDESystemReduced
3028                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
3029                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3030                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3031                                  const double tmp0 = D_p[INDEX2(k, m, numEq)]*w3;                                  const double wD0 = D_p[INDEX2(k, m, numEq)]*w3;
3032                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=wD0;
3033                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=wD0;
3034                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=wD0;
3035                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=wD0;
3036                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=wD0;
3037                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=wD0;
3038                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=wD0;
3039                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=wD0;
3040                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=wD0;
3041                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=wD0;
3042                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=wD0;
3043                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=wD0;
3044                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=wD0;
3045                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=wD0;
3046                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=wD0;
3047                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=wD0;
3048                              }                              }
3049                          }                          }
3050                      }                      }
# Line 3073  void Rectangle::assemblePDESystemReduced Line 3055  void Rectangle::assemblePDESystemReduced
3055                          add_EM_F=true;                          add_EM_F=true;
3056                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
3057                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3058                              const double X_0 = X_p[INDEX2(k, 0, numEq)];                              const double wX0 = 4*X_p[INDEX2(k, 0, numEq)]*w2;
3059                              const double X_1 = X_p[INDEX2(k, 1, numEq)];                              const double wX1 = 4*X_p[INDEX2(k, 1, numEq)]*w1;
3060                              const double tmp0 = 4*X_0*w2;                              EM_F[INDEX2(k,0,numEq)]+=-wX0 - wX1;
3061                              const double tmp1 = 4*X_1*w1;                              EM_F[INDEX2(k,1,numEq)]+= wX0 - wX1;
3062                              EM_F[INDEX2(k,0,numEq)]+=-tmp0 - tmp1;                              EM_F[INDEX2(k,2,numEq)]+=-wX0 + wX1;
3063                              EM_F[INDEX2(k,1,numEq)]+= tmp0 - tmp1;                              EM_F[INDEX2(k,3,numEq)]+= wX0 + wX1;
                             EM_F[INDEX2(k,2,numEq)]+= tmp1 - tmp0;  
                             EM_F[INDEX2(k,3,numEq)]+= tmp0 + tmp1;  
3064                          }                          }
3065                      }                      }
3066                      ///////////////                      ///////////////
# Line 3111  void Rectangle::assemblePDESystemReduced Line 3091  void Rectangle::assemblePDESystemReduced
3091  void Rectangle::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
3092        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3093  {  {
3094      const double w0 = 0.31100423396407310779*m_dx[1];      const double SQRT3 = 1.73205080756887719318;
3095      const double w1 = 0.022329099369260225539*m_dx[1];      const double w5 = m_dx[0]/12;
3096      const double w10 = 0.022329099369260225539*m_dx[0];      const double w6 = w5*(SQRT3 + 2);
3097      const double w11 = 0.16666666666666666667*m_dx[0];      const double w7 = w5*(-SQRT3 + 2);
3098      const double w12 = 0.33333333333333333333*m_dx[0];      const double w8 = w5*(SQRT3 + 3);
3099      const double w13 = 0.39433756729740644113*m_dx[0];      const double w9 = w5*(-SQRT3 + 3);
3100      const double w14 = 0.10566243270259355887*m_dx[0];      const double w2 = m_dx[1]/12;
3101      const double w15 = 0.5*m_dx[0];      const double w0 = w2*(SQRT3 + 2);
3102      const double w2 = 0.083333333333333333333*m_dx[1];      const double w1 = w2*(-SQRT3 + 2);
3103      const double w3 = 0.33333333333333333333*m_dx[1];      const double w3 = w2*(SQRT3 + 3);
3104      const double w4 = 0.16666666666666666667*m_dx[1];      const double w4 = w2*(-SQRT3 + 3);
     const double w5 = 0.39433756729740644113*m_dx[1];  
     const double w6 = 0.10566243270259355887*m_dx[1];  
     const double w7 = 0.5*m_dx[1];  
     const double w8 = 0.083333333333333333333*m_dx[0];  
     const double w9 = 0.31100423396407310779*m_dx[0];  
3105      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3106      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3107      rhs.requireWrite();      rhs.requireWrite();
# Line 3138  void Rectangle::assemblePDEBoundarySingl Line 3113  void Rectangle::assemblePDEBoundarySingl
3113                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
3114                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3115                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
                     const index_t e = k1;  
3116                      ///////////////                      ///////////////
3117                      // process d //                      // process d //
3118                      ///////////////                      ///////////////
3119                      if (add_EM_S) {                      if (add_EM_S) {
3120                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(k1);
3121                          if (d.actsExpanded()) {                          if (d.actsExpanded()) {
3122                              const double d_0 = d_p[0];                              const double d_0 = d_p[0];
3123                              const double d_1 = d_p[1];                              const double d_1 = d_p[1];
3124                              const double tmp0_0 = d_0 + d_1;                              const double tmp0 = w2*(d_0 + d_1);
3125                              const double tmp1_1 = d_1*w1;                              EM_S[INDEX2(0,0,4)]+=d_0*w0 + d_1*w1;
3126                              const double tmp4_1 = d_0*w1;                              EM_S[INDEX2(2,0,4)]+=tmp0;
3127                              const double tmp0_1 = d_0*w0;                              EM_S[INDEX2(0,2,4)]+=tmp0;
3128                              const double tmp3_1 = d_1*w0;                              EM_S[INDEX2(2,2,4)]+=d_0*w1 + d_1*w0;
3129                              const double tmp2_1 = tmp0_0*w2;                          } else { // constant data
3130                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(0,0,4)]+=4*d_p[0]*w2;
3131                              EM_S[INDEX2(2,0,4)]+=tmp2_1;                              EM_S[INDEX2(2,0,4)]+=2*d_p[0]*w2;
3132                              EM_S[INDEX2(0,2,4)]+=tmp2_1;                              EM_S[INDEX2(0,2,4)]+=2*d_p[0]*w2;
3133                              EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(2,2,4)]+=4*d_p[0]*w2;
                         } else { /* constant data */  
                             const double d_0 = d_p[0];  
                             const double tmp1_1 = d_0*w4;  
                             const double tmp0_1 = d_0*w3;  
                             EM_S[INDEX2(0,0,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp1_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp1_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp0_1;  
3134                          }                          }
3135                      }                      }
3136                      ///////////////                      ///////////////
3137                      // process y //                      // process y //
3138                      ///////////////                      ///////////////
3139                      if (add_EM_F) {                      if (add_EM_F) {
3140                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(k1);
3141                          if (y.actsExpanded()) {                          if (y.actsExpanded()) {
3142                              const double y_0 = y_p[0];                              EM_F[0]+=w3*y_p[0] + w4*y_p[1];
3143                              const double y_1 = y_p[1];                              EM_F[2]+=w3*y_p[1] + w4*y_p[0];
3144                              const double tmp3_1 = w5*y_1;                          } else { // constant data
3145                              const double tmp2_1 = w6*y_0;                              EM_F[0]+=6*w2*y_p[0];
3146                              const double tmp0_1 = w6*y_1;                              EM_F[2]+=6*w2*y_p[0];
                             const double tmp1_1 = w5*y_0;  
                             EM_F[0]+=tmp0_1 + tmp1_1;  
                             EM_F[2]+=tmp2_1 + tmp3_1;  
                         } else { /* constant data */  
                             const double y_0 = y_p[0];  
                             const double tmp0_1 = w7*y_0;  
                             EM_F[0]+=tmp0_1;  
                             EM_F[2]+=tmp0_1;  
3147                          }                          }
3148                      }                      }
3149                      const index_t firstNode=m_NN[0]*k1;                      const index_t firstNode=m_NN[0]*k1;
# Line 3209  void Rectangle::assemblePDEBoundarySingl Line 3167  void Rectangle::assemblePDEBoundarySingl
3167                          if (d.actsExpanded()) {                          if (d.actsExpanded()) {
3168                              const double d_0 = d_p[0];                              const double d_0 = d_p[0];
3169                              const double d_1 = d_p[1];                              const double d_1 = d_p[1];
3170                              const double tmp0_0 = d_0 + d_1;                              const double tmp0 = w2*(d_0 + d_1);
3171                              const double tmp4_1 = d_1*w1;                              EM_S[INDEX2(1,1,4)]+=d_0*w0 + d_1*w1;
3172                              const double tmp1_1 = d_0*w1;                              EM_S[INDEX2(3,1,4)]+=tmp0;
3173                              const double tmp3_1 = d_0*w0;                              EM_S[INDEX2(1,3,4)]+=tmp0;
3174                              const double tmp0_1 = d_1*w0;                              EM_S[INDEX2(3,3,4)]+=d_0*w1 + d_1*w0;
3175                              const double tmp2_1 = tmp0_0*w2;                          } else { // constant data
3176                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(1,1,4)]+=4*d_p[0]*w2;
3177                              EM_S[INDEX2(3,1,4)]+=tmp2_1;                              EM_S[INDEX2(3,1,4)]+=2*d_p[0]*w2;
3178                              EM_S[INDEX2(1,3,4)]+=tmp2_1;                              EM_S[INDEX2(1,3,4)]+=2*d_p[0]*w2;
3179                              EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(3,3,4)]+=4*d_p[0]*w2;
                         } else { /* constant data */  
                             const double d_0 = d_p[0];  
                             const double tmp1_1 = d_0*w4;  
                             const double tmp0_1 = d_0*w3;  
                             EM_S[INDEX2(1,1,4)]+=tmp0_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp1_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp1_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp0_1;  
3180                          }                          }
3181                      }                      }
3182                      ///////////////                      ///////////////
# Line 3235  void Rectangle::assemblePDEBoundarySingl Line 3185  void Rectangle::assemblePDEBoundarySingl
3185                      if (add_EM_F) {                      if (add_EM_F) {
3186                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3187                          if (y.actsExpanded()) {                          if (y.actsExpanded()) {
3188                              const double y_0 = y_p[0];                              EM_F[1]+=w3*y_p[0] + w4*y_p[1];
3189                              const double y_1 = y_p[1];                              EM_F[3]+=w3*y_p[1] + w4*y_p[0];
3190                              const double tmp3_1 = w5*y_1;                          } else { // constant data
3191                              const double tmp2_1 = w6*y_0;                              EM_F[1]+=6*w2*y_p[0];
3192                              const double tmp0_1 = w6*y_1;                              EM_F[3]+=6*w2*y_p[0];
                             const double tmp1_1 = w5*y_0;  
                             EM_F[1]+=tmp0_1 + tmp1_1;  
                             EM_F[3]+=tmp2_1 + tmp3_1;  
                         } else { /* constant data */  
                             const double y_0 = y_p[0];  
                             const double tmp0_1 = w7*y_0;  
                             EM_F[1]+=tmp0_1;  
                             EM_F[3]+=tmp0_1;  
3193                          }                          }
3194                      }                      }
3195                      const index_t firstNode=m_NN[0]*(k1+1)-2;                      const index_t firstNode=m_NN[0]*(k1+1)-2;
# Line 3271  void Rectangle::assemblePDEBoundarySingl Line 3213  void Rectangle::assemblePDEBoundarySingl
3213                          if (d.actsExpanded()) {                          if (d.actsExpanded()) {
3214                              const double d_0 = d_p[0];                              const double d_0 = d_p[0];
3215                              const double d_1 = d_p[1];                              const double d_1 = d_p[1];
3216                              const double tmp0_0 = d_0 + d_1;                              const double tmp0 = w5*(d_0 + d_1);
3217                              const double tmp4_1 = d_1*w9;                              EM_S[INDEX2(0,0,4)]+=d_0*w6 + d_1*w7;
3218                              const double tmp2_1 = d_0*w9;                              EM_S[INDEX2(1,0,4)]+=tmp0;
3219                              const double tmp0_1 = tmp0_0*w8;                              EM_S[INDEX2(0,1,4)]+=tmp0;
3220                              const double tmp1_1 = d_1*w10;                              EM_S[INDEX2(1,1,4)]+=d_0*w7 + d_1*w6;
3221                              const double tmp3_1 = d_0*w10;                          } else { // constant data
3222                              EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp2_1;                              EM_S[INDEX2(0,0,4)]+=4*d_p[0]*w5;
3223                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(1,0,4)]+=2*d_p[0]*w5;
3224                              EM_S[INDEX2(0,1,4)]+=tmp0_1;                              EM_S[INDEX2(0,1,4)]+=2*d_p[0]*w5;
3225                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(1,1,4)]+=4*d_p[0]*w5;
                         } else { /* constant data */  
                             const double d_0 = d_p[0];  
                             const double tmp0_1 = d_0*w11;  
                             const double tmp1_1 = d_0*w12;  
                             EM_S[INDEX2(0,0,4)]+=tmp1_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp1_1;  
3226                          }                          }
3227                      }                      }
3228                      ///////////////                      ///////////////
# Line 3297  void Rectangle::assemblePDEBoundarySingl Line 3231  void Rectangle::assemblePDEBoundarySingl
3231                      if (add_EM_F) {                      if (add_EM_F) {
3232                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3233                          if (y.actsExpanded()) {                          if (y.actsExpanded()) {
3234                              const double y_0 = y_p[0];                              EM_F[0]+=w8*y_p[0] + w9*y_p[1];
3235                              const double y_1 = y_p[1];                              EM_F[1]+=w8*y_p[1] + w9*y_p[0];
3236                              const double tmp2_1 = w13*y_1;                          } else { // constant data
3237                              const double tmp1_1 = w14*y_1;                              EM_F[0]+=6*w5*y_p[0];
3238                              const double tmp3_1 = w14*y_0;                              EM_F[1]+=6*w5*y_p[0];
                             const double tmp0_1 = w13*y_0;  
                             EM_F[0]+=tmp0_1 + tmp1_1;  
                             EM_F[1]+=tmp2_1 + tmp3_1;  
                         } else { /* constant data */  
                             const double y_0 = y_p[0];  
                             const double tmp0_1 = w15*y_0;  
                             EM_F[0]+=tmp0_1;  
                             EM_F[1]+=tmp0_1;  
3239                          }                          }
3240                      }                      }
3241                      const index_t firstNode=k0;                      const index_t firstNode=k0;
# Line 3333  void Rectangle::assemblePDEBoundarySingl Line 3259  void Rectangle::assemblePDEBoundarySingl
3259                          if (d.actsExpanded()) {                          if (d.actsExpanded()) {
3260                              const double d_0 = d_p[0];                              const double d_0 = d_p[0];
3261                              const double d_1 = d_p[1];                              const double d_1 = d_p[1];
3262                              const double tmp0_0 = d_0 + d_1;                              const double tmp0 = w5*(d_0 + d_1);
3263                              const double tmp2_1 = d_1*w9;                              EM_S[INDEX2(2,2,4)]+=d_0*w6 + d_1*w7;
3264                              const double tmp4_1 = d_0*w9;                              EM_S[INDEX2(3,2,4)]+=tmp0;
3265                              const double tmp0_1 = tmp0_0*w8;                              EM_S[INDEX2(2,3,4)]+=tmp0;
3266                              const double tmp3_1 = d_1*w10;                              EM_S[INDEX2(3,3,4)]+=d_0*w7 + d_1*w6;
3267                              const double tmp1_1 = d_0*w10;                          } else { // constant data
3268                              EM_S[INDEX2(2,2,4)]+=tmp3_1 + tmp4_1;                              EM_S[INDEX2(2,2,4)]+=4*d_p[0]*w5;
3269                              EM_S[INDEX2(3,2,4)]+=tmp0_1;                              EM_S[INDEX2(3,2,4)]+=2*d_p[0]*w5;
3270                              EM_S[INDEX2(2,3,4)]+=tmp0_1;                              EM_S[INDEX2(2,3,4)]+=2*d_p[0]*w5;
3271                              EM_S[INDEX2(3,3,4)]+=tmp1_1 + tmp2_1;                              EM_S[INDEX2(3,3,4)]+=4*d_p[0]*w5;
                         } else { /* constant data */  
                             const double d_0 = d_p[0];  
                             const double tmp0_1 = d_0*w11;  
                             const double tmp1_1 = d_0*w12;  
                             EM_S[INDEX2(2,2,4)]+=tmp1_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp1_1;  
3272                          }                          }
3273                      }                      }
3274                      ///////////////                      ///////////////
# Line 3359  void Rectangle::assemblePDEBoundarySingl Line 3277  void Rectangle::assemblePDEBoundarySingl
3277                      if (add_EM_F) {                      if (add_EM_F) {
3278                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3279                          if (y.actsExpanded()) {                          if (y.actsExpanded()) {
3280                              const double y_0 = y_p[0];                              EM_F[2]+=w8*y_p[0] + w9*y_p[1];
3281                              const double y_1 = y_p[1];                              EM_F[3]+=w8*y_p[1] + w9*y_p[0];
3282                              const double tmp2_1 = w13*y_1;                          } else { // constant data
3283                              const double tmp1_1 = w14*y_1;                              EM_F[2]+=6*w5*y_p[0];
3284                              const double tmp3_1 = w14*y_0;                              EM_F[3]+=6*w5*y_p[0];
                             const double tmp0_1 = w13*y_0;  
                             EM_F[2]+=tmp0_1 + tmp1_1;  
                             EM_F[3]+=tmp2_1 + tmp3_1;  
                         } else { /* constant data */  
                             const double y_0 = y_p[0];  
                             const double tmp0_1 = w15*y_0;  
                             EM_F[2]+=tmp0_1;  
                             EM_F[3]+=tmp0_1;  
3285                          }                          }
3286                      }                      }
3287                      const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;                      const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;
# Line 3386  void Rectangle::assemblePDEBoundarySingl Line 3296  void Rectangle::assemblePDEBoundarySingl
3296  void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
3297        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3298  {  {
3299      const double w0 = 0.25*m_dx[0];      const double w0 = m_dx[0]/4;
3300      const double w1 = 0.25*m_dx[1];      const double w1 = m_dx[1]/4;
3301      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3302      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3303      rhs.requireWrite();      rhs.requireWrite();
# Line 3404  void Rectangle::assemblePDEBoundarySingl Line 3314  void Rectangle::assemblePDEBoundarySingl
3314                      ///////////////                      ///////////////
3315                      if (add_EM_S) {                      if (add_EM_S) {
3316                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(k1);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(k1);
3317                          const double tmp0 = d_p[0]*w1;                          EM_S[INDEX2(0,0,4)]+=d_p[0]*w1;
3318                          EM_S[INDEX2(0,0,4)]+=tmp0;                          EM_S[INDEX2(2,0,4)]+=d_p[0]*w1;
3319                          EM_S[INDEX2(2,0,4)]+=tmp0;                          EM_S[INDEX2(0,2,4)]+=d_p[0]*w1;
3320                          EM_S[INDEX2(0,2,4)]+=tmp0;                          EM_S[INDEX2(2,2,4)]+=d_p[0]*w1;
                         EM_S[INDEX2(2,2,4)]+=tmp0;  
3321                      }                      }
3322                      ///////////////                      ///////////////
3323                      // process y //                      // process y //
3324                      ///////////////                      ///////////////
3325                      if (add_EM_F) {                      if (add_EM_F) {
3326                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(k1);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(k1);
3327                          const double tmp0 = 2*w1*y_p[0];                          EM_F[0]+=2*w1*y_p[0];
3328                          EM_F[0]+=tmp0;                          EM_F[2]+=2*w1*y_p[0];
                         EM_F[2]+=tmp0;  
3329                      }                      }
3330                      const index_t firstNode=m_NN[0]*k1;                      const index_t firstNode=m_NN[0]*k1;
3331                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
# Line 3437  void Rectangle::assemblePDEBoundarySingl Line 3345  void Rectangle::assemblePDEBoundarySingl
3345                      ///////////////                      ///////////////
3346                      if (add_EM_S) {                      if (add_EM_S) {
3347                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3348                          const double tmp0 = d_p[0]*w1;                          EM_S[INDEX2(1,1,4)]+=d_p[0]*w1;
3349                          EM_S[INDEX2(1,1,4)]+=tmp0;                          EM_S[INDEX2(3,1,4)]+=d_p[0]*w1;
3350                          EM_S[INDEX2(3,1,4)]+=tmp0;                          EM_S[INDEX2(1,3,4)]+=d_p[0]*w1;
3351                          EM_S[INDEX2(1,3,4)]+=tmp0;                          EM_S[INDEX2(3,3,4)]+=d_p[0]*w1;
                         EM_S[INDEX2(3,3,4)]+=tmp0;  
3352                      }                      }
3353                      ///////////////                      ///////////////
3354                      // process y //                      // process y //
3355                      ///////////////                      ///////////////
3356                      if (add_EM_F) {                      if (add_EM_F) {
3357                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3358                          const double tmp0 = 2*w1*y_p[0];                          EM_F[1]+=2*w1*y_p[0];
3359                          EM_F[1]+=tmp0;                          EM_F[3]+=2*w1*y_p[0];
                         EM_F[3]+=tmp0;  
3360                      }                      }
3361                      const index_t firstNode=m_NN[0]*(k1+1)-2;                      const index_t firstNode=m_NN[0]*(k1+1)-2;
3362                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
# Line 3470  void Rectangle::assemblePDEBoundarySingl Line 3376  void Rectangle::assemblePDEBoundarySingl
3376                      ///////////////                      ///////////////
3377                      if (add_EM_S) {                      if (add_EM_S) {
3378                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3379                          const double tmp0 = d_p[0]*w0;                          EM_S[INDEX2(0,0,4)]+=d_p[0]*w0;
3380                          EM_S[INDEX2(0,0,4)]+=tmp0;                          EM_S[INDEX2(1,0,4)]+=d_p[0]*w0;
3381                          EM_S[INDEX2(1,0,4)]+=tmp0;                          EM_S[INDEX2(0,1,4)]+=d_p[0]*w0;
3382                          EM_S[INDEX2(0,1,4)]+=tmp0;                          EM_S[INDEX2(1,1,4)]+=d_p[0]*w0;
                         EM_S[INDEX2(1,1,4)]+=tmp0;  
3383                      }                      }
3384                      ///////////////                      ///////////////
3385                      // process y //                      // process y //
3386                      ///////////////                      ///////////////
3387                      if (add_EM_F) {                      if (add_EM_F) {
3388                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3389                          const double tmp0 = 2*w0*y_p[0];                          EM_F[0]+=2*w0*y_p[0];
3390                          EM_F[0]+=tmp0;                          EM_F[1]+=2*w0*y_p[0];
                         EM_F[1]+=tmp0;  
3391                      }                      }
3392                      const index_t firstNode=k0;                      const index_t firstNode=k0;
3393                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
# Line 3503  void Rectangle::assemblePDEBoundarySingl Line 3407  void Rectangle::assemblePDEBoundarySingl
3407                      ///////////////                      ///////////////
3408                      if (add_EM_S) {                      if (add_EM_S) {
3409                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3410                          const double tmp0 = d_p[0]*w0;                          EM_S[INDEX2(2,2,4)]+=d_p[0]*w0;
3411                          EM_S[INDEX2(2,2,4)]+=tmp0;                          EM_S[INDEX2(3,2,4)]+=d_p[0]*w0;
3412                          EM_S[INDEX2(3,2,4)]+=tmp0;                          EM_S[INDEX2(2,3,4)]+=d_p[0]*w0;
3413                          EM_S[INDEX2(2,3,4)]+=tmp0;                          EM_S[INDEX2(3,3,4)]+=d_p[0]*w0;
                         EM_S[INDEX2(3,3,4)]+=tmp0;  
3414                      }                      }
3415                      ///////////////                      ///////////////
3416                      // process y //                      // process y //
3417                      ///////////////                      ///////////////
3418                      if (add_EM_F) {                      if (add_EM_F) {
3419                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3420                          const double tmp0 = 2*w0*y_p[0];                          EM_F[2]+=2*w0*y_p[0];
3421                          EM_F[2]+=tmp0;                          EM_F[3]+=2*w0*y_p[0];
                         EM_F[3]+=tmp0;  
3422                      }                      }
3423                      const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;                      const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;
3424                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);                      addToMatrixAndRHS(mat, rhs, EM_S, EM_F, add_EM_S, add_EM_F, firstNode);
# Line 3537  void Rectangle::assemblePDEBoundarySyste Line 3439  void Rectangle::assemblePDEBoundarySyste
3439          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
3440          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
3441      }      }
3442      const double w0 = 0.31100423396407310779*m_dx[1];      const double SQRT3 = 1.73205080756887719318;
3443      const double w1 = 0.022329099369260225539*m_dx[1];      const double w5 = m_dx[0]/12;
3444      const double w10 = 0.022329099369260225539*m_dx[0];      const double w6 = w5*(SQRT3 + 2);
3445      const double w11 = 0.16666666666666666667*m_dx[0];      const double w7 = w5*(-SQRT3 + 2);
3446      const double w12 = 0.33333333333333333333*m_dx[0];      const double w8 = w5*(SQRT3 + 3);
3447      const double w13 = 0.39433756729740644113*m_dx[0];      const double w9 = w5*(-SQRT3 + 3);
3448      const double w14 = 0.10566243270259355887*m_dx[0];      const double w2 = m_dx[1]/12;
3449      const double w15 = 0.5*m_dx[0];      const double w0 = w2*(SQRT3 + 2);
3450      const double w2 = 0.083333333333333333333*m_dx[1];      const double w1 = w2*(-SQRT3 + 2);
3451      const double w3 = 0.33333333333333333333*m_dx[1];      const double w3 = w2*(SQRT3 + 3);
3452      const double w4 = 0.16666666666666666667*m_dx[1];      const double w4 = w2*(-SQRT3 + 3);
     const double w5 = 0.39433756729740644113*m_dx[1];  
     const double w6 = 0.10566243270259355887*m_dx[1];  
     const double w7 = 0.5*m_dx[1];  
     const double w8 = 0.083333333333333333333*m_dx[0];  
     const double w9 = 0.31100423396407310779*m_dx[0];  
3453      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3454      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3455      rhs.requireWrite();      rhs.requireWrite();
# Line 3573  void Rectangle::assemblePDEBoundarySyste Line 3470  void Rectangle::assemblePDEBoundarySyste
3470                          if (d.actsExpanded()) {                          if (d.actsExpanded()) {
3471                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3472                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3473                                      const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];                                      const double d_0 = d_p[INDEX3(k,m,0,numEq,numComp)];
3474                                      const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];                                      const double d_1 = d_p[INDEX3(k,m,1,numEq,numComp)];
3475                                      const double tmp0_0 = d_0 + d_1;                                      const double tmp0 = w2*(d_0 + d_1);
3476                                      const double tmp0_1 = d_0*w0;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=d_0*w0 + d_1*w1;
3477                                      const double tmp1_1 = d_1*w1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0;
3478                                      const double tmp2_1 = tmp0_0*w2;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0;
3479                                      const double tmp3_1 = d_1*w0;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=d_0*w1 + d_1*w0;
                                     const double tmp4_1 = d_0*w1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
3480                                  }                                  }
3481                               }                               }
3482                          } else { /* constant data */                          } else { // constant data
3483                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3484                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3485                                      const double d_0 = d_p[INDEX2(k, m, numEq)];                                      const double d_0 = d_p[INDEX2(k, m, numEq)];
3486                                      const double tmp0_1 = d_0*w3;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=4*d_0*w2;
3487                                      const double tmp1_1 = d_0*w4;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=2*d_0*w2;
3488                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=2*d_0*w2;
3489                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=4*d_0*w2;
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;  
3490                                  }                                  }
3491                              }                              }
3492                          }                          }
# Line 3610  void Rectangle::assemblePDEBoundarySyste Line 3500  void Rectangle::assemblePDEBoundarySyste
3500                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3501                                  const double y_0 = y_p[INDEX2(k, 0, numEq)];                                  const double y_0 = y_p[INDEX2(k, 0, numEq)];
3502                                  const double y_1 = y_p[INDEX2(k, 1, numEq)];                                  const double y_1 = y_p[INDEX2(k, 1, numEq)];
3503                                  const double tmp3_1 = w5*y_1;                                  EM_F[INDEX2(k,0,numEq)]+=w3*y_0 + w4*y_1;
3504                                  const double tmp2_1 = w6*y_0;                                  EM_F[INDEX2(k,2,numEq)]+=w3*y_1 + w4*y_0;
3505                                  const double tmp0_1 = w6*y_1;                              }
3506                                  const double tmp1_1 = w5*y_0;                          } else { // constant data
3507                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;                              for (index_t k=0; k<numEq; k++) {
3508                                  EM_F[INDEX2(k,2,numEq)]+=tmp2_1 + tmp3_1;                                  EM_F[INDEX2(k,0,numEq)]+=6*w2*y_p[k];
3509                              }                                  EM_F[INDEX2(k,2,numEq)]+=6*w2*y_p[k];
                         } else { /* constant data */  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double y_0 = y_p[k];  
                                 const double tmp0_1 = w7*y_0;  
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp0_1;  
3510                              }                              }
3511                          }                          }
3512                      }                      }
# Line 3648  void Rectangle::assemblePDEBoundarySyste Line 3532  void Rectangle::assemblePDEBoundarySyste
3532                          if (d.actsExpanded()) {                          if (d.actsExpanded()) {
3533                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3534                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3535                                      const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];                                      const double d_0 = d_p[INDEX3(k,m,0,numEq,numComp)];
3536                                      const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];                                      const double d_1 = d_p[INDEX3(k,m,1,numEq,numComp)];
3537                                      const double tmp0_0 = d_0 + d_1;                                      const double tmp0 = w2*(d_0 + d_1);
3538                                      const double tmp4_1 = d_1*w1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=d_0*w0 + d_1*w1;
3539                                      const double tmp1_1 = d_0*w1;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0;
3540                                      const double tmp3_1 = d_0*w0;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0;
3541                                      const double tmp0_1 = d_1*w0;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=d_0*w1 + d_1*w0;
                                     const double tmp2_1 = tmp0_0*w2;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
3542                                  }                                  }
3543                               }                               }
3544                          } else { /* constant data */                          } else { // constant data
3545                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3546                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3547                                      const double d_0 = d_p[INDEX2(k, m, numEq)];                                      const double d_0 = d_p[INDEX2(k, m, numEq)];
3548                                      const double tmp1_1 = d_0*w4;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=4*d_0*w2;
3549                                      const double tmp0_1 = d_0*w3;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=2*d_0*w2;
3550                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=2*d_0*w2;
3551                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=4*d_0*w2;
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;  
3552                                  }                                  }
3553                              }                              }
3554                          }                          }
# Line 3685  void Rectangle::assemblePDEBoundarySyste Line 3562  void Rectangle::assemblePDEBoundarySyste
3562                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3563                                  const double y_0 = y_p[INDEX2(k, 0, numEq)];                                  const double y_0 = y_p[INDEX2(k, 0, numEq)];
3564                                  const double y_1 = y_p[INDEX2(k, 1, numEq)];                                  const double y_1 = y_p[INDEX2(k, 1, numEq)];
3565                                  const double tmp3_1 = w5*y_1;                                  EM_F[INDEX2(k,1,numEq)]+=w3*y_0 + w4*y_1;
3566                                  const double tmp2_1 = w6*y_0;                                  EM_F[INDEX2(k,3,numEq)]+=w3*y_1 + w4*y_0;
3567                                  const double tmp0_1 = w6*y_1;                              }
3568                                  const double tmp1_1 = w5*y_0;                          } else { // constant data
3569                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp1_1;                              for (index_t k=0; k<numEq; k++) {
3570                                  EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;                                  EM_F[INDEX2(k,1,numEq)]+=6*w2*y_p[k];
3571                              }                                  EM_F[INDEX2(k,3,numEq)]+=6*w2*y_p[k];
                         } else { /* constant data */  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double y_0 = y_p[k];  
                                 const double tmp0_1 = w7*y_0;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp0_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp0_1;  
3572                              }                              }
3573                          }                          }
3574                      }                      }
# Line 3723  void Rectangle::assemblePDEBoundarySyste Line 3594  void Rectangle::assemblePDEBoundarySyste
3594                          if (d.actsExpanded()) {                          if (d.actsExpanded()) {
3595                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3596                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3597                                      const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];                                      const double d_0 = d_p[INDEX3(k,m,0,numEq,numComp)];
3598                                      const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];                                      const double d_1 = d_p[INDEX3(k,m,1,numEq,numComp)];
3599                                      const double tmp0_0 = d_0 + d_1;                                      const double tmp0 = w5*(d_0 + d_1);
3600                                      const double tmp4_1 = d_1*w9;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=d_0*w6 + d_1*w7;
3601                                      const double tmp2_1 = d_0*w9;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0;
3602                                      const double tmp0_1 = tmp0_0*w8;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0;
3603                                      const double tmp1_1 = d_1*w10;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=d_0*w7 + d_1*w6;
                                     const double tmp3_1 = d_0*w10;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
3604                                  }                                  }
3605                               }                               }
3606                          } else { /* constant data */                          } else { // constant data
3607                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3608                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3609                                      const double d_0 = d_p[INDEX2(k, m, numEq)];                                      const double d_0 = d_p[INDEX2(k, m, numEq)];
3610                                      const double tmp0_1 = d_0*w11;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=4*d_0*w5;
3611                                      const double tmp1_1 = d_0*w12;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=2*d_0*w5;
3612                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=2*d_0*w5;
3613                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=4*d_0*w5;
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp1_1;  
3614                                  }                                  }
3615                              }                              }
3616                          }                          }
# Line 3760  void Rectangle::assemblePDEBoundarySyste Line 3624  void Rectangle::assemblePDEBoundarySyste
3624                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3625                                  const double y_0 = y_p[INDEX2(k, 0, numEq)];                                  const double y_0 = y_p[INDEX2(k, 0, numEq)];
3626                                  const double y_1 = y_p[INDEX2(k, 1, numEq)];                                  const double y_1 = y_p[INDEX2(k, 1, numEq)];
3627                                  const double tmp2_1 = w13*y_1;                                  EM_F[INDEX2(k,0,numEq)]+=w8*y_0 + w9*y_1;
3628                                  const double tmp1_1 = w14*y_1;                                  EM_F[INDEX2(k,1,numEq)]+=w8*y_1 + w9*y_0;
3629                                  const double tmp3_1 = w14*y_0;                              }
3630                                  const double tmp0_1 = w13*y_0;                          } else { // constant data
3631                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;                              for (index_t k=0; k<numEq; k++) {
3632                                  EM_F[INDEX2(k,1,numEq)]+=tmp2_1 + tmp3_1;                                  EM_F[INDEX2(k,0,numEq)]+=6*w5*y_p[k];
3633                              }                                  EM_F[INDEX2(k,1,numEq)]+=6*w5*y_p[k];
                         } else { /* constant data */  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double y_0 = y_p[k];  
                                 const double tmp0_1 = w15*y_0;  
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp0_1;  
3634                              }                              }
3635                          }                          }
3636                      }                      }
# Line 3798  void Rectangle::assemblePDEBoundarySyste Line 3656  void Rectangle::assemblePDEBoundarySyste
3656                          if (d.actsExpanded()) {                          if (d.actsExpanded()) {
3657                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3658                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3659                                      const double d_0 = d_p[INDEX3(k, m, 0, numEq, numComp)];                                      const double d_0 = d_p[INDEX3(k,m,0,numEq,numComp)];
3660                                      const double d_1 = d_p[INDEX3(k, m, 1, numEq, numComp)];                                      const double d_1 = d_p[INDEX3(k,m,1,numEq,numComp)];
3661                                      const double tmp0_0 = d_0 + d_1;                                      const double tmp0 = w5*(d_0 + d_1);
3662                                      const double tmp2_1 = d_1*w9;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=d_0*w6 + d_1*w7;
3663                                      const double tmp4_1 = d_0*w9;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0;
3664                                      const double tmp0_1 = tmp0_0*w8;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0;
3665                                      const double tmp3_1 = d_1*w10;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=d_0*w7 + d_1*w6;
                                     const double tmp1_1 = d_0*w10;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
3666                                  }                                  }
3667                               }                               }
3668                          } else { /* constant data */                          } else { // constant data
3669                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3670                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
3671                                      const double d_0 = d_p[INDEX2(k, m, numEq)];                                      const double d_0 = d_p[INDEX2(k, m, numEq)];
3672                                      const double tmp0_1 = d_0*w11;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=4*d_0*w5;
3673                                      const double tmp1_1 = d_0*w12;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=2*d_0*w5;
3674                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=2*d_0*w5;
3675                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=4*d_0*w5;
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp1_1;  
3676                                  }                                  }
3677                              }                              }
3678                          }                          }
# Line 3835  void Rectangle::assemblePDEBoundarySyste Line 3686  void Rectangle::assemblePDEBoundarySyste
3686                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
3687                                  const double y_0 = y_p[INDEX2(k, 0, numEq)];                                  const double y_0 = y_p[INDEX2(k, 0, numEq)];
3688                                  const double y_1 = y_p[INDEX2(k, 1, numEq)];                                  const double y_1 = y_p[INDEX2(k, 1, numEq)];
3689                                  const double tmp2_1 = w13*y_1;                                  EM_F[INDEX2(k,2,numEq)]+=w8*y_0 + w9*y_1;
3690                                  const double tmp1_1 = w14*y_1;                                  EM_F[INDEX2(k,3,numEq)]+=w8*y_1 + w9*y_0;
3691                                  const double tmp3_1 = w14*y_0;                              }
3692                                  const double tmp0_1 = w13*y_0;                          } else { // constant data
3693                                  EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp1_1;                              for (index_t k=0; k<numEq; k++) {
3694                                  EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;                                  EM_F[INDEX2(k,2,numEq)]+=6*w5*y_p[k];
3695                              }                                  EM_F[INDEX2(k,3,numEq)]+=6*w5*y_p[k];
                         } else { /* constant data */  
                             for (index_t k=0; k<numEq; k++) {  
                                 const double y_0 = y_p[k];  
                                 const double tmp0_1 = w15*y_0;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp0_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp0_1;  
3696                              }                              }
3697                          }                          }
3698                      }                      }
# Line 3871  void Rectangle::assemblePDEBoundarySyste Line 3716  void Rectangle::assemblePDEBoundarySyste
3716          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
3717          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
3718      }      }
3719      const double w0 = 0.25*m_dx[0];      const double w0 = m_dx[0]/4;
3720      const double w1 = 0.25*m_dx[1];      const double w1 = m_dx[1]/4;
3721      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3722      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3723      rhs.requireWrite();      rhs.requireWrite();
# Line 3905  void Rectangle::assemblePDEBoundarySyste Line 3750  void Rectangle::assemblePDEBoundarySyste
3750                      if (add_EM_F) {                      if (add_EM_F) {
3751                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(k1);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(k1);
3752                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3753                              const double tmp0 = 2*w1*y_p[k];                              EM_F[INDEX2(k,0,numEq)]+=2*w1*y_p[k];
3754                              EM_F[INDEX2(k,0,numEq)]+=tmp0;                              EM_F[INDEX2(k,2,numEq)]+=2*w1*y_p[k];
                             EM_F[INDEX2(k,2,numEq)]+=tmp0;  
3755                          }                          }
3756                      }                      }
3757                      const index_t firstNode=m_NN[0]*k1;                      const index_t firstNode=m_NN[0]*k1;
# Line 3945  void Rectangle::assemblePDEBoundarySyste Line 3789  void Rectangle::assemblePDEBoundarySyste
3789                      if (add_EM_F) {                      if (add_EM_F) {
3790                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3791                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3792                              const double tmp0 = 2*w1*y_p[k];                              EM_F[INDEX2(k,1,numEq)]+=2*w1*y_p[k];
3793                              EM_F[INDEX2(k,1,numEq)]+=tmp0;                              EM_F[INDEX2(k,3,numEq)]+=2*w1*y_p[k];
                             EM_F[INDEX2(k,3,numEq)]+=tmp0;  
3794                          }                          }
3795                      }                      }
3796                      const index_t firstNode=m_NN[0]*(k1+1)-2;                      const index_t firstNode=m_NN[0]*(k1+1)-2;
# Line 3985  void Rectangle::assemblePDEBoundarySyste Line 3828  void Rectangle::assemblePDEBoundarySyste
3828                      if (add_EM_F) {                      if (add_EM_F) {
3829                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3830                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3831                              const double tmp0 = 2*w0*y_p[k];                              EM_F[INDEX2(k,0,numEq)]+=2*w0*y_p[k];
3832                              EM_F[INDEX2(k,0,numEq)]+=tmp0;                              EM_F[INDEX2(k,1,numEq)]+=2*w0*y_p[k];
                             EM_F[INDEX2(k,1,numEq)]+=tmp0;  
3833                          }                          }
3834                      }                      }
3835                      const index_t firstNode=k0;                      const index_t firstNode=k0;
# Line 4025  void Rectangle::assemblePDEBoundarySyste Line 3867  void Rectangle::assemblePDEBoundarySyste
3867                      if (add_EM_F) {                      if (add_EM_F) {
3868                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3869                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3870                              const double tmp0 = 2*w0*y_p[k];                              EM_F[INDEX2(k,2,numEq)]+=2*w0*y_p[k];
3871                              EM_F[INDEX2(k,2,numEq)]+=tmp0;                              EM_F[INDEX2(k,3,numEq)]+=2*w0*y_p[k];
                             EM_F[INDEX2(k,3,numEq)]+=tmp0;  
3872                          }                          }
3873                      }                      }
3874                      const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;                      const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;

Legend:
Removed from v.4375  
changed lines
  Added in v.4378

  ViewVC Help
Powered by ViewVC 1.1.26