/[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 4368 by jfenwick, Thu Apr 18 02:26:35 2013 UTC revision 4370 by caltinay, Fri Apr 19 06:15:24 2013 UTC
# Line 1819  void Rectangle::assemblePDESingle(Paso_S Line 1819  void Rectangle::assemblePDESingle(Paso_S
1819          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
1820          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
1821  {  {
1822        /* GENERATOR SNIP_PDE_SINGLE_PRE TOP */
1823      const double w0 = -0.1555021169820365539*m_dx[1]/m_dx[0];      const double w0 = -0.1555021169820365539*m_dx[1]/m_dx[0];
1824      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
1825      const double w2 = -0.15550211698203655390;      const double w2 = -0.15550211698203655390;
1826      const double w3 = 0.041666666666666666667*m_dx[0]/m_dx[1];      const double w3 = 0.041666666666666666667*m_dx[0]/m_dx[1];
1827      const double w4 = 0.15550211698203655390;      const double w4 = -0.01116454968463011277*m_dx[1]/m_dx[0];
1828      const double w5 = -0.041666666666666666667;      const double w5 = 0.011164549684630112770;
1829      const double w6 = -0.01116454968463011277*m_dx[1]/m_dx[0];      const double w6 = -0.041666666666666666667*m_dx[1]/m_dx[0];
1830      const double w7 = 0.011164549684630112770;      const double w7 = 0.1555021169820365539*m_dx[0]/m_dx[1];
1831      const double w8 = -0.011164549684630112770;      const double w8 = 0.01116454968463011277*m_dx[0]/m_dx[1];
1832      const double w9 = -0.041666666666666666667*m_dx[1]/m_dx[0];      const double w9 = -0.25000000000000000000;
1833      const double w10 = -0.041666666666666666667*m_dx[0]/m_dx[1];      const double w10 = -0.16666666666666666667*m_dx[0]/m_dx[1];
1834      const double w11 = 0.1555021169820365539*m_dx[1]/m_dx[0];      const double w11 = -0.032861463941450536761*m_dx[1];
1835      const double w12 = 0.1555021169820365539*m_dx[0]/m_dx[1];      const double w12 = -0.032861463941450536761*m_dx[0];
1836      const double w13 = 0.01116454968463011277*m_dx[0]/m_dx[1];      const double w13 = -0.12264065304058601714*m_dx[1];
1837      const double w14 = 0.01116454968463011277*m_dx[1]/m_dx[0];      const double w14 = -0.0023593469594139828636*m_dx[1];
1838      const double w15 = 0.041666666666666666667*m_dx[1]/m_dx[0];      const double w15 = -0.008805202725216129906*m_dx[0];
1839      const double w16 = -0.01116454968463011277*m_dx[0]/m_dx[1];      const double w16 = -0.008805202725216129906*m_dx[1];
1840      const double w17 = -0.1555021169820365539*m_dx[0]/m_dx[1];      const double w17 = -0.12264065304058601714*m_dx[0];
1841      const double w18 = -0.33333333333333333333*m_dx[1]/m_dx[0];      const double w18 = -0.0023593469594139828636*m_dx[0];
1842      const double w19 = 0.25;      const double w19 = -0.16666666666666666667*m_dx[1];
1843      const double w20 = -0.25;      const double w20 = -0.083333333333333333333*m_dx[0];
1844      const double w21 = 0.16666666666666666667*m_dx[0]/m_dx[1];      const double w21 = 0.025917019497006092316*m_dx[0]*m_dx[1];
1845      const double w22 = -0.16666666666666666667*m_dx[1]/m_dx[0];      const double w22 = 0.0018607582807716854616*m_dx[0]*m_dx[1];
1846      const double w23 = -0.16666666666666666667*m_dx[0]/m_dx[1];      const double w23 = 0.0069444444444444444444*m_dx[0]*m_dx[1];
1847      const double w24 = 0.33333333333333333333*m_dx[1]/m_dx[0];      const double w24 = 0.09672363354357992482*m_dx[0]*m_dx[1];
1848      const double w25 = 0.33333333333333333333*m_dx[0]/m_dx[1];      const double w25 = 0.000498588678642297402*m_dx[0]*m_dx[1];
1849      const double w26 = 0.16666666666666666667*m_dx[1]/m_dx[0];      const double w26 = 0.19716878364870322056*m_dx[1];
1850      const double w27 = -0.33333333333333333333*m_dx[0]/m_dx[1];      const double w27 = 0.052831216351296779436*m_dx[1];
1851      const double w28 = -0.032861463941450536761*m_dx[1];      const double w28 = 0.19716878364870322056*m_dx[0];
1852      const double w29 = -0.032861463941450536761*m_dx[0];      const double w29 = 0.052831216351296779436*m_dx[0];
1853      const double w30 = -0.12264065304058601714*m_dx[1];      /* GENERATOR SNIP_PDE_SINGLE_PRE BOTTOM */
     const double w31 = -0.0023593469594139828636*m_dx[1];  
     const double w32 = -0.008805202725216129906*m_dx[0];  
     const double w33 = -0.008805202725216129906*m_dx[1];  
     const double w34 = 0.032861463941450536761*m_dx[1];  
     const double w35 = 0.008805202725216129906*m_dx[1];  
     const double w36 = 0.008805202725216129906*m_dx[0];  
     const double w37 = 0.0023593469594139828636*m_dx[1];  
     const double w38 = 0.12264065304058601714*m_dx[1];  
     const double w39 = 0.032861463941450536761*m_dx[0];  
     const double w40 = -0.12264065304058601714*m_dx[0];  
     const double w41 = -0.0023593469594139828636*m_dx[0];  
     const double w42 = 0.0023593469594139828636*m_dx[0];  
     const double w43 = 0.12264065304058601714*m_dx[0];  
     const double w44 = -0.16666666666666666667*m_dx[1];  
     const double w45 = -0.083333333333333333333*m_dx[0];  
     const double w46 = 0.083333333333333333333*m_dx[1];  
     const double w47 = 0.16666666666666666667*m_dx[1];  
     const double w48 = 0.083333333333333333333*m_dx[0];  
     const double w49 = -0.16666666666666666667*m_dx[0];  
     const double w50 = 0.16666666666666666667*m_dx[0];  
     const double w51 = -0.083333333333333333333*m_dx[1];  
     const double w52 = 0.025917019497006092316*m_dx[0]*m_dx[1];  
     const double w53 = 0.0018607582807716854616*m_dx[0]*m_dx[1];  
     const double w54 = 0.0069444444444444444444*m_dx[0]*m_dx[1];  
     const double w55 = 0.09672363354357992482*m_dx[0]*m_dx[1];  
     const double w56 = 0.00049858867864229740201*m_dx[0]*m_dx[1];  
     const double w57 = 0.055555555555555555556*m_dx[0]*m_dx[1];  
     const double w58 = 0.027777777777777777778*m_dx[0]*m_dx[1];  
     const double w59 = 0.11111111111111111111*m_dx[0]*m_dx[1];  
     const double w60 = -0.19716878364870322056*m_dx[1];  
     const double w61 = -0.19716878364870322056*m_dx[0];  
     const double w62 = -0.052831216351296779436*m_dx[0];  
     const double w63 = -0.052831216351296779436*m_dx[1];  
     const double w64 = 0.19716878364870322056*m_dx[1];  
     const double w65 = 0.052831216351296779436*m_dx[1];  
     const double w66 = 0.19716878364870322056*m_dx[0];  
     const double w67 = 0.052831216351296779436*m_dx[0];  
     const double w68 = -0.5*m_dx[1];  
     const double w69 = -0.5*m_dx[0];  
     const double w70 = 0.5*m_dx[1];  
     const double w71 = 0.5*m_dx[0];  
     const double w72 = 0.1555021169820365539*m_dx[0]*m_dx[1];  
     const double w73 = 0.041666666666666666667*m_dx[0]*m_dx[1];  
     const double w74 = 0.01116454968463011277*m_dx[0]*m_dx[1];  
     const double w75 = 0.25*m_dx[0]*m_dx[1];  
1854    
1855      rhs.requireWrite();      rhs.requireWrite();
1856  #pragma omp parallel  #pragma omp parallel
# Line 1908  void Rectangle::assemblePDESingle(Paso_S Line 1864  void Rectangle::assemblePDESingle(Paso_S
1864                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
1865                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
1866                      const index_t e = k0 + m_NE[0]*k1;                      const index_t e = k0 + m_NE[0]*k1;
1867                        /* GENERATOR SNIP_PDE_SINGLE TOP */
1868                      ///////////////                      ///////////////
1869                      // process A //                      // process A //
1870                      ///////////////                      ///////////////
1871                      if (!A.isEmpty()) {                      if (!A.isEmpty()) {
1872                          add_EM_S=true;                          add_EM_S = true;
1873                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p = const_cast<escript::Data*>(&A)->getSampleDataRO(e);
1874                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
1875                              const double A_00_0 = A_p[INDEX3(0,0,0,2,2)];                              const double A_00_0 = A_p[INDEX3(0,0,0,2,2)];
                             const double A_10_0 = A_p[INDEX3(1,0,0,2,2)];  
1876                              const double A_01_0 = A_p[INDEX3(0,1,0,2,2)];                              const double A_01_0 = A_p[INDEX3(0,1,0,2,2)];
1877                                const double A_10_0 = A_p[INDEX3(1,0,0,2,2)];
1878                              const double A_11_0 = A_p[INDEX3(1,1,0,2,2)];                              const double A_11_0 = A_p[INDEX3(1,1,0,2,2)];
1879                              const double A_00_1 = A_p[INDEX3(0,0,1,2,2)];                              const double A_00_1 = A_p[INDEX3(0,0,1,2,2)];
                             const double A_10_1 = A_p[INDEX3(1,0,1,2,2)];  
1880                              const double A_01_1 = A_p[INDEX3(0,1,1,2,2)];                              const double A_01_1 = A_p[INDEX3(0,1,1,2,2)];
1881                                const double A_10_1 = A_p[INDEX3(1,0,1,2,2)];
1882                              const double A_11_1 = A_p[INDEX3(1,1,1,2,2)];                              const double A_11_1 = A_p[INDEX3(1,1,1,2,2)];
1883                              const double A_00_2 = A_p[INDEX3(0,0,2,2,2)];                              const double A_00_2 = A_p[INDEX3(0,0,2,2,2)];
                             const double A_10_2 = A_p[INDEX3(1,0,2,2,2)];  
1884                              const double A_01_2 = A_p[INDEX3(0,1,2,2,2)];                              const double A_01_2 = A_p[INDEX3(0,1,2,2,2)];
1885                                const double A_10_2 = A_p[INDEX3(1,0,2,2,2)];
1886                              const double A_11_2 = A_p[INDEX3(1,1,2,2,2)];                              const double A_11_2 = A_p[INDEX3(1,1,2,2,2)];
1887                              const double A_00_3 = A_p[INDEX3(0,0,3,2,2)];                              const double A_00_3 = A_p[INDEX3(0,0,3,2,2)];
                             const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];  
1888                              const double A_01_3 = A_p[INDEX3(0,1,3,2,2)];                              const double A_01_3 = A_p[INDEX3(0,1,3,2,2)];
1889                                const double A_10_3 = A_p[INDEX3(1,0,3,2,2)];
1890                              const double A_11_3 = A_p[INDEX3(1,1,3,2,2)];                              const double A_11_3 = A_p[INDEX3(1,1,3,2,2)];
1891                              const double tmp0_0 = A_01_0 + A_01_3;                              const double tmp0 = w3*(A_11_0 + A_11_1 + A_11_2 + A_11_3);
1892                              const double tmp1_0 = A_00_0 + A_00_1;                              const double tmp1 = w1*(A_01_0 + A_01_3 - A_10_1 - A_10_2);
1893                              const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;                              const double tmp2 = w4*(A_00_2 + A_00_3);
1894                              const double tmp3_0 = A_00_2 + A_00_3;                              const double tmp3 = w0*(A_00_0 + A_00_1);
1895                              const double tmp4_0 = A_10_1 + A_10_2;                              const double tmp4 = w5*(A_01_2 - A_10_3);
1896                              const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;                              const double tmp5 = w2*(-A_01_1 + A_10_0);
1897                              const double tmp6_0 = A_01_3 + A_10_0;                              const double tmp6 = w5*(A_01_3 + A_10_0);
1898                              const double tmp7_0 = A_01_0 + A_10_3;                              const double tmp7 = w3*(-A_11_0 - A_11_1 - A_11_2 - A_11_3);
1899                              const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;                              const double tmp8 = w6*(A_00_0 + A_00_1 + A_00_2 + A_00_3);
1900                              const double tmp9_0 = A_01_0 + A_10_0;                              const double tmp9 = w1*(A_01_1 + A_01_2 + A_10_1 + A_10_2);
1901                              const double tmp12_0 = A_11_0 + A_11_2;                              const double tmp10 = w2*(-A_01_0 - A_10_3);
1902                              const double tmp10_0 = A_01_3 + A_10_3;                              const double tmp11 = w4*(A_00_0 + A_00_1);
1903                              const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;                              const double tmp12 = w0*(A_00_2 + A_00_3);
1904                              const double tmp13_0 = A_01_2 + A_10_1;                              const double tmp13 = w5*(A_01_1 - A_10_0);
1905                              const double tmp11_0 = A_11_1 + A_11_3;                              const double tmp14 = w2*(-A_01_2 + A_10_3);
1906                              const double tmp18_0 = A_01_1 + A_10_1;                              const double tmp15 = w7*(A_11_0 + A_11_2);
1907                              const double tmp15_0 = A_01_1 + A_10_2;                              const double tmp16 = w4*(-A_00_2 - A_00_3);
1908                              const double tmp16_0 = A_10_0 + A_10_3;                              const double tmp17 = w0*(-A_00_0 - A_00_1);
1909                              const double tmp17_0 = A_01_1 + A_01_2;                              const double tmp18 = w5*(A_01_3 + A_10_3);
1910                              const double tmp19_0 = A_01_2 + A_10_2;                              const double tmp19 = w8*(A_11_1 + A_11_3);
1911                              const double tmp0_1 = A_10_3*w8;                              const double tmp20 = w2*(-A_01_0 - A_10_0);
1912                              const double tmp1_1 = tmp0_0*w1;                              const double tmp21 = w7*(A_11_1 + A_11_3);
1913                              const double tmp2_1 = A_01_1*w4;                              const double tmp22 = w4*(-A_00_0 - A_00_1);
1914                              const double tmp3_1 = tmp1_0*w0;                              const double tmp23 = w0*(-A_00_2 - A_00_3);
1915                              const double tmp4_1 = A_01_2*w7;                              const double tmp24 = w5*(A_01_0 + A_10_0);
1916                              const double tmp5_1 = tmp2_0*w3;                              const double tmp25 = w8*(A_11_0 + A_11_2);
1917                              const double tmp6_1 = tmp3_0*w6;                              const double tmp26 = w2*(-A_01_3 - A_10_3);
1918                              const double tmp7_1 = A_10_0*w2;                              const double tmp27 = w5*(-A_01_1 - A_10_2);
1919                              const double tmp8_1 = tmp4_0*w5;                              const double tmp28 = w1*(-A_01_0 - A_01_3 - A_10_0 - A_10_3);
1920                              const double tmp9_1 = tmp2_0*w10;                              const double tmp29 = w2*(A_01_2 + A_10_1);
1921                              const double tmp14_1 = A_10_0*w8;                              const double tmp30 = w7*(-A_11_1 - A_11_3);
1922                              const double tmp23_1 = tmp3_0*w14;                              const double tmp31 = w1*(-A_01_1 - A_01_2 + A_10_0 + A_10_3);
1923                              const double tmp35_1 = A_01_0*w8;                              const double tmp32 = w5*(-A_01_0 + A_10_2);
1924                              const double tmp54_1 = tmp13_0*w8;                              const double tmp33 = w8*(-A_11_0 - A_11_2);
1925                              const double tmp20_1 = tmp9_0*w4;                              const double tmp34 = w6*(-A_00_0 - A_00_1 - A_00_2 - A_00_3);
1926                              const double tmp25_1 = tmp12_0*w12;                              const double tmp35 = w2*(A_01_3 - A_10_1);
1927                              const double tmp44_1 = tmp7_0*w7;                              const double tmp36 = w5*(A_01_0 + A_10_3);
1928                              const double tmp26_1 = tmp10_0*w4;                              const double tmp37 = w2*(-A_01_3 - A_10_0);
1929                              const double tmp52_1 = tmp18_0*w8;                              const double tmp38 = w7*(-A_11_0 - A_11_2);
1930                              const double tmp48_1 = A_10_1*w7;                              const double tmp39 = w5*(-A_01_3 + A_10_1);
1931                              const double tmp46_1 = A_01_3*w8;                              const double tmp40 = w8*(-A_11_1 - A_11_3);
1932                              const double tmp50_1 = A_01_0*w2;                              const double tmp41 = w2*(A_01_0 - A_10_2);
1933                              const double tmp56_1 = tmp19_0*w8;                              const double tmp42 = w5*(A_01_1 - A_10_3);
1934                              const double tmp19_1 = A_10_3*w2;                              const double tmp43 = w2*(-A_01_2 + A_10_0);
1935                              const double tmp47_1 = A_10_2*w4;                              const double tmp44 = w5*(A_01_2 - A_10_0);
1936                              const double tmp16_1 = tmp3_0*w0;                              const double tmp45 = w2*(-A_01_1 + A_10_3);
1937                              const double tmp18_1 = tmp1_0*w6;                              const double tmp46 = w5*(-A_01_0 + A_10_1);
1938                              const double tmp31_1 = tmp11_0*w12;                              const double tmp47 = w2*(A_01_3 - A_10_2);
1939                              const double tmp55_1 = tmp15_0*w2;                              const double tmp48 = w5*(-A_01_1 - A_10_1);
1940                              const double tmp39_1 = A_10_2*w7;                              const double tmp49 = w2*(A_01_2 + A_10_2);
1941                              const double tmp11_1 = tmp6_0*w7;                              const double tmp50 = w5*(-A_01_3 + A_10_2);
1942                              const double tmp40_1 = tmp11_0*w17;                              const double tmp51 = w2*(A_01_0 - A_10_1);
1943                              const double tmp34_1 = tmp15_0*w8;                              const double tmp52 = w5*(-A_01_2 - A_10_1);
1944                              const double tmp33_1 = tmp14_0*w5;                              const double tmp53 = w2*(A_01_1 + A_10_2);
1945                              const double tmp24_1 = tmp11_0*w13;                              const double tmp54 = w5*(-A_01_2 - A_10_2);
1946                              const double tmp43_1 = tmp17_0*w5;                              const double tmp55 = w2*(A_01_1 + A_10_1);
1947                              const double tmp15_1 = A_01_2*w4;                              EM_S[INDEX2(0,0,4)]+=tmp15 + tmp16 + tmp17 + tmp18 + tmp19 + tmp20 + tmp9;
1948                              const double tmp53_1 = tmp19_0*w2;                              EM_S[INDEX2(0,1,4)]+=tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5;
1949                              const double tmp27_1 = tmp3_0*w11;                              EM_S[INDEX2(0,2,4)]+=tmp31 + tmp34 + tmp38 + tmp39 + tmp40 + tmp41;
1950                              const double tmp32_1 = tmp13_0*w2;                              EM_S[INDEX2(0,3,4)]+=tmp28 + tmp52 + tmp53 + tmp7 + tmp8;
1951                              const double tmp10_1 = tmp5_0*w9;                              EM_S[INDEX2(1,0,4)]+=tmp0 + tmp2 + tmp3 + tmp31 + tmp50 + tmp51;
1952                              const double tmp37_1 = A_10_1*w4;                              EM_S[INDEX2(1,1,4)]+=tmp16 + tmp17 + tmp21 + tmp25 + tmp28 + tmp54 + tmp55;
1953                              const double tmp38_1 = tmp5_0*w15;                              EM_S[INDEX2(1,2,4)]+=tmp10 + tmp6 + tmp7 + tmp8 + tmp9;
1954                              const double tmp17_1 = A_01_1*w7;                              EM_S[INDEX2(1,3,4)]+=tmp1 + tmp30 + tmp33 + tmp34 + tmp44 + tmp45;
1955                              const double tmp12_1 = tmp7_0*w4;                              EM_S[INDEX2(2,0,4)]+=tmp1 + tmp34 + tmp38 + tmp40 + tmp42 + tmp43;
1956                              const double tmp22_1 = tmp10_0*w7;                              EM_S[INDEX2(2,1,4)]+=tmp36 + tmp37 + tmp7 + tmp8 + tmp9;
1957                              const double tmp57_1 = tmp18_0*w2;                              EM_S[INDEX2(2,2,4)]+=tmp15 + tmp19 + tmp22 + tmp23 + tmp28 + tmp48 + tmp49;
1958                              const double tmp28_1 = tmp9_0*w7;                              EM_S[INDEX2(2,3,4)]+=tmp0 + tmp11 + tmp12 + tmp31 + tmp46 + tmp47;
1959                              const double tmp29_1 = tmp1_0*w14;                              EM_S[INDEX2(3,0,4)]+=tmp27 + tmp28 + tmp29 + tmp7 + tmp8;
1960                              const double tmp51_1 = tmp11_0*w16;                              EM_S[INDEX2(3,1,4)]+=tmp30 + tmp31 + tmp32 + tmp33 + tmp34 + tmp35;
1961                              const double tmp42_1 = tmp12_0*w16;                              EM_S[INDEX2(3,2,4)]+=tmp0 + tmp1 + tmp11 + tmp12 + tmp13 + tmp14;
1962                              const double tmp49_1 = tmp12_0*w17;                              EM_S[INDEX2(3,3,4)]+=tmp21 + tmp22 + tmp23 + tmp24 + tmp25 + tmp26 + tmp9;
                             const double tmp21_1 = tmp1_0*w11;  
                             const double tmp45_1 = tmp6_0*w4;  
                             const double tmp13_1 = tmp8_0*w1;  
                             const double tmp36_1 = tmp16_0*w1;  
                             const double tmp41_1 = A_01_3*w2;  
                             const double tmp30_1 = tmp12_0*w13;  
                             EM_S[INDEX2(0,0,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
1963                          } else { // constant data                          } else { // constant data
1964                              const double A_00 = A_p[INDEX2(0,0,2)];                              const double A_00 = A_p[INDEX2(0,0,2)];
                             const double A_10 = A_p[INDEX2(1,0,2)];  
1965                              const double A_01 = A_p[INDEX2(0,1,2)];                              const double A_01 = A_p[INDEX2(0,1,2)];
1966                                const double A_10 = A_p[INDEX2(1,0,2)];
1967                              const double A_11 = A_p[INDEX2(1,1,2)];                              const double A_11 = A_p[INDEX2(1,1,2)];
1968                              const double tmp0_0 = A_01 + A_10;                              const double tmp0 = w9*(-A_01 - A_10);
1969                              const double tmp0_1 = A_00*w18;                              const double tmp1 = w9*(-A_01 + A_10);
1970                              const double tmp1_1 = A_01*w19;                              const double tmp2 = w9*(A_01 + A_10);
1971                              const double tmp2_1 = A_10*w20;                              const double tmp3 = w9*(A_01 - A_10);
1972                              const double tmp3_1 = A_11*w21;                              EM_S[INDEX2(0,0,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;
1973                              const double tmp4_1 = A_00*w22;                              EM_S[INDEX2(0,1,4)]+=8*A_00*w6 + 6*A_01*w1 + A_10*w9 + 4*A_11*w3;
1974                              const double tmp5_1 = tmp0_0*w19;                              EM_S[INDEX2(0,2,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;
1975                              const double tmp6_1 = A_11*w23;                              EM_S[INDEX2(0,3,4)]+=4*A_00*w6 + A_11*w10 + tmp2;
1976                              const double tmp7_1 = A_11*w25;                              EM_S[INDEX2(1,0,4)]+=8*A_00*w6 - A_11*w10 + tmp3;
1977                              const double tmp8_1 = A_00*w24;                              EM_S[INDEX2(1,1,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;
1978                              const double tmp9_1 = tmp0_0*w20;                              EM_S[INDEX2(1,2,4)]+=4*A_00*w6 + A_11*w10 + tmp0;
1979                              const double tmp10_1 = A_01*w20;                              EM_S[INDEX2(1,3,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;
1980                              const double tmp11_1 = A_11*w27;                              EM_S[INDEX2(2,0,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;
1981                              const double tmp12_1 = A_00*w26;                              EM_S[INDEX2(2,1,4)]+=4*A_00*w6 + A_11*w10 + tmp0;
1982                              const double tmp13_1 = A_10*w19;                              EM_S[INDEX2(2,2,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;
1983                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                              EM_S[INDEX2(2,3,4)]+=8*A_00*w6 - A_11*w10 + tmp3;
1984                              EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;                              EM_S[INDEX2(3,0,4)]+=4*A_00*w6 + A_11*w10 + tmp2;
1985                              EM_S[INDEX2(2,0,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;                              EM_S[INDEX2(3,1,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;
1986                              EM_S[INDEX2(3,0,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                              EM_S[INDEX2(3,2,4)]+=8*A_00*w6 - A_11*w10 + tmp1;
1987                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                              EM_S[INDEX2(3,3,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;
                             EM_S[INDEX2(1,1,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp4_1 + tmp6_1 + tmp9_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp7_1 + tmp8_1;  
1988                          }                          }
1989                      }                      }
1990                      ///////////////                      ///////////////
# Line 2078  void Rectangle::assemblePDESingle(Paso_S Line 2002  void Rectangle::assemblePDESingle(Paso_S
2002                              const double B_1_2 = B_p[INDEX2(1,2,2)];                              const double B_1_2 = B_p[INDEX2(1,2,2)];
2003                              const double B_0_3 = B_p[INDEX2(0,3,2)];                              const double B_0_3 = B_p[INDEX2(0,3,2)];
2004                              const double B_1_3 = B_p[INDEX2(1,3,2)];                              const double B_1_3 = B_p[INDEX2(1,3,2)];
2005                              const double tmp0_0 = B_1_0 + B_1_1;                              const double tmp0 = w15*(B_1_2 + B_1_3);
2006                              const double tmp1_0 = B_1_2 + B_1_3;                              const double tmp1 = w12*(B_1_0 + B_1_1);
2007                              const double tmp2_0 = B_0_1 + B_0_3;                              const double tmp2 = w15*(B_1_0 + B_1_1);
2008                              const double tmp3_0 = B_0_0 + B_0_2;                              const double tmp3 = w16*(-B_0_1 - B_0_3);
2009                              const double tmp63_1 = B_1_1*w42;                              const double tmp4 = w11*(-B_0_0 - B_0_2);
2010                              const double tmp79_1 = B_1_1*w40;                              const double tmp5 = w12*(B_1_2 + B_1_3);
2011                              const double tmp37_1 = tmp3_0*w35;                              const double tmp6 = w15*(-B_1_0 - B_1_1);
2012                              const double tmp8_1 = tmp0_0*w32;                              const double tmp7 = w12*(-B_1_2 - B_1_3);
2013                              const double tmp71_1 = B_0_1*w34;                              const double tmp8 = w15*(-B_1_2 - B_1_3);
2014                              const double tmp19_1 = B_0_3*w31;                              const double tmp9 = w12*(-B_1_0 - B_1_1);
2015                              const double tmp15_1 = B_0_3*w34;                              const double tmp10 = w11*(-B_0_1 - B_0_3);
2016                              const double tmp9_1 = tmp3_0*w34;                              const double tmp11 = w16*(-B_0_0 - B_0_2);
2017                              const double tmp35_1 = B_1_0*w36;                              const double tmp12 = w16*(B_0_0 + B_0_2);
2018                              const double tmp66_1 = B_0_3*w28;                              const double tmp13 = w11*(B_0_1 + B_0_3);
2019                              const double tmp28_1 = B_1_0*w42;                              const double tmp14 = w11*(B_0_0 + B_0_2);
2020                              const double tmp22_1 = B_1_0*w40;                              const double tmp15 = w16*(B_0_1 + B_0_3);
2021                              const double tmp16_1 = B_1_2*w29;                              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;
2022                              const double tmp6_1 = tmp2_0*w35;                              EM_S[INDEX2(0,1,4)]+=B_0_0*w11 + B_0_1*w13 + B_0_2*w14 + B_0_3*w16 + tmp0 + tmp1;
2023                              const double tmp55_1 = B_1_3*w40;                              EM_S[INDEX2(0,2,4)]+=B_1_0*w12 + B_1_1*w18 + B_1_2*w17 + B_1_3*w15 + tmp14 + tmp15;
2024                              const double tmp50_1 = B_1_3*w42;                              EM_S[INDEX2(0,3,4)]+=tmp12 + tmp13 + tmp2 + tmp5;
2025                              const double tmp7_1 = tmp1_0*w29;                              EM_S[INDEX2(1,0,4)]+=-B_0_0*w13 - B_0_1*w11 - B_0_2*w16 - B_0_3*w14 + tmp0 + tmp1;
2026                              const double tmp1_1 = tmp1_0*w32;                              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;
2027                              const double tmp57_1 = B_0_3*w30;                              EM_S[INDEX2(1,2,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2028                              const double tmp18_1 = B_1_1*w32;                              EM_S[INDEX2(1,3,4)]+=B_1_0*w18 + B_1_1*w12 + B_1_2*w15 + B_1_3*w17 + tmp10 + tmp11;
2029                              const double tmp53_1 = B_1_0*w41;                              EM_S[INDEX2(2,0,4)]+=-B_1_0*w17 - B_1_1*w15 - B_1_2*w12 - B_1_3*w18 + tmp14 + tmp15;
2030                              const double tmp61_1 = B_1_3*w36;                              EM_S[INDEX2(2,1,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2031                              const double tmp27_1 = B_0_3*w38;                              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;
2032                              const double tmp64_1 = B_0_2*w30;                              EM_S[INDEX2(2,3,4)]+=B_0_0*w14 + B_0_1*w16 + B_0_2*w11 + B_0_3*w13 + tmp6 + tmp7;
2033                              const double tmp76_1 = B_0_1*w38;                              EM_S[INDEX2(3,0,4)]+=tmp3 + tmp4 + tmp8 + tmp9;
2034                              const double tmp39_1 = tmp2_0*w34;                              EM_S[INDEX2(3,1,4)]+=-B_1_0*w15 - B_1_1*w17 - B_1_2*w18 - B_1_3*w12 + tmp10 + tmp11;
2035                              const double tmp62_1 = B_0_1*w31;                              EM_S[INDEX2(3,2,4)]+=-B_0_0*w16 - B_0_1*w14 - B_0_2*w13 - B_0_3*w11 + tmp6 + tmp7;
2036                              const double tmp56_1 = B_0_0*w31;                              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;
                             const double tmp49_1 = B_1_1*w36;  
                             const double tmp2_1 = B_0_2*w31;  
                             const double tmp23_1 = B_0_2*w33;  
                             const double tmp38_1 = B_1_1*w43;  
                             const double tmp74_1 = B_1_2*w41;  
                             const double tmp43_1 = B_1_1*w41;  
                             const double tmp58_1 = B_0_2*w28;  
                             const double tmp67_1 = B_0_0*w33;  
                             const double tmp33_1 = tmp0_0*w39;  
                             const double tmp4_1 = B_0_0*w28;  
                             const double tmp20_1 = B_0_0*w30;  
                             const double tmp13_1 = B_0_2*w38;  
                             const double tmp65_1 = B_1_2*w43;  
                             const double tmp0_1 = tmp0_0*w29;  
                             const double tmp41_1 = tmp3_0*w33;  
                             const double tmp73_1 = B_0_2*w37;  
                             const double tmp69_1 = B_0_0*w38;  
                             const double tmp48_1 = B_1_2*w39;  
                             const double tmp59_1 = B_0_1*w33;  
                             const double tmp17_1 = B_1_3*w41;  
                             const double tmp5_1 = B_0_3*w33;  
                             const double tmp3_1 = B_0_1*w30;  
                             const double tmp21_1 = B_0_1*w28;  
                             const double tmp42_1 = B_1_0*w29;  
                             const double tmp54_1 = B_1_2*w32;  
                             const double tmp60_1 = B_1_0*w39;  
                             const double tmp32_1 = tmp1_0*w36;  
                             const double tmp10_1 = B_0_1*w37;  
                             const double tmp14_1 = B_0_0*w35;  
                             const double tmp29_1 = B_0_1*w35;  
                             const double tmp26_1 = B_1_2*w36;  
                             const double tmp30_1 = B_1_3*w43;  
                             const double tmp70_1 = B_0_2*w35;  
                             const double tmp34_1 = B_1_3*w39;  
                             const double tmp51_1 = B_1_0*w43;  
                             const double tmp31_1 = B_0_2*w34;  
                             const double tmp45_1 = tmp3_0*w28;  
                             const double tmp11_1 = tmp1_0*w39;  
                             const double tmp52_1 = B_1_1*w29;  
                             const double tmp44_1 = B_1_3*w32;  
                             const double tmp25_1 = B_1_1*w39;  
                             const double tmp47_1 = tmp2_0*w33;  
                             const double tmp72_1 = B_1_3*w29;  
                             const double tmp40_1 = tmp2_0*w28;  
                             const double tmp46_1 = B_1_2*w40;  
                             const double tmp36_1 = B_1_2*w42;  
                             const double tmp24_1 = B_0_0*w37;  
                             const double tmp77_1 = B_0_3*w35;  
                             const double tmp68_1 = B_0_3*w37;  
                             const double tmp78_1 = B_0_0*w34;  
                             const double tmp12_1 = tmp0_0*w36;  
                             const double tmp75_1 = B_1_0*w32;  
                             EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
2037                          } else { // constant data                          } else { // constant data
2038                              const double B_0 = B_p[0];                              const double B_0 = B_p[0];
2039                              const double B_1 = B_p[1];                              const double B_1 = B_p[1];
2040                              const double tmp0_1 = B_0*w44;                              EM_S[INDEX2(0,0,4)]+=B_0*w19 + 2*B_1*w20;
2041                              const double tmp1_1 = B_1*w45;                              EM_S[INDEX2(0,1,4)]+=B_0*w19 + B_1*w20;
2042                              const double tmp2_1 = B_0*w46;                              EM_S[INDEX2(0,2,4)]+=B_0*w19/2 + 2*B_1*w20;
2043                              const double tmp3_1 = B_0*w47;                              EM_S[INDEX2(0,3,4)]+=B_0*w19/2 + B_1*w20;
2044                              const double tmp4_1 = B_1*w48;                              EM_S[INDEX2(1,0,4)]+=-B_0*w19 + B_1*w20;
2045                              const double tmp5_1 = B_1*w49;                              EM_S[INDEX2(1,1,4)]+=-B_0*w19 + 2*B_1*w20;
2046                              const double tmp6_1 = B_1*w50;                              EM_S[INDEX2(1,2,4)]+=-B_0*w19/2 + B_1*w20;
2047                              const double tmp7_1 = B_0*w51;                              EM_S[INDEX2(1,3,4)]+=-B_0*w19/2 + 2*B_1*w20;
2048                              EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(2,0,4)]+=B_0*w19/2 - 2*B_1*w20;
2049                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(2,1,4)]+=B_0*w19/2 - B_1*w20;
2050                              EM_S[INDEX2(2,0,4)]+=tmp6_1 + tmp7_1;                              EM_S[INDEX2(2,2,4)]+=B_0*w19 - 2*B_1*w20;
2051                              EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp4_1;                              EM_S[INDEX2(2,3,4)]+=B_0*w19 - B_1*w20;
2052                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(3,0,4)]+=-B_0*w19/2 - B_1*w20;
2053                              EM_S[INDEX2(1,1,4)]+=tmp3_1 + tmp5_1;                              EM_S[INDEX2(3,1,4)]+=-B_0*w19/2 - 2*B_1*w20;
2054                              EM_S[INDEX2(2,1,4)]+=tmp4_1 + tmp7_1;                              EM_S[INDEX2(3,2,4)]+=-B_0*w19 - B_1*w20;
2055                              EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp6_1;                              EM_S[INDEX2(3,3,4)]+=-B_0*w19 - 2*B_1*w20;
                             EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp7_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp7_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp5_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp4_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp3_1 + tmp6_1;  
2056                          }                          }
2057                      }                      }
2058                      ///////////////                      ///////////////
# Line 2222  void Rectangle::assemblePDESingle(Paso_S Line 2070  void Rectangle::assemblePDESingle(Paso_S
2070                              const double C_1_2 = C_p[INDEX2(1,2,2)];                              const double C_1_2 = C_p[INDEX2(1,2,2)];
2071                              const double C_0_3 = C_p[INDEX2(0,3,2)];                              const double C_0_3 = C_p[INDEX2(0,3,2)];
2072                              const double C_1_3 = C_p[INDEX2(1,3,2)];                              const double C_1_3 = C_p[INDEX2(1,3,2)];
2073                              const double tmp0_0 = C_1_0 + C_1_1;                              const double tmp0 = w15*(C_1_2 + C_1_3);
2074                              const double tmp1_0 = C_1_2 + C_1_3;                              const double tmp1 = w12*(C_1_0 + C_1_1);
2075                              const double tmp2_0 = C_0_1 + C_0_3;                              const double tmp2 = w15*(-C_1_2 - C_1_3);
2076                              const double tmp3_0 = C_0_0 + C_0_2;                              const double tmp3 = w16*(C_0_0 + C_0_2);
2077                              const double tmp64_1 = C_0_2*w30;                              const double tmp4 = w11*(C_0_1 + C_0_3);
2078                              const double tmp14_1 = C_0_2*w28;                              const double tmp5 = w12*(-C_1_0 - C_1_1);
2079                              const double tmp19_1 = C_0_3*w31;                              const double tmp6 = w15*(-C_1_0 - C_1_1);
2080                              const double tmp22_1 = C_1_0*w40;                              const double tmp7 = w12*(-C_1_2 - C_1_3);
2081                              const double tmp37_1 = tmp3_0*w35;                              const double tmp8 = w15*(C_1_0 + C_1_1);
2082                              const double tmp29_1 = C_0_1*w35;                              const double tmp9 = w12*(C_1_2 + C_1_3);
2083                              const double tmp73_1 = C_0_2*w37;                              const double tmp10 = w11*(-C_0_1 - C_0_3);
2084                              const double tmp74_1 = C_1_2*w41;                              const double tmp11 = w16*(-C_0_0 - C_0_2);
2085                              const double tmp52_1 = C_1_3*w39;                              const double tmp12 = w16*(-C_0_1 - C_0_3);
2086                              const double tmp25_1 = C_1_1*w39;                              const double tmp13 = w11*(-C_0_0 - C_0_2);
2087                              const double tmp62_1 = C_0_1*w31;                              const double tmp14 = w11*(C_0_0 + C_0_2);
2088                              const double tmp79_1 = C_1_1*w40;                              const double tmp15 = w16*(C_0_1 + C_0_3);
2089                              const double tmp43_1 = C_1_1*w36;                              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;
2090                              const double tmp27_1 = C_0_3*w38;                              EM_S[INDEX2(0,1,4)]+=-C_0_0*w13 - C_0_1*w11 - C_0_2*w16 - C_0_3*w14 + tmp0 + tmp1;
2091                              const double tmp28_1 = C_1_0*w42;                              EM_S[INDEX2(0,2,4)]+=-C_1_0*w17 - C_1_1*w15 - C_1_2*w12 - C_1_3*w18 + tmp14 + tmp15;
2092                              const double tmp63_1 = C_1_1*w42;                              EM_S[INDEX2(0,3,4)]+=tmp12 + tmp13 + tmp2 + tmp5;
2093                              const double tmp59_1 = C_0_3*w34;                              EM_S[INDEX2(1,0,4)]+=C_0_0*w11 + C_0_1*w13 + C_0_2*w14 + C_0_3*w16 + tmp0 + tmp1;
2094                              const double tmp72_1 = C_1_3*w29;                              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;
2095                              const double tmp40_1 = tmp2_0*w35;                              EM_S[INDEX2(1,2,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2096                              const double tmp13_1 = C_0_3*w30;                              EM_S[INDEX2(1,3,4)]+=-C_1_0*w15 - C_1_1*w17 - C_1_2*w18 - C_1_3*w12 + tmp10 + tmp11;
2097                              const double tmp51_1 = C_1_2*w40;                              EM_S[INDEX2(2,0,4)]+=C_1_0*w12 + C_1_1*w18 + C_1_2*w17 + C_1_3*w15 + tmp14 + tmp15;
2098                              const double tmp54_1 = C_1_2*w42;                              EM_S[INDEX2(2,1,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2099                              const double tmp12_1 = C_0_0*w31;                              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;
2100                              const double tmp2_1 = tmp1_0*w32;                              EM_S[INDEX2(2,3,4)]+=-C_0_0*w16 - C_0_1*w14 - C_0_2*w13 - C_0_3*w11 + tmp6 + tmp7;
2101                              const double tmp68_1 = C_0_2*w31;                              EM_S[INDEX2(3,0,4)]+=tmp3 + tmp4 + tmp8 + tmp9;
2102                              const double tmp75_1 = C_1_0*w32;                              EM_S[INDEX2(3,1,4)]+=C_1_0*w18 + C_1_1*w12 + C_1_2*w15 + C_1_3*w17 + tmp10 + tmp11;
2103                              const double tmp49_1 = C_1_1*w41;                              EM_S[INDEX2(3,2,4)]+=C_0_0*w14 + C_0_1*w16 + C_0_2*w11 + C_0_3*w13 + tmp6 + tmp7;
2104                              const double tmp4_1 = C_0_2*w35;                              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;
                             const double tmp66_1 = C_0_3*w28;  
                             const double tmp56_1 = C_0_1*w37;  
                             const double tmp5_1 = C_0_1*w34;  
                             const double tmp38_1 = tmp2_0*w34;  
                             const double tmp76_1 = C_0_1*w38;  
                             const double tmp21_1 = C_0_1*w28;  
                             const double tmp69_1 = C_0_1*w30;  
                             const double tmp53_1 = C_1_0*w36;  
                             const double tmp42_1 = C_1_2*w39;  
                             const double tmp32_1 = tmp1_0*w29;  
                             const double tmp45_1 = C_1_0*w43;  
                             const double tmp33_1 = tmp0_0*w32;  
                             const double tmp35_1 = C_1_0*w41;  
                             const double tmp26_1 = C_1_2*w36;  
                             const double tmp67_1 = C_0_0*w33;  
                             const double tmp31_1 = C_0_2*w34;  
                             const double tmp20_1 = C_0_0*w30;  
                             const double tmp70_1 = C_0_0*w28;  
                             const double tmp8_1 = tmp0_0*w39;  
                             const double tmp30_1 = C_1_3*w43;  
                             const double tmp0_1 = tmp0_0*w29;  
                             const double tmp17_1 = C_1_3*w41;  
                             const double tmp58_1 = C_0_0*w35;  
                             const double tmp9_1 = tmp3_0*w33;  
                             const double tmp61_1 = C_1_3*w36;  
                             const double tmp41_1 = tmp3_0*w34;  
                             const double tmp50_1 = C_1_3*w32;  
                             const double tmp18_1 = C_1_1*w32;  
                             const double tmp6_1 = tmp1_0*w36;  
                             const double tmp3_1 = C_0_0*w38;  
                             const double tmp34_1 = C_1_1*w29;  
                             const double tmp77_1 = C_0_3*w35;  
                             const double tmp65_1 = C_1_2*w43;  
                             const double tmp71_1 = C_0_3*w33;  
                             const double tmp55_1 = C_1_1*w43;  
                             const double tmp46_1 = tmp3_0*w28;  
                             const double tmp24_1 = C_0_0*w37;  
                             const double tmp10_1 = tmp1_0*w39;  
                             const double tmp48_1 = C_1_0*w29;  
                             const double tmp15_1 = C_0_1*w33;  
                             const double tmp36_1 = C_1_2*w32;  
                             const double tmp60_1 = C_1_0*w39;  
                             const double tmp47_1 = tmp2_0*w33;  
                             const double tmp16_1 = C_1_2*w29;  
                             const double tmp1_1 = C_0_3*w37;  
                             const double tmp7_1 = tmp2_0*w28;  
                             const double tmp39_1 = C_1_3*w40;  
                             const double tmp44_1 = C_1_3*w42;  
                             const double tmp57_1 = C_0_2*w38;  
                             const double tmp78_1 = C_0_0*w34;  
                             const double tmp11_1 = tmp0_0*w36;  
                             const double tmp23_1 = C_0_2*w33;  
                             EM_S[INDEX2(0,0,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
2105                          } else { // constant data                          } else { // constant data
2106                              const double C_0 = C_p[0];                              const double C_0 = C_p[0];
2107                              const double C_1 = C_p[1];                              const double C_1 = C_p[1];
2108                              const double tmp0_1 = C_0*w47;                              EM_S[INDEX2(0,0,4)]+=C_0*w19 + 2*C_1*w20;
2109                              const double tmp1_1 = C_1*w45;                              EM_S[INDEX2(0,1,4)]+=-C_0*w19 + C_1*w20;
2110                              const double tmp2_1 = C_1*w48;                              EM_S[INDEX2(0,2,4)]+=C_0*w19/2 - 2*C_1*w20;
2111                              const double tmp3_1 = C_0*w51;                              EM_S[INDEX2(0,3,4)]+=-C_0*w19/2 - C_1*w20;
2112                              const double tmp4_1 = C_0*w44;                              EM_S[INDEX2(1,0,4)]+=C_0*w19 + C_1*w20;
2113                              const double tmp5_1 = C_1*w49;                              EM_S[INDEX2(1,1,4)]+=-C_0*w19 + 2*C_1*w20;
2114                              const double tmp6_1 = C_1*w50;                              EM_S[INDEX2(1,2,4)]+=C_0*w19/2 - C_1*w20;
2115                              const double tmp7_1 = C_0*w46;                              EM_S[INDEX2(1,3,4)]+=-C_0*w19/2 - 2*C_1*w20;
2116                              EM_S[INDEX2(0,0,4)]+=tmp4_1 + tmp5_1;                              EM_S[INDEX2(2,0,4)]+=C_0*w19/2 + 2*C_1*w20;
2117                              EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp4_1;                              EM_S[INDEX2(2,1,4)]+=-C_0*w19/2 + C_1*w20;
2118                              EM_S[INDEX2(2,0,4)]+=tmp3_1 + tmp5_1;                              EM_S[INDEX2(2,2,4)]+=C_0*w19 - 2*C_1*w20;
2119                              EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;                              EM_S[INDEX2(2,3,4)]+=-C_0*w19 - C_1*w20;
2120                              EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                              EM_S[INDEX2(3,0,4)]+=C_0*w19/2 + C_1*w20;
2121                              EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp5_1;                              EM_S[INDEX2(3,1,4)]+=-C_0*w19/2 + 2*C_1*w20;
2122                              EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp7_1;                              EM_S[INDEX2(3,2,4)]+=C_0*w19 - C_1*w20;
2123                              EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp7_1;                              EM_S[INDEX2(3,3,4)]+=-C_0*w19 - 2*C_1*w20;
                             EM_S[INDEX2(0,2,4)]+=tmp3_1 + tmp6_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp4_1 + tmp6_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp4_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp2_1 + tmp7_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp6_1 + tmp7_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp6_1;  
2124                          }                          }
2125                      }                      }
2126                      ///////////////                      ///////////////
# Line 2362  void Rectangle::assemblePDESingle(Paso_S Line 2134  void Rectangle::assemblePDESingle(Paso_S
2134                              const double D_1 = D_p[1];                              const double D_1 = D_p[1];
2135                              const double D_2 = D_p[2];                              const double D_2 = D_p[2];
2136                              const double D_3 = D_p[3];                              const double D_3 = D_p[3];
2137                              const double tmp0_0 = D_0 + D_1;                              const double tmp0 = w21*(D_0 + D_1);
2138                              const double tmp1_0 = D_2 + D_3;                              const double tmp1 = w22*(D_2 + D_3);
2139                              const double tmp2_0 = D_0 + D_1 + D_2 + D_3;                              const double tmp2 = w23*(D_0 + D_1 + D_2 + D_3);
2140                              const double tmp3_0 = D_1 + D_2;                              const double tmp3 = w21*(D_2 + D_3);
2141                              const double tmp4_0 = D_1 + D_3;                              const double tmp4 = w22*(D_0 + D_1);
2142                              const double tmp5_0 = D_0 + D_2;                              const double tmp5 = w23*(D_1 + D_2);
2143                              const double tmp6_0 = D_0 + D_3;                              const double tmp6 = w21*(D_1 + D_3);
2144                              const double tmp0_1 = tmp0_0*w52;                              const double tmp7 = w22*(D_0 + D_2);
2145                              const double tmp1_1 = tmp1_0*w53;                              const double tmp8 = w21*(D_0 + D_2);
2146                              const double tmp2_1 = tmp2_0*w54;                              const double tmp9 = w22*(D_1 + D_3);
2147                              const double tmp3_1 = tmp1_0*w52;                              const double tmp10 = w23*(D_0 + D_3);
2148                              const double tmp4_1 = tmp0_0*w53;                              EM_S[INDEX2(0,0,4)]+=D_0*w24 + D_3*w25 + tmp5;
2149                              const double tmp5_1 = tmp3_0*w54;                              EM_S[INDEX2(0,1,4)]+=tmp0 + tmp1;
2150                              const double tmp6_1 = D_0*w55;                              EM_S[INDEX2(0,2,4)]+=tmp8 + tmp9;
2151                              const double tmp7_1 = D_3*w56;                              EM_S[INDEX2(0,3,4)]+=tmp2;
2152                              const double tmp8_1 = D_3*w55;                              EM_S[INDEX2(1,0,4)]+=tmp0 + tmp1;
2153                              const double tmp9_1 = D_0*w56;                              EM_S[INDEX2(1,1,4)]+=D_1*w24 + D_2*w25 + tmp10;
2154                              const double tmp10_1 = tmp4_0*w52;                              EM_S[INDEX2(1,2,4)]+=tmp2;
2155                              const double tmp11_1 = tmp5_0*w53;                              EM_S[INDEX2(1,3,4)]+=tmp6 + tmp7;
2156                              const double tmp12_1 = tmp5_0*w52;                              EM_S[INDEX2(2,0,4)]+=tmp8 + tmp9;
2157                              const double tmp13_1 = tmp4_0*w53;                              EM_S[INDEX2(2,1,4)]+=tmp2;
2158                              const double tmp14_1 = tmp6_0*w54;                              EM_S[INDEX2(2,2,4)]+=D_1*w25 + D_2*w24 + tmp10;
2159                              const double tmp15_1 = D_2*w55;                              EM_S[INDEX2(2,3,4)]+=tmp3 + tmp4;
2160                              const double tmp16_1 = D_1*w56;                              EM_S[INDEX2(3,0,4)]+=tmp2;
2161                              const double tmp17_1 = D_1*w55;                              EM_S[INDEX2(3,1,4)]+=tmp6 + tmp7;
2162                              const double tmp18_1 = D_2*w56;                              EM_S[INDEX2(3,2,4)]+=tmp3 + tmp4;
2163                              EM_S[INDEX2(0,0,4)]+=tmp5_1 + tmp6_1 + tmp7_1;                              EM_S[INDEX2(3,3,4)]+=D_0*w25 + D_3*w24 + tmp5;
                             EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(2,0,4)]+=tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(3,0,4)]+=tmp2_1;  
                             EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;  
                             EM_S[INDEX2(1,1,4)]+=tmp14_1 + tmp17_1 + tmp18_1;  
                             EM_S[INDEX2(2,1,4)]+=tmp2_1;  
                             EM_S[INDEX2(3,1,4)]+=tmp10_1 + tmp11_1;  
                             EM_S[INDEX2(0,2,4)]+=tmp12_1 + tmp13_1;  
                             EM_S[INDEX2(1,2,4)]+=tmp2_1;  
                             EM_S[INDEX2(2,2,4)]+=tmp14_1 + tmp15_1 + tmp16_1;  
                             EM_S[INDEX2(3,2,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(0,3,4)]+=tmp2_1;  
                             EM_S[INDEX2(1,3,4)]+=tmp10_1 + tmp11_1;  
                             EM_S[INDEX2(2,3,4)]+=tmp3_1 + tmp4_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp5_1 + tmp8_1 + tmp9_1;  
2164                          } else { // constant data                          } else { // constant data
2165                              const double tmp0_1 = D_p[0]*w57;                              const double D_0 = D_p[0];
2166                              const double tmp1_1 = D_p[0]*w58;                              EM_S[INDEX2(0,0,4)]+=16*D_0*w23;
2167                              const double tmp2_1 = D_p[0]*w59;                              EM_S[INDEX2(0,1,4)]+=8*D_0*w23;
2168                              EM_S[INDEX2(0,0,4)]+=tmp2_1;                              EM_S[INDEX2(0,2,4)]+=8*D_0*w23;
2169                              EM_S[INDEX2(1,0,4)]+=tmp0_1;                              EM_S[INDEX2(0,3,4)]+=4*D_0*w23;
2170                              EM_S[INDEX2(2,0,4)]+=tmp0_1;                              EM_S[INDEX2(1,0,4)]+=8*D_0*w23;
2171                              EM_S[INDEX2(3,0,4)]+=tmp1_1;                              EM_S[INDEX2(1,1,4)]+=16*D_0*w23;
2172                              EM_S[INDEX2(0,1,4)]+=tmp0_1;                              EM_S[INDEX2(1,2,4)]+=4*D_0*w23;
2173                              EM_S[INDEX2(1,1,4)]+=tmp2_1;                              EM_S[INDEX2(1,3,4)]+=8*D_0*w23;
2174                              EM_S[INDEX2(2,1,4)]+=tmp1_1;                              EM_S[INDEX2(2,0,4)]+=8*D_0*w23;
2175                              EM_S[INDEX2(3,1,4)]+=tmp0_1;                              EM_S[INDEX2(2,1,4)]+=4*D_0*w23;
2176                              EM_S[INDEX2(0,2,4)]+=tmp0_1;                              EM_S[INDEX2(2,2,4)]+=16*D_0*w23;
2177                              EM_S[INDEX2(1,2,4)]+=tmp1_1;                              EM_S[INDEX2(2,3,4)]+=8*D_0*w23;
2178                              EM_S[INDEX2(2,2,4)]+=tmp2_1;                              EM_S[INDEX2(3,0,4)]+=4*D_0*w23;
2179                              EM_S[INDEX2(3,2,4)]+=tmp0_1;                              EM_S[INDEX2(3,1,4)]+=8*D_0*w23;
2180                              EM_S[INDEX2(0,3,4)]+=tmp1_1;                              EM_S[INDEX2(3,2,4)]+=8*D_0*w23;
2181                              EM_S[INDEX2(1,3,4)]+=tmp0_1;                              EM_S[INDEX2(3,3,4)]+=16*D_0*w23;
                             EM_S[INDEX2(2,3,4)]+=tmp0_1;  
                             EM_S[INDEX2(3,3,4)]+=tmp2_1;  
2182                          }                          }
2183                      }                      }
2184                      ///////////////                      ///////////////
# Line 2441  void Rectangle::assemblePDESingle(Paso_S Line 2196  void Rectangle::assemblePDESingle(Paso_S
2196                              const double X_1_2 = X_p[INDEX2(1,2,2)];                              const double X_1_2 = X_p[INDEX2(1,2,2)];
2197                              const double X_0_3 = X_p[INDEX2(0,3,2)];                              const double X_0_3 = X_p[INDEX2(0,3,2)];
2198                              const double X_1_3 = X_p[INDEX2(1,3,2)];                              const double X_1_3 = X_p[INDEX2(1,3,2)];
2199                              const double tmp0_0 = X_0_2 + X_0_3;                              const double tmp0 = 6*w15*(X_1_1 + X_1_3);
2200                              const double tmp1_0 = X_1_1 + X_1_3;                              const double tmp1 = 6*w16*(X_0_2 + X_0_3);
2201                              const double tmp2_0 = X_1_0 + X_1_2;                              const double tmp2 = 6*w11*(X_0_0 + X_0_1);
2202                              const double tmp3_0 = X_0_0 + X_0_1;                              const double tmp3 = 6*w12*(X_1_0 + X_1_2);
2203                              const double tmp0_1 = tmp0_0*w63;                              const double tmp4 = 6*w15*(X_1_0 + X_1_2);
2204                              const double tmp1_1 = tmp1_0*w62;                              const double tmp5 = w27*(X_0_2 + X_0_3);
2205                              const double tmp2_1 = tmp2_0*w61;                              const double tmp6 = w26*(X_0_0 + X_0_1);
2206                              const double tmp3_1 = tmp3_0*w60;                              const double tmp7 = 6*w12*(X_1_1 + X_1_3);
2207                              const double tmp4_1 = tmp0_0*w65;                              const double tmp8 = w29*(X_1_1 + X_1_3);
2208                              const double tmp5_1 = tmp3_0*w64;                              const double tmp9 = w27*(-X_0_0 - X_0_1);
2209                              const double tmp6_1 = tmp2_0*w62;                              const double tmp10 = w28*(X_1_0 + X_1_2);
2210                              const double tmp7_1 = tmp1_0*w61;                              const double tmp11 = w26*(-X_0_2 - X_0_3);
2211                              const double tmp8_1 = tmp2_0*w66;                              const double tmp12 = w29*(X_1_0 + X_1_2);
2212                              const double tmp9_1 = tmp3_0*w63;                              const double tmp13 = w27*(X_0_0 + X_0_1);
2213                              const double tmp10_1 = tmp0_0*w60;                              const double tmp14 = w28*(X_1_1 + X_1_3);
2214                              const double tmp11_1 = tmp1_0*w67;                              const double tmp15 = w26*(X_0_2 + X_0_3);
2215                              const double tmp12_1 = tmp1_0*w66;                              EM_F[0]+=tmp0 + tmp1 + tmp2 + tmp3;
2216                              const double tmp13_1 = tmp3_0*w65;                              EM_F[1]+=tmp4 + tmp5 + tmp6 + tmp7;
2217                              const double tmp14_1 = tmp0_0*w64;                              EM_F[2]+=tmp10 + tmp11 + tmp8 + tmp9;
2218                              const double tmp15_1 = tmp2_0*w67;                              EM_F[3]+=tmp12 + tmp13 + tmp14 + tmp15;
                             EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                             EM_F[1]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;  
                             EM_F[2]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;  
                             EM_F[3]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
2219                          } else { // constant data                          } else { // constant data
2220                              const double tmp0_1 = X_p[1]*w69;                              const double X_0 = X_p[0];
2221                              const double tmp1_1 = X_p[0]*w68;                              const double X_1 = X_p[1];
2222                              const double tmp2_1 = X_p[0]*w70;                              EM_F[0]+=3*X_0*w19 + 6*X_1*w20;
2223                              const double tmp3_1 = X_p[1]*w71;                              EM_F[1]+=-3*X_0*w19 + 6*X_1*w20;
2224                              EM_F[0]+=tmp0_1 + tmp1_1;                              EM_F[2]+=3*X_0*w19 - 6*X_1*w20;
2225                              EM_F[1]+=tmp0_1 + tmp2_1;                              EM_F[3]+=-3*X_0*w19 - 6*X_1*w20;
                             EM_F[2]+=tmp1_1 + tmp3_1;  
                             EM_F[3]+=tmp2_1 + tmp3_1;  
2226                          }                          }
2227                      }                      }
2228                      ///////////////                      ///////////////
# Line 2487  void Rectangle::assemblePDESingle(Paso_S Line 2236  void Rectangle::assemblePDESingle(Paso_S
2236                              const double Y_1 = Y_p[1];                              const double Y_1 = Y_p[1];
2237                              const double Y_2 = Y_p[2];                              const double Y_2 = Y_p[2];
2238                              const double Y_3 = Y_p[3];                              const double Y_3 = Y_p[3];
2239                              const double tmp0_0 = Y_1 + Y_2;                              const double tmp0 = 6*w23*(Y_1 + Y_2);
2240                              const double tmp1_0 = Y_0 + Y_3;                              const double tmp1 = 6*w23*(Y_0 + Y_3);
2241                              const double tmp0_1 = Y_0*w72;                              EM_F[0]+=6*Y_0*w21 + 6*Y_3*w22 + tmp0;
2242                              const double tmp1_1 = tmp0_0*w73;                              EM_F[1]+=6*Y_1*w21 + 6*Y_2*w22 + tmp1;
2243                              const double tmp2_1 = Y_3*w74;                              EM_F[2]+=6*Y_1*w22 + 6*Y_2*w21 + tmp1;
2244                              const double tmp3_1 = Y_1*w72;                              EM_F[3]+=6*Y_0*w22 + 6*Y_3*w21 + tmp0;
                             const double tmp4_1 = tmp1_0*w73;  
                             const double tmp5_1 = Y_2*w74;  
                             const double tmp6_1 = Y_2*w72;  
                             const double tmp7_1 = Y_1*w74;  
                             const double tmp8_1 = Y_3*w72;  
                             const double tmp9_1 = Y_0*w74;  
                             EM_F[0]+=tmp0_1 + tmp1_1 + tmp2_1;  
                             EM_F[1]+=tmp3_1 + tmp4_1 + tmp5_1;  
                             EM_F[2]+=tmp4_1 + tmp6_1 + tmp7_1;  
                             EM_F[3]+=tmp1_1 + tmp8_1 + tmp9_1;  
2245                          } else { // constant data                          } else { // constant data
2246                              const double tmp0_1 = Y_p[0]*w75;                              const double Y_0 = Y_p[0];
2247                              EM_F[0]+=tmp0_1;                              EM_F[0]+=36*Y_0*w23;
2248                              EM_F[1]+=tmp0_1;                              EM_F[1]+=36*Y_0*w23;
2249                              EM_F[2]+=tmp0_1;                              EM_F[2]+=36*Y_0*w23;
2250                              EM_F[3]+=tmp0_1;                              EM_F[3]+=36*Y_0*w23;
2251                          }                          }
2252                      }                      }
2253                        /* GENERATOR SNIP_PDE_SINGLE BOTTOM */
2254    
2255                      // 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)
2256                      const index_t firstNode=m_NN[0]*k1+k0;                      const index_t firstNode=m_NN[0]*k1+k0;
# Line 2527  void Rectangle::assemblePDESingleReduced Line 2267  void Rectangle::assemblePDESingleReduced
2267          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
2268          const escript::Data& X, const escript::Data& Y) const          const escript::Data& X, const escript::Data& Y) const
2269  {  {
2270      const double w0 = -.25*m_dx[1]/m_dx[0];      const double w0 = .25;
2271      const double w1 = .25;      const double w1 = .125*m_dx[0];
2272      const double w2 = -.25;      const double w2 = .125*m_dx[1];
2273      const double w3 = .25*m_dx[0]/m_dx[1];      const double w3 = .0625*m_dx[0]*m_dx[1];
2274      const double w4 = -.25*m_dx[0]/m_dx[1];      const double w4 = .25*m_dx[0]/m_dx[1];
2275      const double w5 = .25*m_dx[1]/m_dx[0];      const double w5 = .25*m_dx[1]/m_dx[0];
     const double w6 = -.125*m_dx[1];  
     const double w7 = -.125*m_dx[0];  
     const double w8 = .125*m_dx[1];  
     const double w9 = .125*m_dx[0];  
     const double w10 = .0625*m_dx[0]*m_dx[1];  
     const double w11 = -.5*m_dx[1];  
     const double w12 = -.5*m_dx[0];  
     const double w13 = .5*m_dx[1];  
     const double w14 = .5*m_dx[0];  
     const double w15 = .25*m_dx[0]*m_dx[1];  
2276    
2277      rhs.requireWrite();      rhs.requireWrite();
2278  #pragma omp parallel  #pragma omp parallel
# Line 2566  void Rectangle::assemblePDESingleReduced Line 2296  void Rectangle::assemblePDESingleReduced
2296                          const double A_10 = A_p[INDEX2(1,0,2)];                          const double A_10 = A_p[INDEX2(1,0,2)];
2297                          const double A_01 = A_p[INDEX2(0,1,2)];                          const double A_01 = A_p[INDEX2(0,1,2)];
2298                          const double A_11 = A_p[INDEX2(1,1,2)];                          const double A_11 = A_p[INDEX2(1,1,2)];
2299                          const double tmp0_0 = A_01 + A_10;                          const double tmp0 = (A_01 + A_10)*w0;
2300                          const double tmp0_1 = A_11*w3;                          const double tmp1 = A_00*w5;
2301                          const double tmp1_1 = A_00*w0;                          const double tmp2 = A_01*w0;
2302                          const double tmp2_1 = A_01*w1;                          const double tmp3 = A_10*w0;
2303                          const double tmp3_1 = A_10*w2;                          const double tmp4 = A_11*w4;
2304                          const double tmp4_1 = tmp0_0*w1;                          EM_S[INDEX2(0,0,4)]+=tmp4 + tmp0 + tmp1;
2305                          const double tmp5_1 = A_11*w4;                          EM_S[INDEX2(1,0,4)]+=tmp4 - tmp1 + tmp3 - tmp2;
2306                          const double tmp6_1 = A_00*w5;                          EM_S[INDEX2(2,0,4)]+=tmp2 - tmp3 - tmp4 + tmp1;
2307                          const double tmp7_1 = tmp0_0*w2;                          EM_S[INDEX2(3,0,4)]+=-tmp1 - tmp4 - tmp0;
2308                          const double tmp8_1 = A_10*w1;                          EM_S[INDEX2(0,1,4)]+=tmp4 - tmp1 + tmp2 - tmp3;
2309                          const double tmp9_1 = A_01*w2;                          EM_S[INDEX2(1,1,4)]+=tmp4 + tmp1 - tmp0;
2310                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;                          EM_S[INDEX2(2,1,4)]+=-tmp1 + tmp0 - tmp4;
2311                          EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;                          EM_S[INDEX2(3,1,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;
2312                          EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;                          EM_S[INDEX2(0,2,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;
2313                          EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;                          EM_S[INDEX2(1,2,4)]+=-tmp1 + tmp0 - tmp4;
2314                          EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                          EM_S[INDEX2(2,2,4)]+=tmp4 + tmp1 - tmp0;
2315                          EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;                          EM_S[INDEX2(3,2,4)]+=tmp4 - tmp1 + tmp2 - tmp3;
2316                          EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;                          EM_S[INDEX2(0,3,4)]+=-tmp1 - tmp4 - tmp0;
2317                          EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;                          EM_S[INDEX2(1,3,4)]+=tmp2 - tmp3 - tmp4 + tmp1;
2318                          EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;                          EM_S[INDEX2(2,3,4)]+=tmp4 - tmp1 + tmp3 - tmp2;
2319                          EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;                          EM_S[INDEX2(3,3,4)]+=tmp4 + tmp0 + tmp1;
                         EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;  
                         EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                         EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;  
                         EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;  
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;  
2320                      }                      }
2321                      ///////////////                      ///////////////
2322                      // process B //                      // process B //
# Line 2600  void Rectangle::assemblePDESingleReduced Line 2324  void Rectangle::assemblePDESingleReduced
2324                      if (!B.isEmpty()) {                      if (!B.isEmpty()) {
2325                          add_EM_S=true;                          add_EM_S=true;
2326                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2327                          const double tmp2_1 = B_p[0]*w8;                          const double tmp0 = B_p[0]*w2;
2328                          const double tmp0_1 = B_p[0]*w6;                          const double tmp1 = B_p[1]*w1;
2329                          const double tmp3_1 = B_p[1]*w9;                          EM_S[INDEX2(0,0,4)]+=-tmp0 - tmp1;
2330                          const double tmp1_1 = B_p[1]*w7;                          EM_S[INDEX2(1,0,4)]+= tmp0 - tmp1;
2331                          EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(2,0,4)]+= tmp1 - tmp0;
2332                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(3,0,4)]+= tmp0 + tmp1;
2333                          EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;                          EM_S[INDEX2(0,1,4)]+=-tmp0 - tmp1;
2334                          EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(1,1,4)]+= tmp0 - tmp1;
2335                          EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(2,1,4)]+= tmp1 - tmp0;
2336                          EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(3,1,4)]+= tmp0 + tmp1;
2337                          EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;                          EM_S[INDEX2(0,2,4)]+=-tmp0 - tmp1;
2338                          EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(1,2,4)]+= tmp0 - tmp1;
2339                          EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(2,2,4)]+= tmp1 - tmp0;
2340                          EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(3,2,4)]+= tmp0 + tmp1;
2341                          EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;                          EM_S[INDEX2(0,3,4)]+=-tmp0 - tmp1;
2342                          EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(1,3,4)]+= tmp0 - tmp1;
2343                          EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(2,3,4)]+= tmp1 - tmp0;
2344                          EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;                          EM_S[INDEX2(3,3,4)]+= tmp0 + tmp1;
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;  
2345                      }                      }
2346                      ///////////////                      ///////////////
2347                      // process C //                      // process C //
# Line 2627  void Rectangle::assemblePDESingleReduced Line 2349  void Rectangle::assemblePDESingleReduced
2349                      if (!C.isEmpty()) {                      if (!C.isEmpty()) {
2350                          add_EM_S=true;                          add_EM_S=true;
2351                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2352                          const double tmp1_1 = C_p[1]*w7;                          const double tmp0 = C_p[0]*w2;
2353                          const double tmp0_1 = C_p[0]*w8;                          const double tmp1 = C_p[1]*w1;
2354                          const double tmp3_1 = C_p[0]*w6;                          EM_S[INDEX2(0,0,4)]+=-tmp1 - tmp0;
2355                          const double tmp2_1 = C_p[1]*w9;                          EM_S[INDEX2(1,0,4)]+=-tmp1 - tmp0;
2356                          EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(2,0,4)]+=-tmp1 - tmp0;
2357                          EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(3,0,4)]+=-tmp1 - tmp0;
2358                          EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(0,1,4)]+= tmp0 - tmp1;
2359                          EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;                          EM_S[INDEX2(1,1,4)]+= tmp0 - tmp1;
2360                          EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(2,1,4)]+= tmp0 - tmp1;
2361                          EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(3,1,4)]+= tmp0 - tmp1;
2362                          EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(0,2,4)]+= tmp1 - tmp0;
2363                          EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;                          EM_S[INDEX2(1,2,4)]+= tmp1 - tmp0;
2364                          EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(2,2,4)]+= tmp1 - tmp0;
2365                          EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(3,2,4)]+= tmp1 - tmp0;
2366                          EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(0,3,4)]+= tmp0 + tmp1;
2367                          EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;                          EM_S[INDEX2(1,3,4)]+= tmp0 + tmp1;
2368                          EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;                          EM_S[INDEX2(2,3,4)]+= tmp0 + tmp1;
2369                          EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;                          EM_S[INDEX2(3,3,4)]+= tmp0 + tmp1;
                         EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;  
                         EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;  
2370                      }                      }
2371                      ///////////////                      ///////////////
2372                      // process D //                      // process D //
# Line 2654  void Rectangle::assemblePDESingleReduced Line 2374  void Rectangle::assemblePDESingleReduced
2374                      if (!D.isEmpty()) {                      if (!D.isEmpty()) {
2375                          add_EM_S=true;                          add_EM_S=true;
2376                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2377                          const double tmp0_1 = D_p[0]*w10;                          EM_S[INDEX2(0,0,4)]+=D_p[0]*w3;
2378                          EM_S[INDEX2(0,0,4)]+=tmp0_1;                          EM_S[INDEX2(1,0,4)]+=D_p[0]*w3;
2379                          EM_S[INDEX2(1,0,4)]+=tmp0_1;                          EM_S[INDEX2(2,0,4)]+=D_p[0]*w3;
2380                          EM_S[INDEX2(2,0,4)]+=tmp0_1;                          EM_S[INDEX2(3,0,4)]+=D_p[0]*w3;
2381                          EM_S[INDEX2(3,0,4)]+=tmp0_1;                          EM_S[INDEX2(0,1,4)]+=D_p[0]*w3;
2382                          EM_S[INDEX2(0,1,4)]+=tmp0_1;                          EM_S[INDEX2(1,1,4)]+=D_p[0]*w3;
2383                          EM_S[INDEX2(1,1,4)]+=tmp0_1;                          EM_S[INDEX2(2,1,4)]+=D_p[0]*w3;
2384                          EM_S[INDEX2(2,1,4)]+=tmp0_1;                          EM_S[INDEX2(3,1,4)]+=D_p[0]*w3;
2385                          EM_S[INDEX2(3,1,4)]+=tmp0_1;                          EM_S[INDEX2(0,2,4)]+=D_p[0]*w3;
2386                          EM_S[INDEX2(0,2,4)]+=tmp0_1;                          EM_S[INDEX2(1,2,4)]+=D_p[0]*w3;
2387                          EM_S[INDEX2(1,2,4)]+=tmp0_1;                          EM_S[INDEX2(2,2,4)]+=D_p[0]*w3;
2388                          EM_S[INDEX2(2,2,4)]+=tmp0_1;                          EM_S[INDEX2(3,2,4)]+=D_p[0]*w3;
2389                          EM_S[INDEX2(3,2,4)]+=tmp0_1;                          EM_S[INDEX2(0,3,4)]+=D_p[0]*w3;
2390                          EM_S[INDEX2(0,3,4)]+=tmp0_1;                          EM_S[INDEX2(1,3,4)]+=D_p[0]*w3;
2391                          EM_S[INDEX2(1,3,4)]+=tmp0_1;                          EM_S[INDEX2(2,3,4)]+=D_p[0]*w3;
2392                          EM_S[INDEX2(2,3,4)]+=tmp0_1;                          EM_S[INDEX2(3,3,4)]+=D_p[0]*w3;
                         EM_S[INDEX2(3,3,4)]+=tmp0_1;  
2393                      }                      }
2394                      ///////////////                      ///////////////
2395                      // process X //                      // process X //
# Line 2678  void Rectangle::assemblePDESingleReduced Line 2397  void Rectangle::assemblePDESingleReduced
2397                      if (!X.isEmpty()) {                      if (!X.isEmpty()) {
2398                          add_EM_F=true;                          add_EM_F=true;
2399                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2400                          const double tmp0_1 = X_p[0]*w11;                          const double tmp0 = 4*X_p[0]*w2;
2401                          const double tmp2_1 = X_p[0]*w13;                          const double tmp1 = 4*X_p[1]*w1;
2402                          const double tmp1_1 = X_p[1]*w12;                          EM_F[0]+=-tmp0 - tmp1;
2403                          const double tmp3_1 = X_p[1]*w14;                          EM_F[1]+=-tmp1 + tmp0;
2404                          EM_F[0]+=tmp0_1 + tmp1_1;                          EM_F[2]+=-tmp0 + tmp1;
2405                          EM_F[1]+=tmp1_1 + tmp2_1;                          EM_F[3]+= tmp0 + tmp1;
                         EM_F[2]+=tmp0_1 + tmp3_1;  
                         EM_F[3]+=tmp2_1 + tmp3_1;  
2406                      }                      }
2407                      ///////////////                      ///////////////
2408                      // process Y //                      // process Y //
# Line 2693  void Rectangle::assemblePDESingleReduced Line 2410  void Rectangle::assemblePDESingleReduced
2410                      if (!Y.isEmpty()) {                      if (!Y.isEmpty()) {
2411                          add_EM_F=true;                          add_EM_F=true;
2412                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2413                          const double tmp0_1 = Y_p[0]*w15;                          EM_F[0]+=4*Y_p[0]*w3;
2414                          EM_F[0]+=tmp0_1;                          EM_F[1]+=4*Y_p[0]*w3;
2415                          EM_F[1]+=tmp0_1;                          EM_F[2]+=4*Y_p[0]*w3;
2416                          EM_F[2]+=tmp0_1;                          EM_F[3]+=4*Y_p[0]*w3;
                         EM_F[3]+=tmp0_1;  
2417                      }                      }
2418    
2419                      // 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)
# Line 2723  void Rectangle::assemblePDESystem(Paso_S Line 2439  void Rectangle::assemblePDESystem(Paso_S
2439          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
2440      }      }
2441    
2442        /* GENERATOR SNIP_PDE_SYSTEM_PRE TOP */
2443      const double w0 = -0.1555021169820365539*m_dx[1]/m_dx[0];      const double w0 = -0.1555021169820365539*m_dx[1]/m_dx[0];
2444      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
2445      const double w2 = -0.15550211698203655390;      const double w2 = -0.15550211698203655390;
2446      const double w3 = 0.041666666666666666667*m_dx[0]/m_dx[1];      const double w3 = 0.041666666666666666667*m_dx[0]/m_dx[1];
2447      const double w4 = 0.15550211698203655390;      const double w4 = -0.01116454968463011277*m_dx[1]/m_dx[0];
2448      const double w5 = -0.041666666666666666667;      const double w5 = 0.011164549684630112770;
2449      const double w6 = -0.01116454968463011277*m_dx[1]/m_dx[0];      const double w6 = -0.041666666666666666667*m_dx[1]/m_dx[0];
2450      const double w7 = 0.011164549684630112770;      const double w7 = 0.1555021169820365539*m_dx[0]/m_dx[1];
2451      const double w8 = -0.011164549684630112770;      const double w8 = 0.01116454968463011277*m_dx[0]/m_dx[1];
2452      const double w9 = -0.041666666666666666667*m_dx[1]/m_dx[0];      const double w9 = -0.25;
2453      const double w10 = -0.041666666666666666667*m_dx[0]/m_dx[1];      const double w10 = -0.16666666666666666667*m_dx[0]/m_dx[1];
2454      const double w11 = 0.1555021169820365539*m_dx[1]/m_dx[0];      const double w11 = -0.032861463941450536761*m_dx[1];
2455      const double w12 = 0.1555021169820365539*m_dx[0]/m_dx[1];      const double w12 = -0.032861463941450536761*m_dx[0];
2456      const double w13 = 0.01116454968463011277*m_dx[0]/m_dx[1];      const double w13 = -0.12264065304058601714*m_dx[1];
2457      const double w14 = 0.01116454968463011277*m_dx[1]/m_dx[0];      const double w14 = -0.0023593469594139828636*m_dx[1];
2458      const double w15 = 0.041666666666666666667*m_dx[1]/m_dx[0];      const double w15 = -0.008805202725216129906*m_dx[0];
2459      const double w16 = -0.01116454968463011277*m_dx[0]/m_dx[1];      const double w16 = -0.008805202725216129906*m_dx[1];
2460      const double w17 = -0.1555021169820365539*m_dx[0]/m_dx[1];      const double w17 = -0.12264065304058601714*m_dx[0];
2461      const double w18 = -0.33333333333333333333*m_dx[1]/m_dx[0];      const double w18 = -0.0023593469594139828636*m_dx[0];
2462      const double w19 = 0.25000000000000000000;      const double w19 = -0.16666666666666666667*m_dx[1];
2463      const double w20 = -0.25000000000000000000;      const double w20 = -0.083333333333333333333*m_dx[0];
2464      const double w21 = 0.16666666666666666667*m_dx[0]/m_dx[1];      const double w21 = 0.025917019497006092316*m_dx[0]*m_dx[1];
2465      const double w22 = -0.16666666666666666667*m_dx[1]/m_dx[0];      const double w22 = 0.0018607582807716854616*m_dx[0]*m_dx[1];
2466      const double w23 = -0.16666666666666666667*m_dx[0]/m_dx[1];      const double w23 = 0.0069444444444444444444*m_dx[0]*m_dx[1];
2467      const double w24 = 0.33333333333333333333*m_dx[1]/m_dx[0];      const double w24 = 0.09672363354357992482*m_dx[0]*m_dx[1];
2468      const double w25 = 0.33333333333333333333*m_dx[0]/m_dx[1];      const double w25 = 0.000498588678642297402*m_dx[0]*m_dx[1];
2469      const double w26 = 0.16666666666666666667*m_dx[1]/m_dx[0];      const double w26 = 0.19716878364870322056*m_dx[1];
2470      const double w27 = -0.33333333333333333333*m_dx[0]/m_dx[1];      const double w27 = 0.052831216351296779436*m_dx[1];
2471      const double w28 = -0.032861463941450536761*m_dx[1];      const double w28 = 0.19716878364870322056*m_dx[0];
2472      const double w29 = -0.032861463941450536761*m_dx[0];      const double w29 = 0.052831216351296779436*m_dx[0];
2473      const double w30 = -0.12264065304058601714*m_dx[1];      /* GENERATOR SNIP_PDE_SYSTEM_PRE BOTTOM */
     const double w31 = -0.0023593469594139828636*m_dx[1];  
     const double w32 = -0.008805202725216129906*m_dx[0];  
     const double w33 = -0.008805202725216129906*m_dx[1];  
     const double w34 = 0.032861463941450536761*m_dx[1];  
     const double w35 = 0.008805202725216129906*m_dx[1];  
     const double w36 = 0.008805202725216129906*m_dx[0];  
     const double w37 = 0.0023593469594139828636*m_dx[1];  
     const double w38 = 0.12264065304058601714*m_dx[1];  
     const double w39 = 0.032861463941450536761*m_dx[0];  
     const double w40 = -0.12264065304058601714*m_dx[0];  
     const double w41 = -0.0023593469594139828636*m_dx[0];  
     const double w42 = 0.0023593469594139828636*m_dx[0];  
     const double w43 = 0.12264065304058601714*m_dx[0];  
     const double w44 = -0.16666666666666666667*m_dx[1];  
     const double w45 = -0.083333333333333333333*m_dx[0];  
     const double w46 = 0.083333333333333333333*m_dx[1];  
     const double w47 = 0.16666666666666666667*m_dx[1];  
     const double w48 = 0.083333333333333333333*m_dx[0];  
     const double w49 = -0.16666666666666666667*m_dx[0];  
     const double w50 = 0.16666666666666666667*m_dx[0];  
     const double w51 = -0.083333333333333333333*m_dx[1];  
     const double w52 = 0.025917019497006092316*m_dx[0]*m_dx[1];  
     const double w53 = 0.0018607582807716854616*m_dx[0]*m_dx[1];  
     const double w54 = 0.0069444444444444444444*m_dx[0]*m_dx[1];  
     const double w55 = 0.09672363354357992482*m_dx[0]*m_dx[1];  
     const double w56 = 0.00049858867864229740201*m_dx[0]*m_dx[1];  
     const double w57 = 0.055555555555555555556*m_dx[0]*m_dx[1];  
     const double w58 = 0.027777777777777777778*m_dx[0]*m_dx[1];  
     const double w59 = 0.11111111111111111111*m_dx[0]*m_dx[1];  
     const double w60 = -0.19716878364870322056*m_dx[1];  
     const double w61 = -0.19716878364870322056*m_dx[0];  
     const double w62 = -0.052831216351296779436*m_dx[0];  
     const double w63 = -0.052831216351296779436*m_dx[1];  
     const double w64 = 0.19716878364870322056*m_dx[1];  
     const double w65 = 0.052831216351296779436*m_dx[1];  
     const double w66 = 0.19716878364870322056*m_dx[0];  
     const double w67 = 0.052831216351296779436*m_dx[0];  
     const double w68 = -0.5*m_dx[1];  
     const double w69 = -0.5*m_dx[0];  
     const double w70 = 0.5*m_dx[1];  
     const double w71 = 0.5*m_dx[0];  
     const double w72 = 0.1555021169820365539*m_dx[0]*m_dx[1];  
     const double w73 = 0.041666666666666666667*m_dx[0]*m_dx[1];  
     const double w74 = 0.01116454968463011277*m_dx[0]*m_dx[1];  
     const double w75 = 0.25*m_dx[0]*m_dx[1];  
2474    
2475      rhs.requireWrite();      rhs.requireWrite();
2476  #pragma omp parallel  #pragma omp parallel
# Line 2812  void Rectangle::assemblePDESystem(Paso_S Line 2484  void Rectangle::assemblePDESystem(Paso_S
2484                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
2485                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
2486                      const index_t e = k0 + m_NE[0]*k1;                      const index_t e = k0 + m_NE[0]*k1;
2487                        /* GENERATOR SNIP_PDE_SYSTEM TOP */
2488                      ///////////////                      ///////////////
2489                      // process A //                      // process A //
2490                      ///////////////                      ///////////////
2491                      if (!A.isEmpty()) {                      if (!A.isEmpty()) {
2492                          add_EM_S=true;                          add_EM_S = true;
2493                          const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);                          const double* A_p = const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2494                          if (A.actsExpanded()) {                          if (A.actsExpanded()) {
2495                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2496                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2497                                      const double A_00_0 = A_p[INDEX5(k,0,m,0,0, numEq,2,numComp,2)];                                      const double A_00_0 = A_p[INDEX5(k,0,m,0,0,numEq,2,numComp,2)];
2498                                      const double A_01_0 = A_p[INDEX5(k,0,m,1,0, numEq,2,numComp,2)];                                      const double A_01_0 = A_p[INDEX5(k,0,m,1,0,numEq,2,numComp,2)];
2499                                      const double A_10_0 = A_p[INDEX5(k,1,m,0,0, numEq,2,numComp,2)];                                      const double A_10_0 = A_p[INDEX5(k,1,m,0,0,numEq,2,numComp,2)];
2500                                      const double A_11_0 = A_p[INDEX5(k,1,m,1,0, numEq,2,numComp,2)];                                      const double A_11_0 = A_p[INDEX5(k,1,m,1,0,numEq,2,numComp,2)];
2501                                      const double A_00_1 = A_p[INDEX5(k,0,m,0,1, numEq,2,numComp,2)];                                      const double A_00_1 = A_p[INDEX5(k,0,m,0,1,numEq,2,numComp,2)];
2502                                      const double A_01_1 = A_p[INDEX5(k,0,m,1,1, numEq,2,numComp,2)];                                      const double A_01_1 = A_p[INDEX5(k,0,m,1,1,numEq,2,numComp,2)];
2503                                      const double A_10_1 = A_p[INDEX5(k,1,m,0,1, numEq,2,numComp,2)];                                      const double A_10_1 = A_p[INDEX5(k,1,m,0,1,numEq,2,numComp,2)];
2504                                      const double A_11_1 = A_p[INDEX5(k,1,m,1,1, numEq,2,numComp,2)];                                      const double A_11_1 = A_p[INDEX5(k,1,m,1,1,numEq,2,numComp,2)];
2505                                      const double A_00_2 = A_p[INDEX5(k,0,m,0,2, numEq,2,numComp,2)];                                      const double A_00_2 = A_p[INDEX5(k,0,m,0,2,numEq,2,numComp,2)];
2506                                      const double A_01_2 = A_p[INDEX5(k,0,m,1,2, numEq,2,numComp,2)];                                      const double A_01_2 = A_p[INDEX5(k,0,m,1,2,numEq,2,numComp,2)];
2507                                      const double A_10_2 = A_p[INDEX5(k,1,m,0,2, numEq,2,numComp,2)];                                      const double A_10_2 = A_p[INDEX5(k,1,m,0,2,numEq,2,numComp,2)];
2508                                      const double A_11_2 = A_p[INDEX5(k,1,m,1,2, numEq,2,numComp,2)];                                      const double A_11_2 = A_p[INDEX5(k,1,m,1,2,numEq,2,numComp,2)];
2509                                      const double A_00_3 = A_p[INDEX5(k,0,m,0,3, numEq,2,numComp,2)];                                      const double A_00_3 = A_p[INDEX5(k,0,m,0,3,numEq,2,numComp,2)];
2510                                      const double A_01_3 = A_p[INDEX5(k,0,m,1,3, numEq,2,numComp,2)];                                      const double A_01_3 = A_p[INDEX5(k,0,m,1,3,numEq,2,numComp,2)];
2511                                      const double A_10_3 = A_p[INDEX5(k,1,m,0,3, numEq,2,numComp,2)];                                      const double A_10_3 = A_p[INDEX5(k,1,m,0,3,numEq,2,numComp,2)];
2512                                      const double A_11_3 = A_p[INDEX5(k,1,m,1,3, numEq,2,numComp,2)];                                      const double A_11_3 = A_p[INDEX5(k,1,m,1,3,numEq,2,numComp,2)];
2513                                      const double tmp0_0 = A_01_0 + A_01_3;                                      const double tmp0 = w3*(A_11_0 + A_11_1 + A_11_2 + A_11_3);
2514                                      const double tmp1_0 = A_00_0 + A_00_1;                                      const double tmp1 = w1*(A_01_0 + A_01_3 - A_10_1 - A_10_2);
2515                                      const double tmp2_0 = A_11_0 + A_11_1 + A_11_2 + A_11_3;                                      const double tmp2 = w4*(A_00_2 + A_00_3);
2516                                      const double tmp3_0 = A_00_2 + A_00_3;                                      const double tmp3 = w0*(A_00_0 + A_00_1);
2517                                      const double tmp4_0 = A_10_1 + A_10_2;                                      const double tmp4 = w5*(A_01_2 - A_10_3);
2518                                      const double tmp5_0 = A_00_0 + A_00_1 + A_00_2 + A_00_3;                                      const double tmp5 = w2*(-A_01_1 + A_10_0);
2519                                      const double tmp6_0 = A_01_3 + A_10_0;                                      const double tmp6 = w5*(A_01_3 + A_10_0);
2520                                      const double tmp7_0 = A_01_0 + A_10_3;                                      const double tmp7 = w3*(-A_11_0 - A_11_1 - A_11_2 - A_11_3);
2521                                      const double tmp8_0 = A_01_1 + A_01_2 + A_10_1 + A_10_2;                                      const double tmp8 = w6*(A_00_0 + A_00_1 + A_00_2 + A_00_3);
2522                                      const double tmp9_0 = A_01_0 + A_10_0;                                      const double tmp9 = w1*(A_01_1 + A_01_2 + A_10_1 + A_10_2);
2523                                      const double tmp10_0 = A_01_3 + A_10_3;                                      const double tmp10 = w2*(-A_01_0 - A_10_3);
2524                                      const double tmp11_0 = A_11_1 + A_11_3;                                      const double tmp11 = w4*(A_00_0 + A_00_1);
2525                                      const double tmp12_0 = A_11_0 + A_11_2;                                      const double tmp12 = w0*(A_00_2 + A_00_3);
2526                                      const double tmp13_0 = A_01_2 + A_10_1;                                      const double tmp13 = w5*(A_01_1 - A_10_0);
2527                                      const double tmp14_0 = A_01_0 + A_01_3 + A_10_0 + A_10_3;                                      const double tmp14 = w2*(-A_01_2 + A_10_3);
2528                                      const double tmp15_0 = A_01_1 + A_10_2;                                      const double tmp15 = w7*(A_11_0 + A_11_2);
2529                                      const double tmp16_0 = A_10_0 + A_10_3;                                      const double tmp16 = w4*(-A_00_2 - A_00_3);
2530                                      const double tmp17_0 = A_01_1 + A_01_2;                                      const double tmp17 = w0*(-A_00_0 - A_00_1);
2531                                      const double tmp18_0 = A_01_1 + A_10_1;                                      const double tmp18 = w5*(A_01_3 + A_10_3);
2532                                      const double tmp19_0 = A_01_2 + A_10_2;                                      const double tmp19 = w8*(A_11_1 + A_11_3);
2533                                      const double tmp0_1 = A_10_3*w8;                                      const double tmp20 = w2*(-A_01_0 - A_10_0);
2534                                      const double tmp1_1 = tmp0_0*w1;                                      const double tmp21 = w7*(A_11_1 + A_11_3);
2535                                      const double tmp2_1 = A_01_1*w4;                                      const double tmp22 = w4*(-A_00_0 - A_00_1);
2536                                      const double tmp3_1 = tmp1_0*w0;                                      const double tmp23 = w0*(-A_00_2 - A_00_3);
2537                                      const double tmp4_1 = A_01_2*w7;                                      const double tmp24 = w5*(A_01_0 + A_10_0);
2538                                      const double tmp5_1 = tmp2_0*w3;                                      const double tmp25 = w8*(A_11_0 + A_11_2);
2539                                      const double tmp6_1 = tmp3_0*w6;                                      const double tmp26 = w2*(-A_01_3 - A_10_3);
2540                                      const double tmp7_1 = A_10_0*w2;                                      const double tmp27 = w5*(-A_01_1 - A_10_2);
2541                                      const double tmp8_1 = tmp4_0*w5;                                      const double tmp28 = w1*(-A_01_0 - A_01_3 - A_10_0 - A_10_3);
2542                                      const double tmp9_1 = tmp2_0*w10;                                      const double tmp29 = w2*(A_01_2 + A_10_1);
2543                                      const double tmp10_1 = tmp5_0*w9;                                      const double tmp30 = w7*(-A_11_1 - A_11_3);
2544                                      const double tmp11_1 = tmp6_0*w7;                                      const double tmp31 = w1*(-A_01_1 - A_01_2 + A_10_0 + A_10_3);
2545                                      const double tmp12_1 = tmp7_0*w4;                                      const double tmp32 = w5*(-A_01_0 + A_10_2);
2546                                      const double tmp13_1 = tmp8_0*w1;                                      const double tmp33 = w8*(-A_11_0 - A_11_2);
2547                                      const double tmp14_1 = A_10_0*w8;                                      const double tmp34 = w6*(-A_00_0 - A_00_1 - A_00_2 - A_00_3);
2548                                      const double tmp15_1 = A_01_2*w4;                                      const double tmp35 = w2*(A_01_3 - A_10_1);
2549                                      const double tmp16_1 = tmp3_0*w0;                                      const double tmp36 = w5*(A_01_0 + A_10_3);
2550                                      const double tmp17_1 = A_01_1*w7;                                      const double tmp37 = w2*(-A_01_3 - A_10_0);
2551                                      const double tmp18_1 = tmp1_0*w6;                                      const double tmp38 = w7*(-A_11_0 - A_11_2);
2552                                      const double tmp19_1 = A_10_3*w2;                                      const double tmp39 = w5*(-A_01_3 + A_10_1);
2553                                      const double tmp20_1 = tmp9_0*w4;                                      const double tmp40 = w8*(-A_11_1 - A_11_3);
2554                                      const double tmp21_1 = tmp1_0*w11;                                      const double tmp41 = w2*(A_01_0 - A_10_2);
2555                                      const double tmp22_1 = tmp10_0*w7;                                      const double tmp42 = w5*(A_01_1 - A_10_3);
2556                                      const double tmp23_1 = tmp3_0*w14;                                      const double tmp43 = w2*(-A_01_2 + A_10_0);
2557                                      const double tmp24_1 = tmp11_0*w13;                                      const double tmp44 = w5*(A_01_2 - A_10_0);
2558                                      const double tmp25_1 = tmp12_0*w12;                                      const double tmp45 = w2*(-A_01_1 + A_10_3);
2559                                      const double tmp26_1 = tmp10_0*w4;                                      const double tmp46 = w5*(-A_01_0 + A_10_1);
2560                                      const double tmp27_1 = tmp3_0*w11;                                      const double tmp47 = w2*(A_01_3 - A_10_2);
2561                                      const double tmp28_1 = tmp9_0*w7;                                      const double tmp48 = w5*(-A_01_1 - A_10_1);
2562                                      const double tmp29_1 = tmp1_0*w14;                                      const double tmp49 = w2*(A_01_2 + A_10_2);
2563                                      const double tmp30_1 = tmp12_0*w13;                                      const double tmp50 = w5*(-A_01_3 + A_10_2);
2564                                      const double tmp31_1 = tmp11_0*w12;                                      const double tmp51 = w2*(A_01_0 - A_10_1);
2565                                      const double tmp32_1 = tmp13_0*w2;                                      const double tmp52 = w5*(-A_01_2 - A_10_1);
2566                                      const double tmp33_1 = tmp14_0*w5;                                      const double tmp53 = w2*(A_01_1 + A_10_2);
2567                                      const double tmp34_1 = tmp15_0*w8;                                      const double tmp54 = w5*(-A_01_2 - A_10_2);
2568                                      const double tmp35_1 = A_01_0*w8;                                      const double tmp55 = w2*(A_01_1 + A_10_1);
2569                                      const double tmp36_1 = tmp16_0*w1;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp15 + tmp16 + tmp17 + tmp18 + tmp19 + tmp20 + tmp9;
2570                                      const double tmp37_1 = A_10_1*w4;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5;
2571                                      const double tmp38_1 = tmp5_0*w15;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp31 + tmp34 + tmp38 + tmp39 + tmp40 + tmp41;
2572                                      const double tmp39_1 = A_10_2*w7;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp28 + tmp52 + tmp53 + tmp7 + tmp8;
2573                                      const double tmp40_1 = tmp11_0*w17;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0 + tmp2 + tmp3 + tmp31 + tmp50 + tmp51;
2574                                      const double tmp41_1 = A_01_3*w2;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp16 + tmp17 + tmp21 + tmp25 + tmp28 + tmp54 + tmp55;
2575                                      const double tmp42_1 = tmp12_0*w16;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10 + tmp6 + tmp7 + tmp8 + tmp9;
2576                                      const double tmp43_1 = tmp17_0*w5;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1 + tmp30 + tmp33 + tmp34 + tmp44 + tmp45;
2577                                      const double tmp44_1 = tmp7_0*w7;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1 + tmp34 + tmp38 + tmp40 + tmp42 + tmp43;
2578                                      const double tmp45_1 = tmp6_0*w4;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp36 + tmp37 + tmp7 + tmp8 + tmp9;
2579                                      const double tmp46_1 = A_01_3*w8;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp15 + tmp19 + tmp22 + tmp23 + tmp28 + tmp48 + tmp49;
2580                                      const double tmp47_1 = A_10_2*w4;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0 + tmp11 + tmp12 + tmp31 + tmp46 + tmp47;
2581                                      const double tmp48_1 = A_10_1*w7;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp27 + tmp28 + tmp29 + tmp7 + tmp8;
2582                                      const double tmp49_1 = tmp12_0*w17;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp30 + tmp31 + tmp32 + tmp33 + tmp34 + tmp35;
2583                                      const double tmp50_1 = A_01_0*w2;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0 + tmp1 + tmp11 + tmp12 + tmp13 + tmp14;
2584                                      const double tmp51_1 = tmp11_0*w16;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp21 + tmp22 + tmp23 + tmp24 + tmp25 + tmp26 + tmp9;
                                     const double tmp52_1 = tmp18_0*w8;  
                                     const double tmp53_1 = tmp19_0*w2;  
                                     const double tmp54_1 = tmp13_0*w8;  
                                     const double tmp55_1 = tmp15_0*w2;  
                                     const double tmp56_1 = tmp19_0*w8;  
                                     const double tmp57_1 = tmp18_0*w2;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp13_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp36_1 + tmp37_1 + tmp39_1 + tmp3_1 + tmp43_1 + tmp46_1 + tmp50_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp15_1 + tmp17_1 + tmp1_1 + tmp38_1 + tmp49_1 + tmp51_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp10_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp21_1 + tmp23_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp56_1 + tmp57_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp10_1 + tmp13_1 + tmp44_1 + tmp45_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp36_1 + tmp38_1 + tmp43_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp27_1 + tmp29_1 + tmp33_1 + tmp52_1 + tmp53_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp5_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp10_1 + tmp33_1 + tmp54_1 + tmp55_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp14_1 + tmp19_1 + tmp1_1 + tmp2_1 + tmp38_1 + tmp40_1 + tmp42_1 + tmp4_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp16_1 + tmp18_1 + tmp35_1 + tmp36_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp48_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp13_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
2585                                  }                                  }
2586                              }                              }
2587                          } else { // constant data                          } else { // constant data
# Line 2940  void Rectangle::assemblePDESystem(Paso_S Line 2591  void Rectangle::assemblePDESystem(Paso_S
2591                                      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)];
2592                                      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)];
2593                                      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)];
2594                                      const double tmp0_0 = A_01 + A_10;                                      const double tmp0 = w9*(-A_01 - A_10);
2595                                      const double tmp0_1 = A_00*w18;                                      const double tmp1 = w9*(-A_01 + A_10);
2596                                      const double tmp1_1 = A_01*w19;                                      const double tmp2 = w9*(A_01 + A_10);
2597                                      const double tmp2_1 = A_10*w20;                                      const double tmp3 = w9*(A_01 - A_10);
2598                                      const double tmp3_1 = A_11*w21;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;
2599                                      const double tmp4_1 = A_00*w22;                                      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;
2600                                      const double tmp5_1 = tmp0_0*w19;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;
2601                                      const double tmp6_1 = A_11*w23;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp2;
2602                                      const double tmp7_1 = A_11*w25;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=8*A_00*w6 - A_11*w10 + tmp3;
2603                                      const double tmp8_1 = A_00*w24;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;
2604                                      const double tmp9_1 = tmp0_0*w20;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp0;
2605                                      const double tmp10_1 = A_01*w20;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;
2606                                      const double tmp11_1 = A_11*w27;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp1;
2607                                      const double tmp12_1 = A_00*w26;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp0;
2608                                      const double tmp13_1 = A_10*w19;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp2;
2609                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=8*A_00*w6 - A_11*w10 + tmp3;
2610                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=4*A_00*w6 + A_11*w10 + tmp2;
2611                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-4*A_00*w6 + 2*A_11*w10 + tmp3;
2612                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=8*A_00*w6 - A_11*w10 + tmp1;
2613                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-8*A_00*w6 - 2*A_11*w10 + tmp0;
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp4_1 + tmp6_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp3_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp7_1 + tmp8_1;  
2614                                  }                                  }
2615                              }                              }
2616                          }                          }
# Line 2992  void Rectangle::assemblePDESystem(Paso_S Line 2632  void Rectangle::assemblePDESystem(Paso_S
2632                                      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)];
2633                                      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)];
2634                                      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)];
2635                                      const double tmp0_0 = B_1_0 + B_1_1;                                      const double tmp0 = w15*(B_1_2 + B_1_3);
2636                                      const double tmp1_0 = B_1_2 + B_1_3;                                      const double tmp1 = w12*(B_1_0 + B_1_1);
2637                                      const double tmp2_0 = B_0_1 + B_0_3;                                      const double tmp2 = w15*(B_1_0 + B_1_1);
2638                                      const double tmp3_0 = B_0_0 + B_0_2;                                      const double tmp3 = w16*(-B_0_1 - B_0_3);
2639                                      const double tmp63_1 = B_1_1*w42;                                      const double tmp4 = w11*(-B_0_0 - B_0_2);
2640                                      const double tmp79_1 = B_1_1*w40;                                      const double tmp5 = w12*(B_1_2 + B_1_3);
2641                                      const double tmp37_1 = tmp3_0*w35;                                      const double tmp6 = w15*(-B_1_0 - B_1_1);
2642                                      const double tmp8_1 = tmp0_0*w32;                                      const double tmp7 = w12*(-B_1_2 - B_1_3);
2643                                      const double tmp71_1 = B_0_1*w34;                                      const double tmp8 = w15*(-B_1_2 - B_1_3);
2644                                      const double tmp19_1 = B_0_3*w31;                                      const double tmp9 = w12*(-B_1_0 - B_1_1);
2645                                      const double tmp15_1 = B_0_3*w34;                                      const double tmp10 = w11*(-B_0_1 - B_0_3);
2646                                      const double tmp9_1 = tmp3_0*w34;                                      const double tmp11 = w16*(-B_0_0 - B_0_2);
2647                                      const double tmp35_1 = B_1_0*w36;                                      const double tmp12 = w16*(B_0_0 + B_0_2);
2648                                      const double tmp66_1 = B_0_3*w28;                                      const double tmp13 = w11*(B_0_1 + B_0_3);
2649                                      const double tmp28_1 = B_1_0*w42;                                      const double tmp14 = w11*(B_0_0 + B_0_2);
2650                                      const double tmp22_1 = B_1_0*w40;                                      const double tmp15 = w16*(B_0_1 + B_0_3);
2651                                      const double tmp16_1 = B_1_2*w29;                                      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;
2652                                      const double tmp6_1 = tmp2_0*w35;                                      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;
2653                                      const double tmp55_1 = B_1_3*w40;                                      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;
2654                                      const double tmp50_1 = B_1_3*w42;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp12 + tmp13 + tmp2 + tmp5;
2655                                      const double tmp7_1 = tmp1_0*w29;                                      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;
2656                                      const double tmp1_1 = tmp1_0*w32;                                      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;
2657                                      const double tmp57_1 = B_0_3*w30;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2658                                      const double tmp18_1 = B_1_1*w32;                                      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;
2659                                      const double tmp53_1 = B_1_0*w41;                                      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;
2660                                      const double tmp61_1 = B_1_3*w36;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2661                                      const double tmp27_1 = B_0_3*w38;                                      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;
2662                                      const double tmp64_1 = B_0_2*w30;                                      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;
2663                                      const double tmp76_1 = B_0_1*w38;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp3 + tmp4 + tmp8 + tmp9;
2664                                      const double tmp39_1 = tmp2_0*w34;                                      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;
2665                                      const double tmp62_1 = B_0_1*w31;                                      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;
2666                                      const double tmp56_1 = B_0_0*w31;                                      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;
                                     const double tmp49_1 = B_1_1*w36;  
                                     const double tmp2_1 = B_0_2*w31;  
                                     const double tmp23_1 = B_0_2*w33;  
                                     const double tmp38_1 = B_1_1*w43;  
                                     const double tmp74_1 = B_1_2*w41;  
                                     const double tmp43_1 = B_1_1*w41;  
                                     const double tmp58_1 = B_0_2*w28;  
                                     const double tmp67_1 = B_0_0*w33;  
                                     const double tmp33_1 = tmp0_0*w39;  
                                     const double tmp4_1 = B_0_0*w28;  
                                     const double tmp20_1 = B_0_0*w30;  
                                     const double tmp13_1 = B_0_2*w38;  
                                     const double tmp65_1 = B_1_2*w43;  
                                     const double tmp0_1 = tmp0_0*w29;  
                                     const double tmp41_1 = tmp3_0*w33;  
                                     const double tmp73_1 = B_0_2*w37;  
                                     const double tmp69_1 = B_0_0*w38;  
                                     const double tmp48_1 = B_1_2*w39;  
                                     const double tmp59_1 = B_0_1*w33;  
                                     const double tmp17_1 = B_1_3*w41;  
                                     const double tmp5_1 = B_0_3*w33;  
                                     const double tmp3_1 = B_0_1*w30;  
                                     const double tmp21_1 = B_0_1*w28;  
                                     const double tmp42_1 = B_1_0*w29;  
                                     const double tmp54_1 = B_1_2*w32;  
                                     const double tmp60_1 = B_1_0*w39;  
                                     const double tmp32_1 = tmp1_0*w36;  
                                     const double tmp10_1 = B_0_1*w37;  
                                     const double tmp14_1 = B_0_0*w35;  
                                     const double tmp29_1 = B_0_1*w35;  
                                     const double tmp26_1 = B_1_2*w36;  
                                     const double tmp30_1 = B_1_3*w43;  
                                     const double tmp70_1 = B_0_2*w35;  
                                     const double tmp34_1 = B_1_3*w39;  
                                     const double tmp51_1 = B_1_0*w43;  
                                     const double tmp31_1 = B_0_2*w34;  
                                     const double tmp45_1 = tmp3_0*w28;  
                                     const double tmp11_1 = tmp1_0*w39;  
                                     const double tmp52_1 = B_1_1*w29;  
                                     const double tmp44_1 = B_1_3*w32;  
                                     const double tmp25_1 = B_1_1*w39;  
                                     const double tmp47_1 = tmp2_0*w33;  
                                     const double tmp72_1 = B_1_3*w29;  
                                     const double tmp40_1 = tmp2_0*w28;  
                                     const double tmp46_1 = B_1_2*w40;  
                                     const double tmp36_1 = B_1_2*w42;  
                                     const double tmp24_1 = B_0_0*w37;  
                                     const double tmp77_1 = B_0_3*w35;  
                                     const double tmp68_1 = B_0_3*w37;  
                                     const double tmp78_1 = B_0_0*w34;  
                                     const double tmp12_1 = tmp0_0*w36;  
                                     const double tmp75_1 = B_1_0*w32;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp45_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp6_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp7_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp39_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp11_1 + tmp12_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
2667                                  }                                  }
2668                              }                              }
2669                          } else { // constant data                          } else { // constant data
2670                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2671                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2672                                      const double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];                                      const double B_0 = B_p[INDEX3(k,0,m,numEq,2)];
2673                                      const double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];                                      const double B_1 = B_p[INDEX3(k,1,m,numEq,2)];
2674                                      const double tmp6_1 = B_1*w50;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=B_0*w19 + 2*B_1*w20;
2675                                      const double tmp1_1 = B_1*w45;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=B_0*w19 + B_1*w20;
2676                                      const double tmp5_1 = B_1*w49;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=B_0*w19/2 + 2*B_1*w20;
2677                                      const double tmp4_1 = B_1*w48;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=B_0*w19/2 + B_1*w20;
2678                                      const double tmp0_1 = B_0*w44;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-B_0*w19 + B_1*w20;
2679                                      const double tmp2_1 = B_0*w46;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-B_0*w19 + 2*B_1*w20;
2680                                      const double tmp7_1 = B_0*w51;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=-B_0*w19/2 + B_1*w20;
2681                                      const double tmp3_1 = B_0*w47;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-B_0*w19/2 + 2*B_1*w20;
2682                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp5_1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=B_0*w19/2 - 2*B_1*w20;
2683                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=B_0*w19/2 - B_1*w20;
2684                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp6_1 + tmp7_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=B_0*w19 - 2*B_1*w20;
2685                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=B_0*w19 - B_1*w20;
2686                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=-B_0*w19/2 - B_1*w20;
2687                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp3_1 + tmp5_1;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-B_0*w19/2 - 2*B_1*w20;
2688                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp4_1 + tmp7_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=-B_0*w19 - B_1*w20;
2689                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp6_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-B_0*w19 - 2*B_1*w20;
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp3_1 + tmp6_1;  
2690                                  }                                  }
2691                              }                              }
2692                          }                          }
# Line 3144  void Rectangle::assemblePDESystem(Paso_S Line 2708  void Rectangle::assemblePDESystem(Paso_S
2708                                      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)];
2709                                      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)];
2710                                      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)];
2711                                      const double tmp2_0 = C_0_1 + C_0_3;                                      const double tmp0 = w15*(C_1_2 + C_1_3);
2712                                      const double tmp1_0 = C_1_2 + C_1_3;                                      const double tmp1 = w12*(C_1_0 + C_1_1);
2713                                      const double tmp3_0 = C_0_0 + C_0_2;                                      const double tmp2 = w15*(-C_1_2 - C_1_3);
2714                                      const double tmp0_0 = C_1_0 + C_1_1;                                      const double tmp3 = w16*(C_0_0 + C_0_2);
2715                                      const double tmp64_1 = C_0_2*w30;                                      const double tmp4 = w11*(C_0_1 + C_0_3);
2716                                      const double tmp14_1 = C_0_2*w28;                                      const double tmp5 = w12*(-C_1_0 - C_1_1);
2717                                      const double tmp19_1 = C_0_3*w31;                                      const double tmp6 = w15*(-C_1_0 - C_1_1);
2718                                      const double tmp22_1 = C_1_0*w40;                                      const double tmp7 = w12*(-C_1_2 - C_1_3);
2719                                      const double tmp37_1 = tmp3_0*w35;                                      const double tmp8 = w15*(C_1_0 + C_1_1);
2720                                      const double tmp29_1 = C_0_1*w35;                                      const double tmp9 = w12*(C_1_2 + C_1_3);
2721                                      const double tmp73_1 = C_0_2*w37;                                      const double tmp10 = w11*(-C_0_1 - C_0_3);
2722                                      const double tmp74_1 = C_1_2*w41;                                      const double tmp11 = w16*(-C_0_0 - C_0_2);
2723                                      const double tmp52_1 = C_1_3*w39;                                      const double tmp12 = w16*(-C_0_1 - C_0_3);
2724                                      const double tmp25_1 = C_1_1*w39;                                      const double tmp13 = w11*(-C_0_0 - C_0_2);
2725                                      const double tmp62_1 = C_0_1*w31;                                      const double tmp14 = w11*(C_0_0 + C_0_2);
2726                                      const double tmp79_1 = C_1_1*w40;                                      const double tmp15 = w16*(C_0_1 + C_0_3);
2727                                      const double tmp43_1 = C_1_1*w36;                                      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;
2728                                      const double tmp27_1 = C_0_3*w38;                                      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;
2729                                      const double tmp28_1 = C_1_0*w42;                                      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;
2730                                      const double tmp63_1 = C_1_1*w42;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp12 + tmp13 + tmp2 + tmp5;
2731                                      const double tmp59_1 = C_0_3*w34;                                      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;
2732                                      const double tmp72_1 = C_1_3*w29;                                      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;
2733                                      const double tmp40_1 = tmp2_0*w35;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2 + tmp3 + tmp4 + tmp5;
2734                                      const double tmp13_1 = C_0_3*w30;                                      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;
2735                                      const double tmp51_1 = C_1_2*w40;                                      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;
2736                                      const double tmp54_1 = C_1_2*w42;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp12 + tmp13 + tmp8 + tmp9;
2737                                      const double tmp12_1 = C_0_0*w31;                                      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;
2738                                      const double tmp2_1 = tmp1_0*w32;                                      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;
2739                                      const double tmp68_1 = C_0_2*w31;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp3 + tmp4 + tmp8 + tmp9;
2740                                      const double tmp75_1 = C_1_0*w32;                                      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;
2741                                      const double tmp49_1 = C_1_1*w41;                                      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;
2742                                      const double tmp4_1 = C_0_2*w35;                                      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;
                                     const double tmp66_1 = C_0_3*w28;  
                                     const double tmp56_1 = C_0_1*w37;  
                                     const double tmp5_1 = C_0_1*w34;  
                                     const double tmp38_1 = tmp2_0*w34;  
                                     const double tmp76_1 = C_0_1*w38;  
                                     const double tmp21_1 = C_0_1*w28;  
                                     const double tmp69_1 = C_0_1*w30;  
                                     const double tmp53_1 = C_1_0*w36;  
                                     const double tmp42_1 = C_1_2*w39;  
                                     const double tmp32_1 = tmp1_0*w29;  
                                     const double tmp45_1 = C_1_0*w43;  
                                     const double tmp33_1 = tmp0_0*w32;  
                                     const double tmp35_1 = C_1_0*w41;  
                                     const double tmp26_1 = C_1_2*w36;  
                                     const double tmp67_1 = C_0_0*w33;  
                                     const double tmp31_1 = C_0_2*w34;  
                                     const double tmp20_1 = C_0_0*w30;  
                                     const double tmp70_1 = C_0_0*w28;  
                                     const double tmp8_1 = tmp0_0*w39;  
                                     const double tmp30_1 = C_1_3*w43;  
                                     const double tmp0_1 = tmp0_0*w29;  
                                     const double tmp17_1 = C_1_3*w41;  
                                     const double tmp58_1 = C_0_0*w35;  
                                     const double tmp9_1 = tmp3_0*w33;  
                                     const double tmp61_1 = C_1_3*w36;  
                                     const double tmp41_1 = tmp3_0*w34;  
                                     const double tmp50_1 = C_1_3*w32;  
                                     const double tmp18_1 = C_1_1*w32;  
                                     const double tmp6_1 = tmp1_0*w36;  
                                     const double tmp3_1 = C_0_0*w38;  
                                     const double tmp34_1 = C_1_1*w29;  
                                     const double tmp77_1 = C_0_3*w35;  
                                     const double tmp65_1 = C_1_2*w43;  
                                     const double tmp71_1 = C_0_3*w33;  
                                     const double tmp55_1 = C_1_1*w43;  
                                     const double tmp46_1 = tmp3_0*w28;  
                                     const double tmp24_1 = C_0_0*w37;  
                                     const double tmp10_1 = tmp1_0*w39;  
                                     const double tmp48_1 = C_1_0*w29;  
                                     const double tmp15_1 = C_0_1*w33;  
                                     const double tmp36_1 = C_1_2*w32;  
                                     const double tmp60_1 = C_1_0*w39;  
                                     const double tmp47_1 = tmp2_0*w33;  
                                     const double tmp16_1 = C_1_2*w29;  
                                     const double tmp1_1 = C_0_3*w37;  
                                     const double tmp7_1 = tmp2_0*w28;  
                                     const double tmp39_1 = C_1_3*w40;  
                                     const double tmp44_1 = C_1_3*w42;  
                                     const double tmp57_1 = C_0_2*w38;  
                                     const double tmp78_1 = C_0_0*w34;  
                                     const double tmp11_1 = tmp0_0*w36;  
                                     const double tmp23_1 = C_0_2*w33;  
                                     EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp7_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp32_1 + tmp33_1 + tmp40_1 + tmp41_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp37_1 + tmp38_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp10_1 + tmp11_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp2_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp40_1 + tmp41_1 + tmp6_1 + tmp8_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1;  
2743                                  }                                  }
2744                              }                              }
2745                          } else { // constant data                          } else { // constant data
2746                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2747                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2748                                      const double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];                                      const double C_0 = C_p[INDEX3(k,m,0,numEq,numComp)];
2749                                      const double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];                                      const double C_1 = C_p[INDEX3(k,m,1,numEq,numComp)];
2750                                      const double tmp1_1 = C_1*w45;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=C_0*w19 + 2*C_1*w20;
2751                                      const double tmp3_1 = C_0*w51;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-C_0*w19 + C_1*w20;
2752                                      const double tmp4_1 = C_0*w44;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=C_0*w19/2 - 2*C_1*w20;
2753                                      const double tmp7_1 = C_0*w46;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=-C_0*w19/2 - C_1*w20;
2754                                      const double tmp5_1 = C_1*w49;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=C_0*w19 + C_1*w20;
2755                                      const double tmp2_1 = C_1*w48;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-C_0*w19 + 2*C_1*w20;
2756                                      const double tmp0_1 = C_0*w47;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=C_0*w19/2 - C_1*w20;
2757                                      const double tmp6_1 = C_1*w50;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-C_0*w19/2 - 2*C_1*w20;
2758                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=C_0*w19/2 + 2*C_1*w20;
2759                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=-C_0*w19/2 + C_1*w20;
2760                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp4_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=C_0*w19 - 2*C_1*w20;
2761                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp4_1 + tmp5_1;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-C_0*w19 - C_1*w20;
2762                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp6_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=C_0*w19/2 + C_1*w20;
2763                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-C_0*w19/2 + 2*C_1*w20;
2764                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp5_1 + tmp7_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=C_0*w19 - C_1*w20;
2765                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1 + tmp7_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=-C_0*w19 - 2*C_1*w20;
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp3_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp3_1 + tmp5_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp6_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp4_1 + tmp6_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp2_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp5_1;  
2766                                  }                                  }
2767                              }                              }
2768                          }                          }
# Line 3288  void Rectangle::assemblePDESystem(Paso_S Line 2776  void Rectangle::assemblePDESystem(Paso_S
2776                          if (D.actsExpanded()) {                          if (D.actsExpanded()) {
2777                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2778                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2779                                      const double D_0 = D_p[INDEX3(k, m, 0, numEq, numComp)];                                      const double D_0 = D_p[INDEX3(k,m,0,numEq,numComp)];
2780                                      const double D_1 = D_p[INDEX3(k, m, 1, numEq, numComp)];                                      const double D_1 = D_p[INDEX3(k,m,1,numEq,numComp)];
2781                                      const double D_2 = D_p[INDEX3(k, m, 2, numEq, numComp)];                                      const double D_2 = D_p[INDEX3(k,m,2,numEq,numComp)];
2782                                      const double D_3 = D_p[INDEX3(k, m, 3, numEq, numComp)];                                      const double D_3 = D_p[INDEX3(k,m,3,numEq,numComp)];
2783                                      const double tmp4_0 = D_1 + D_3;                                      const double tmp0 = w21*(D_0 + D_1);
2784                                      const double tmp2_0 = D_0 + D_1 + D_2 + D_3;                                      const double tmp1 = w22*(D_2 + D_3);
2785                                      const double tmp5_0 = D_0 + D_2;                                      const double tmp2 = w23*(D_0 + D_1 + D_2 + D_3);
2786                                      const double tmp0_0 = D_0 + D_1;                                      const double tmp3 = w21*(D_2 + D_3);
2787                                      const double tmp6_0 = D_0 + D_3;                                      const double tmp4 = w22*(D_0 + D_1);
2788                                      const double tmp1_0 = D_2 + D_3;                                      const double tmp5 = w23*(D_1 + D_2);
2789                                      const double tmp3_0 = D_1 + D_2;                                      const double tmp6 = w21*(D_1 + D_3);
2790                                      const double tmp16_1 = D_1*w56;                                      const double tmp7 = w22*(D_0 + D_2);
2791                                      const double tmp14_1 = tmp6_0*w54;                                      const double tmp8 = w21*(D_0 + D_2);
2792                                      const double tmp8_1 = D_3*w55;                                      const double tmp9 = w22*(D_1 + D_3);
2793                                      const double tmp2_1 = tmp2_0*w54;                                      const double tmp10 = w23*(D_0 + D_3);
2794                                      const double tmp12_1 = tmp5_0*w52;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=D_0*w24 + D_3*w25 + tmp5;
2795                                      const double tmp4_1 = tmp0_0*w53;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0 + tmp1;
2796                                      const double tmp3_1 = tmp1_0*w52;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp8 + tmp9;
2797                                      const double tmp13_1 = tmp4_0*w53;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp2;
2798                                      const double tmp10_1 = tmp4_0*w52;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0 + tmp1;
2799                                      const double tmp1_1 = tmp1_0*w53;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=D_1*w24 + D_2*w25 + tmp10;
2800                                      const double tmp7_1 = D_3*w56;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2;
2801                                      const double tmp0_1 = tmp0_0*w52;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp6 + tmp7;
2802                                      const double tmp11_1 = tmp5_0*w53;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp8 + tmp9;
2803                                      const double tmp9_1 = D_0*w56;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp2;
2804                                      const double tmp5_1 = tmp3_0*w54;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=D_1*w25 + D_2*w24 + tmp10;
2805                                      const double tmp18_1 = D_2*w56;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp3 + tmp4;
2806                                      const double tmp17_1 = D_1*w55;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2;
2807                                      const double tmp6_1 = D_0*w55;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp6 + tmp7;
2808                                      const double tmp15_1 = D_2*w55;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3 + tmp4;
2809                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=D_0*w25 + D_3*w24 + tmp5;
                                     EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp7_1;  
                                     EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp5_1 + tmp8_1 + tmp9_1;  
                                     EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp10_1 + tmp11_1;  
                                     EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp12_1 + tmp13_1;  
                                     EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp12_1 + tmp13_1;  
                                     EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp10_1 + tmp11_1;  
                                     EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp3_1 + tmp4_1;  
                                     EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp14_1 + tmp15_1 + tmp16_1;  
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp2_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp14_1 + tmp17_1 + tmp18_1;  
2810                                  }                                  }
2811                               }                               }
2812                          } else { // constant data                          } else { // constant data
2813                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2814                                  for (index_t m=0; m<numComp; m++) {                                  for (index_t m=0; m<numComp; m++) {
2815                                      const double D_0 = D_p[INDEX2(k, m, numEq)];                                      const double D_0 = D_p[INDEX2(k, m, numEq)];
2816                                      const double tmp2_1 = D_0*w59;                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=16*D_0*w23;
2817                                      const double tmp1_1 = D_0*w58;                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=8*D_0*w23;
2818                                      const double tmp0_1 = D_0*w57;                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=8*D_0*w23;
2819                                      EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=4*D_0*w23;
2820                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=8*D_0*w23;
2821                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=16*D_0*w23;
2822                                      EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp2_1;                                      EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=4*D_0*w23;
2823                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp2_1;                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=8*D_0*w23;
2824                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=8*D_0*w23;
2825                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=4*D_0*w23;
2826                                      EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1;                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=16*D_0*w23;
2827                                      EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=8*D_0*w23;
2828                                      EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=4*D_0*w23;
2829                                      EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=8*D_0*w23;
2830                                      EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;                                      EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=8*D_0*w23;
2831                                      EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp2_1;                                      EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=16*D_0*w23;
                                     EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;  
                                     EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1;  
                                     EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp2_1;  
2832                                  }                                  }
2833                              }                              }
2834                          }                          }
# Line 3371  void Rectangle::assemblePDESystem(Paso_S Line 2841  void Rectangle::assemblePDESystem(Paso_S
2841                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);                          const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2842                          if (X.actsExpanded()) {                          if (X.actsExpanded()) {
2843                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2844                                  const double X_0_0 = X_p[INDEX3(k, 0, 0, numEq, 2)];                                  const double X_0_0 = X_p[INDEX3(k,0,0,numEq,2)];
2845                                  const double X_1_0 = X_p[INDEX3(k, 1, 0, numEq, 2)];                                  const double X_1_0 = X_p[INDEX3(k,1,0,numEq,2)];
2846                                  const double X_0_1 = X_p[INDEX3(k, 0, 1, numEq, 2)];                                  const double X_0_1 = X_p[INDEX3(k,0,1,numEq,2)];
2847                                  const double X_1_1 = X_p[INDEX3(k, 1, 1, numEq, 2)];                                  const double X_1_1 = X_p[INDEX3(k,1,1,numEq,2)];
2848                                  const double X_0_2 = X_p[INDEX3(k, 0, 2, numEq, 2)];                                  const double X_0_2 = X_p[INDEX3(k,0,2,numEq,2)];
2849                                  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)];
2850                                  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)];
2851                                  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)];
2852                                  const double tmp1_0 = X_1_1 + X_1_3;                                  const double tmp0 = 6*w15*(X_1_1 + X_1_3);
2853                                  const double tmp3_0 = X_0_0 + X_0_1;                                  const double tmp1 = 6*w16*(X_0_2 + X_0_3);
2854                                  const double tmp2_0 = X_1_0 + X_1_2;                                  const double tmp2 = 6*w11*(X_0_0 + X_0_1);
2855                                  const double tmp0_0 = X_0_2 + X_0_3;                                  const double tmp3 = 6*w12*(X_1_0 + X_1_2);
2856                                  const double tmp8_1 = tmp2_0*w66;                                  const double tmp4 = 6*w15*(X_1_0 + X_1_2);
2857                                  const double tmp5_1 = tmp3_0*w64;                                  const double tmp5 = w27*(X_0_2 + X_0_3);
2858                                  const double tmp14_1 = tmp0_0*w64;                                  const double tmp6 = w26*(X_0_0 + X_0_1);
2859                                  const double tmp3_1 = tmp3_0*w60;                                  const double tmp7 = 6*w12*(X_1_1 + X_1_3);
2860                                  const double tmp9_1 = tmp3_0*w63;                                  const double tmp8 = w29*(X_1_1 + X_1_3);
2861                                  const double tmp13_1 = tmp3_0*w65;                                  const double tmp9 = w27*(-X_0_0 - X_0_1);
2862                                  const double tmp12_1 = tmp1_0*w66;                                  const double tmp10 = w28*(X_1_0 + X_1_2);
2863                                  const double tmp10_1 = tmp0_0*w60;                                  const double tmp11 = w26*(-X_0_2 - X_0_3);
2864                                  const double tmp2_1 = tmp2_0*w61;                                  const double tmp12 = w29*(X_1_0 + X_1_2);
2865                                  const double tmp6_1 = tmp2_0*w62;                                  const double tmp13 = w27*(X_0_0 + X_0_1);
2866                                  const double tmp4_1 = tmp0_0*w65;                                  const double tmp14 = w28*(X_1_1 + X_1_3);
2867                                  const double tmp11_1 = tmp1_0*w67;                                  const double tmp15 = w26*(X_0_2 + X_0_3);
2868                                  const double tmp1_1 = tmp1_0*w62;                                  EM_F[INDEX2(k,0,numEq)]+=tmp0 + tmp1 + tmp2 + tmp3;
2869                                  const double tmp7_1 = tmp1_0*w61;                                  EM_F[INDEX2(k,1,numEq)]+=tmp4 + tmp5 + tmp6 + tmp7;
2870                                  const double tmp0_1 = tmp0_0*w63;                                  EM_F[INDEX2(k,2,numEq)]+=tmp10 + tmp11 + tmp8 + tmp9;
2871                                  const double tmp15_1 = tmp2_0*w67;                                  EM_F[INDEX2(k,3,numEq)]+=tmp12 + tmp13 + tmp14 + tmp15;
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp10_1 + tmp11_1 + tmp8_1 + tmp9_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1;  
2872                              }                              }
2873                          } else { // constant data                          } else { // constant data
2874                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2875                                  const double X_0 = X_p[INDEX2(k, 0, numEq)];                                  const double X_0 = X_p[INDEX2(k, 0, numEq)];
2876                                  const double X_1 = X_p[INDEX2(k, 1, numEq)];                                  const double X_1 = X_p[INDEX2(k, 1, numEq)];
2877                                  const double tmp0_1 = X_1*w69;                                  EM_F[INDEX2(k,0,numEq)]+=3*X_0*w19 + 6*X_1*w20;
2878                                  const double tmp1_1 = X_0*w68;                                  EM_F[INDEX2(k,1,numEq)]+=-3*X_0*w19 + 6*X_1*w20;
2879                                  const double tmp2_1 = X_0*w70;                                  EM_F[INDEX2(k,2,numEq)]+=3*X_0*w19 - 6*X_1*w20;
2880                                  const double tmp3_1 = X_1*w71;                                  EM_F[INDEX2(k,3,numEq)]+=-3*X_0*w19 - 6*X_1*w20;
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp0_1 + tmp2_1;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp1_1 + tmp3_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;  
2881                              }                              }
2882                          }                          }
2883                      }                      }
# Line 3431  void Rectangle::assemblePDESystem(Paso_S Line 2893  void Rectangle::assemblePDESystem(Paso_S
2893                                  const double Y_1 = Y_p[INDEX2(k, 1, numEq)];                                  const double Y_1 = Y_p[INDEX2(k, 1, numEq)];
2894                                  const double Y_2 = Y_p[INDEX2(k, 2, numEq)];                                  const double Y_2 = Y_p[INDEX2(k, 2, numEq)];
2895                                  const double Y_3 = Y_p[INDEX2(k, 3, numEq)];                                  const double Y_3 = Y_p[INDEX2(k, 3, numEq)];
2896                                  const double tmp0_0 = Y_1 + Y_2;                                  const double tmp0 = 6*w23*(Y_1 + Y_2);
2897                                  const double tmp1_0 = Y_0 + Y_3;                                  const double tmp1 = 6*w23*(Y_0 + Y_3);
2898                                  const double tmp0_1 = Y_0*w72;                                  EM_F[INDEX2(k,0,numEq)]+=6*Y_0*w21 + 6*Y_3*w22 + tmp0;
2899                                  const double tmp1_1 = tmp0_0*w73;                                  EM_F[INDEX2(k,1,numEq)]+=6*Y_1*w21 + 6*Y_2*w22 + tmp1;
2900                                  const double tmp2_1 = Y_3*w74;                                  EM_F[INDEX2(k,2,numEq)]+=6*Y_1*w22 + 6*Y_2*w21 + tmp1;
2901                                  const double tmp3_1 = Y_1*w72;                                  EM_F[INDEX2(k,3,numEq)]+=6*Y_0*w22 + 6*Y_3*w21 + tmp0;
                                 const double tmp4_1 = tmp1_0*w73;  
                                 const double tmp5_1 = Y_2*w74;  
                                 const double tmp6_1 = Y_2*w72;  
                                 const double tmp7_1 = Y_1*w74;  
                                 const double tmp8_1 = Y_3*w72;  
                                 const double tmp9_1 = Y_0*w74;  
                                 EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1 + tmp2_1;  
                                 EM_F[INDEX2(k,1,numEq)]+=tmp3_1 + tmp4_1 + tmp5_1;  
                                 EM_F[INDEX2(k,2,numEq)]+=tmp4_1 + tmp6_1 + tmp7_1;  
                                 EM_F[INDEX2(k,3,numEq)]+=tmp1_1 + tmp8_1 + tmp9_1;  
2902                              }                              }
2903                          } else { // constant data                          } else { // constant data
2904                              for (index_t k=0; k<numEq; k++) {                              for (index_t k=0; k<numEq; k++) {
2905                                  const double tmp0_1 = Y_p[k]*w75;                                  const double Y_0 = Y_p[k];
2906                                  EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,0,numEq)]+=36*Y_0*w23;
2907                                  EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,1,numEq)]+=36*Y_0*w23;
2908                                  EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,2,numEq)]+=36*Y_0*w23;
2909                                  EM_F[INDEX2(k,3,numEq)]+=tmp0_1;                                  EM_F[INDEX2(k,3,numEq)]+=36*Y_0*w23;
2910                              }                              }
2911                          }                          }
2912                      }                      }
2913                        /* GENERATOR SNIP_PDE_SYSTEM BOTTOM */
2914    
2915                      // 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)
2916                      const index_t firstNode=m_NN[0]*k1+k0;                      const index_t firstNode=m_NN[0]*k1+k0;
# Line 3483  void Rectangle::assemblePDESystemReduced Line 2936  void Rectangle::assemblePDESystemReduced
2936          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
2937      }      }
2938    
2939      const double w0 = -.25*m_dx[1]/m_dx[0];      const double w0 = .25;
2940      const double w1 = .25;      const double w1 = .125*m_dx[0];
2941      const double w2 = -.25;      const double w2 = .125*m_dx[1];
2942      const double w3 = .25*m_dx[0]/m_dx[1];      const double w3 = .0625*m_dx[0]*m_dx[1];
2943      const double w4 = -.25*m_dx[0]/m_dx[1];      const double w4 = .25*m_dx[0]/m_dx[1];
2944      const double w5 = .25*m_dx[1]/m_dx[0];      const double w5 = .25*m_dx[1]/m_dx[0];
     const double w6 = -.125*m_dx[1];  
     const double w7 = -.125*m_dx[0];  
     const double w8 = .125*m_dx[1];  
     const double w9 = .125*m_dx[0];  
     const double w10 = .0625*m_dx[0]*m_dx[1];  
     const double w11 = -.5*m_dx[1];  
     const double w12 = -.5*m_dx[0];  
     const double w13 = .5*m_dx[1];  
     const double w14 = .5*m_dx[0];  
     const double w15 = .25*m_dx[0]*m_dx[1];  
2945    
2946      rhs.requireWrite();      rhs.requireWrite();
2947  #pragma omp parallel  #pragma omp parallel
# Line 3524  void Rectangle::assemblePDESystemReduced Line 2967  void Rectangle::assemblePDESystemReduced
2967                                  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)];
2968                                  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)];
2969                                  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)];
2970                                  const double tmp0_0 = A_01 + A_10;                                  const double tmp0 = (A_01 + A_10)*w0;
2971                                  const double tmp0_1 = A_11*w3;                                  const double tmp1 = A_00*w5;
2972                                  const double tmp1_1 = A_00*w0;                                  const double tmp2 = A_01*w0;
2973                                  const double tmp2_1 = A_01*w1;                                  const double tmp3 = A_10*w0;
2974                                  const double tmp3_1 = A_10*w2;                                  const double tmp4 = A_11*w4;
2975                                  const double tmp4_1 = tmp0_0*w1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+= tmp4 + tmp0 + tmp1;
2976                                  const double tmp5_1 = A_11*w4;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+= tmp4 - tmp1 + tmp3 - tmp2;
2977                                  const double tmp6_1 = A_00*w5;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+= tmp2 - tmp3 - tmp4 + tmp1;
2978                                  const double tmp7_1 = tmp0_0*w2;                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=-tmp1 - tmp4 - tmp0;
2979                                  const double tmp8_1 = A_10*w1;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+= tmp4 - tmp1 + tmp2 - tmp3;
2980                                  const double tmp9_1 = A_01*w2;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+= tmp4 + tmp1 - tmp0;
2981                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=-tmp1 + tmp0 - tmp4;
2982                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;
2983                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-tmp4 + tmp1 + tmp3 - tmp2;
2984                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=-tmp1 + tmp0 - tmp4;
2985                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+= tmp4 + tmp1 - tmp0;
2986                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+= tmp4 - tmp1 + tmp2 - tmp3;
2987                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=-tmp1 - tmp4 - tmp0;
2988                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+= tmp2 - tmp3 - tmp4 + tmp1;
2989                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+= tmp4 - tmp1 + tmp3 - tmp2;
2990                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+= tmp4 + tmp0 + tmp1;
                                 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;  
                                 EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;  
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;  
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;  
2991                              }                              }
2992                          }                          }
2993                      }                      }
# Line 3562  void Rectangle::assemblePDESystemReduced Line 2999  void Rectangle::assemblePDESystemReduced
2999                          const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);                          const double* B_p=const_cast<escript::Data*>(&B)->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 B_0 = B_p[INDEX3(k,0,m, numEq, 2)];                                  const double tmp0 = B_p[INDEX3(k,0,m, numEq, 2)]*w2;
3003                                  const double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];                                  const double tmp1 = B_p[INDEX3(k,1,m, numEq, 2)]*w1;
3004                                  const double tmp0_1 = B_0*w6;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-tmp0 - tmp1;
3005                                  const double tmp1_1 = B_1*w7;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-tmp1 + tmp0;
3006                                  const double tmp2_1 = B_0*w8;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-tmp0 + tmp1;
3007                                  const double tmp3_1 = B_1*w9;                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+= tmp0 + tmp1;
3008                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=-tmp0 - tmp1;
3009                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=-tmp1 + tmp0;
3010                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=-tmp0 + tmp1;
3011                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+= tmp0 + tmp1;
3012                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=-tmp0 - tmp1;
3013                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp1_1 + tmp2_1;                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=-tmp1 + tmp0;
3014                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=-tmp0 + tmp1;
3015                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+= tmp0 + tmp1;
3016                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=-tmp0 - tmp1;
3017                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=-tmp1 + tmp0;
3018                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=-tmp0 + tmp1;
3019                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+= tmp0 + tmp1;
                                 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;  
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp3_1;  
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1;  
3020                              }                              }
3021                          }                          }
3022                      }                      }
# Line 3595  void Rectangle::assemblePDESystemReduced Line 3028  void Rectangle::assemblePDESystemReduced
3028                          const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);                          const double* C_p=const_cast<escript::Data*>(&C)->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 C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];                                  const double tmp0 = C_p[INDEX3(k, m, 0, numEq, numComp)]*w2;
3032                                  const double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];                                  const double tmp1 = C_p[INDEX3(k, m, 1, numEq, numComp)]*w1;
3033                                  const double tmp0_1 = C_0*w8;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=-tmp1 - tmp0;
3034                                  const double tmp1_1 = C_1*w7;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=-tmp1 - tmp0;
3035                                  const double tmp2_1 = C_1*w9;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=-tmp1 - tmp0;
3036                                  const double tmp3_1 = C_0*w6;                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=-tmp1 - tmp0;
3037                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+= tmp0 - tmp1;
3038                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+= tmp0 - tmp1;
3039                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+= tmp0 - tmp1;
3040                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+= tmp0 - tmp1;
3041                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+= tmp1 - tmp0;
3042                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+= tmp1 - tmp0;
3043                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+= tmp1 - tmp0;
3044                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+= tmp1 - tmp0;
3045                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+= tmp0 + tmp1;
3046                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+= tmp0 + tmp1;
3047                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+= tmp0 + tmp1;
3048                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+= tmp0 + tmp1;
                                 EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;  
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;  
3049                              }                              }
3050                          }                          }
3051                      }                      }
# Line 3628  void Rectangle::assemblePDESystemReduced Line 3057  void Rectangle::assemblePDESystemReduced
3057                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);                          const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
3058                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3059                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3060                                  const double tmp0_1 = D_p[INDEX2(k, m, numEq)]*w10;                                  const double tmp0 = D_p[INDEX2(k, m, numEq)]*w3;
3061                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0;
3062                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0;
3063                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0;
3064                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0;
3065                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0;
3066                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0;
3067                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0;
3068                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0;
3069                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0;
3070                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0;
3071                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0;
3072                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0;
3073                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0;
3074                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0;
3075                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0;
3076                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0;
3077                              }                              }
3078                          }                          }
3079                      }                      }
# Line 3657  void Rectangle::assemblePDESystemReduced Line 3086  void Rectangle::assemblePDESystemReduced
3086                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3087                              const double X_0 = X_p[INDEX2(k, 0, numEq)];                              const double X_0 = X_p[INDEX2(k, 0, numEq)];
3088                              const double X_1 = X_p[INDEX2(k, 1, numEq)];                              const double X_1 = X_p[INDEX2(k, 1, numEq)];
3089                              const double tmp0_1 = X_0*w11;                              const double tmp0 = 4*X_0*w2;
3090                              const double tmp1_1 = X_1*w12;                              const double tmp1 = 4*X_1*w1;
3091                              const double tmp2_1 = X_0*w13;                              EM_F[INDEX2(k,0,numEq)]+=-tmp0 - tmp1;
3092                              const double tmp3_1 = X_1*w14;                              EM_F[INDEX2(k,1,numEq)]+= tmp0 - tmp1;
3093                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;                              EM_F[INDEX2(k,2,numEq)]+= tmp1 - tmp0;
3094                              EM_F[INDEX2(k,1,numEq)]+=tmp1_1 + tmp2_1;                              EM_F[INDEX2(k,3,numEq)]+= tmp0 + tmp1;
                             EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp3_1;  
                             EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;  
3095                          }                          }
3096                      }                      }
3097                      ///////////////                      ///////////////
# Line 3674  void Rectangle::assemblePDESystemReduced Line 3101  void Rectangle::assemblePDESystemReduced
3101                          add_EM_F=true;                          add_EM_F=true;
3102                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);                          const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
3103                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3104                              const double tmp0_1 = Y_p[k]*w15;                              EM_F[INDEX2(k,0,numEq)]+=4*Y_p[k]*w3;
3105                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,1,numEq)]+=4*Y_p[k]*w3;
3106                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,2,numEq)]+=4*Y_p[k]*w3;
3107                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,3,numEq)]+=4*Y_p[k]*w3;
                             EM_F[INDEX2(k,3,numEq)]+=tmp0_1;  
3108                          }                          }
3109                      }                      }
3110    
# Line 3971  void Rectangle::assemblePDEBoundarySingl Line 3397  void Rectangle::assemblePDEBoundarySingl
3397  void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,  void Rectangle::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
3398        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
3399  {  {
3400      const double w0 = 0.25*m_dx[1];      const double w0 = 0.25*m_dx[0];
3401      const double w1 = 0.5*m_dx[1];      const double w1 = 0.25*m_dx[1];
     const double w2 = 0.25*m_dx[0];  
     const double w3 = 0.5*m_dx[0];  
3402      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3403      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3404      rhs.requireWrite();      rhs.requireWrite();
# Line 3986  void Rectangle::assemblePDEBoundarySingl Line 3410  void Rectangle::assemblePDEBoundarySingl
3410                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
3411                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3412                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
                     const index_t e = k1;  
3413                      ///////////////                      ///////////////
3414                      // process d //                      // process d //
3415                      ///////////////                      ///////////////
3416                      if (add_EM_S) {                      if (add_EM_S) {
3417                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(k1);
3418                          const double d_0 = d_p[0];                          const double tmp0 = d_p[0]*w1;
3419                          const double tmp0_1 = d_0*w0;                          EM_S[INDEX2(0,0,4)]+=tmp0;
3420                          EM_S[INDEX2(0,0,4)]+=tmp0_1;                          EM_S[INDEX2(2,0,4)]+=tmp0;
3421                          EM_S[INDEX2(2,0,4)]+=tmp0_1;                          EM_S[INDEX2(0,2,4)]+=tmp0;
3422                          EM_S[INDEX2(0,2,4)]+=tmp0_1;                          EM_S[INDEX2(2,2,4)]+=tmp0;
                         EM_S[INDEX2(2,2,4)]+=tmp0_1;  
3423                      }                      }
3424                      ///////////////                      ///////////////
3425                      // process y //                      // process y //
3426                      ///////////////                      ///////////////
3427                      if (add_EM_F) {                      if (add_EM_F) {
3428                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(k1);
3429                          const double tmp0_1 = w1*y_p[0];                          const double tmp0 = 2*w1*y_p[0];
3430                          EM_F[0]+=tmp0_1;                          EM_F[0]+=tmp0;
3431                          EM_F[2]+=tmp0_1;                          EM_F[2]+=tmp0;
3432                      }                      }
3433                      const index_t firstNode=m_NN[0]*k1;                      const index_t firstNode=m_NN[0]*k1;
3434                      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 4026  void Rectangle::assemblePDEBoundarySingl Line 3448  void Rectangle::assemblePDEBoundarySingl
3448                      ///////////////                      ///////////////
3449                      if (add_EM_S) {                      if (add_EM_S) {
3450                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3451                          const double d_0 = d_p[0];                          const double tmp0 = d_p[0]*w1;
3452                          const double tmp0_1 = d_0*w0;                          EM_S[INDEX2(1,1,4)]+=tmp0;
3453                          EM_S[INDEX2(1,1,4)]+=tmp0_1;                          EM_S[INDEX2(3,1,4)]+=tmp0;
3454                          EM_S[INDEX2(3,1,4)]+=tmp0_1;                          EM_S[INDEX2(1,3,4)]+=tmp0;
3455                          EM_S[INDEX2(1,3,4)]+=tmp0_1;                          EM_S[INDEX2(3,3,4)]+=tmp0;
                         EM_S[INDEX2(3,3,4)]+=tmp0_1;  
3456                      }                      }
3457                      ///////////////                      ///////////////
3458                      // process y //                      // process y //
3459                      ///////////////                      ///////////////
3460                      if (add_EM_F) {                      if (add_EM_F) {
3461                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3462                          const double tmp0_1 = w1*y_p[0];                          const double tmp0 = 2*w1*y_p[0];
3463                          EM_F[1]+=tmp0_1;                          EM_F[1]+=tmp0;
3464                          EM_F[3]+=tmp0_1;                          EM_F[3]+=tmp0;
3465                      }                      }
3466                      const index_t firstNode=m_NN[0]*(k1+1)-2;                      const index_t firstNode=m_NN[0]*(k1+1)-2;
3467                      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 4060  void Rectangle::assemblePDEBoundarySingl Line 3481  void Rectangle::assemblePDEBoundarySingl
3481                      ///////////////                      ///////////////
3482                      if (add_EM_S) {                      if (add_EM_S) {
3483                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3484                          const double tmp0_1 = d_p[0]*w2;                          const double tmp0 = d_p[0]*w0;
3485                          EM_S[INDEX2(0,0,4)]+=tmp0_1;                          EM_S[INDEX2(0,0,4)]+=tmp0;
3486                          EM_S[INDEX2(1,0,4)]+=tmp0_1;                          EM_S[INDEX2(1,0,4)]+=tmp0;
3487                          EM_S[INDEX2(0,1,4)]+=tmp0_1;                          EM_S[INDEX2(0,1,4)]+=tmp0;
3488                          EM_S[INDEX2(1,1,4)]+=tmp0_1;                          EM_S[INDEX2(1,1,4)]+=tmp0;
3489                      }                      }
3490                      ///////////////                      ///////////////
3491                      // process y //                      // process y //
3492                      ///////////////                      ///////////////
3493                      if (add_EM_F) {                      if (add_EM_F) {
3494                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3495                          const double tmp0_1 = w3*y_p[0];                          const double tmp0 = 2*w0*y_p[0];
3496                          EM_F[0]+=tmp0_1;                          EM_F[0]+=tmp0;
3497                          EM_F[1]+=tmp0_1;                          EM_F[1]+=tmp0;
3498                      }                      }
3499                      const index_t firstNode=k0;                      const index_t firstNode=k0;
3500                      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 4093  void Rectangle::assemblePDEBoundarySingl Line 3514  void Rectangle::assemblePDEBoundarySingl
3514                      ///////////////                      ///////////////
3515                      if (add_EM_S) {                      if (add_EM_S) {
3516                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3517                          const double tmp0_1 = d_p[0]*w2;                          const double tmp0 = d_p[0]*w0;
3518                          EM_S[INDEX2(2,2,4)]+=tmp0_1;                          EM_S[INDEX2(2,2,4)]+=tmp0;
3519                          EM_S[INDEX2(3,2,4)]+=tmp0_1;                          EM_S[INDEX2(3,2,4)]+=tmp0;
3520                          EM_S[INDEX2(2,3,4)]+=tmp0_1;                          EM_S[INDEX2(2,3,4)]+=tmp0;
3521                          EM_S[INDEX2(3,3,4)]+=tmp0_1;                          EM_S[INDEX2(3,3,4)]+=tmp0;
3522                      }                      }
3523                      ///////////////                      ///////////////
3524                      // process y //                      // process y //
3525                      ///////////////                      ///////////////
3526                      if (add_EM_F) {                      if (add_EM_F) {
3527                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3528                          const double tmp0_1 = w3*y_p[0];                          const double tmp0 = 2*w0*y_p[0];
3529                          EM_F[2]+=tmp0_1;                          EM_F[2]+=tmp0;
3530                          EM_F[3]+=tmp0_1;                          EM_F[3]+=tmp0;
3531                      }                      }
3532                      const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;                      const index_t firstNode=m_NN[0]*(m_NN[1]-2)+k0;
3533                      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 4461  void Rectangle::assemblePDEBoundarySyste Line 3882  void Rectangle::assemblePDEBoundarySyste
3882          numEq=mat->logical_row_block_size;          numEq=mat->logical_row_block_size;
3883          numComp=mat->logical_col_block_size;          numComp=mat->logical_col_block_size;
3884      }      }
3885      const double w0 = 0.25*m_dx[1];      const double w0 = 0.25*m_dx[0];
3886      const double w1 = 0.5*m_dx[1];      const double w1 = 0.25*m_dx[1];
     const double w2 = 0.25*m_dx[0];  
     const double w3 = 0.5*m_dx[0];  
3887      const bool add_EM_S=!d.isEmpty();      const bool add_EM_S=!d.isEmpty();
3888      const bool add_EM_F=!y.isEmpty();      const bool add_EM_F=!y.isEmpty();
3889      rhs.requireWrite();      rhs.requireWrite();
# Line 4476  void Rectangle::assemblePDEBoundarySyste Line 3895  void Rectangle::assemblePDEBoundarySyste
3895                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE[1]; k1+=2) {
3896                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
3897                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
                     const index_t e = k1;  
3898                      ///////////////                      ///////////////
3899                      // process d //                      // process d //
3900                      ///////////////                      ///////////////
3901                      if (add_EM_S) {                      if (add_EM_S) {
3902                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(k1);
3903                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3904                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3905                                  const double d_0 = d_p[INDEX2(k, m, numEq)];                                  const double tmp0 = d_p[INDEX2(k, m, numEq)]*w1;
3906                                  const double tmp0_1 = d_0*w0;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0;
3907                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0;
3908                                  EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0;
3909                                  EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0;
                                 EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;  
3910                              }                              }
3911                          }                          }
3912                      }                      }
# Line 4497  void Rectangle::assemblePDEBoundarySyste Line 3914  void Rectangle::assemblePDEBoundarySyste
3914                      // process y //                      // process y //
3915                      ///////////////                      ///////////////
3916                      if (add_EM_F) {                      if (add_EM_F) {
3917                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(k1);
3918                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3919                              const double tmp0_1 = w1*y_p[k];                              const double tmp0 = 2*w1*y_p[k];
3920                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,0,numEq)]+=tmp0;
3921                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,2,numEq)]+=tmp0;
3922                          }                          }
3923                      }                      }
3924                      const index_t firstNode=m_NN[0]*k1;                      const index_t firstNode=m_NN[0]*k1;
# Line 4525  void Rectangle::assemblePDEBoundarySyste Line 3942  void Rectangle::assemblePDEBoundarySyste
3942                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3943                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3944                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3945                                  const double d_0 = d_p[INDEX2(k, m, numEq)];                                  const double tmp0 = d_p[INDEX2(k, m, numEq)]*w1;
3946                                  const double tmp0_1 = d_0*w0;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0;
3947                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0;
3948                                  EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0;
3949                                  EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0;
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;  
3950                              }                              }
3951                          }                          }
3952                      }                      }
# Line 4540  void Rectangle::assemblePDEBoundarySyste Line 3956  void Rectangle::assemblePDEBoundarySyste
3956                      if (add_EM_F) {                      if (add_EM_F) {
3957                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3958                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3959                              const double tmp0_1 = w1*y_p[k];                              const double tmp0 = 2*w1*y_p[k];
3960                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,1,numEq)]+=tmp0;
3961                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,3,numEq)]+=tmp0;
3962                          }                          }
3963                      }                      }
3964                      const index_t firstNode=m_NN[0]*(k1+1)-2;                      const index_t firstNode=m_NN[0]*(k1+1)-2;
# Line 4566  void Rectangle::assemblePDEBoundarySyste Line 3982  void Rectangle::assemblePDEBoundarySyste
3982                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
3983                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3984                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
3985                                  const double d_0 = d_p[INDEX2(k, m, numEq)];                                  const double tmp0 = d_p[INDEX2(k, m, numEq)]*w0;
3986                                  const double tmp0_1 = d_0*w2;                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0;
3987                                  EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0;
3988                                  EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0;
3989                                  EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0;
                                 EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;  
3990                              }                              }
3991                          }                          }
3992                      }                      }
# Line 4581  void Rectangle::assemblePDEBoundarySyste Line 3996  void Rectangle::assemblePDEBoundarySyste
3996                      if (add_EM_F) {                      if (add_EM_F) {
3997                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
3998                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
3999                              const double tmp0_1 = w3*y_p[k];                              const double tmp0 = 2*w0*y_p[k];
4000                              EM_F[INDEX2(k,0,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,0,numEq)]+=tmp0;
4001                              EM_F[INDEX2(k,1,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,1,numEq)]+=tmp0;
4002                          }                          }
4003                      }                      }
4004                      const index_t firstNode=k0;                      const index_t firstNode=k0;
# Line 4607  void Rectangle::assemblePDEBoundarySyste Line 4022  void Rectangle::assemblePDEBoundarySyste
4022                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);                          const double* d_p=const_cast<escript::Data*>(&d)->getSampleDataRO(e);
4023                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
4024                              for (index_t m=0; m<numComp; m++) {                              for (index_t m=0; m<numComp; m++) {
4025                                  const double d_0 = d_p[INDEX2(k, m, numEq)];                                  const double tmp0 = d_p[INDEX2(k, m, numEq)]*w0;
4026                                  const double tmp0_1 = d_0*w2;                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0;
4027                                  EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0;
4028                                  EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0;
4029                                  EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;                                  EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0;
                                 EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;  
4030                              }                              }
4031                          }                          }
4032                      }                      }
# Line 4622  void Rectangle::assemblePDEBoundarySyste Line 4036  void Rectangle::assemblePDEBoundarySyste
4036                      if (add_EM_F) {                      if (add_EM_F) {
4037                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);                          const double* y_p=const_cast<escript::Data*>(&y)->getSampleDataRO(e);
4038                          for (index_t k=0; k<numEq; k++) {                          for (index_t k=0; k<numEq; k++) {
4039                              const double tmp0_1 = w3*y_p[k];                              const double tmp0 = 2*w0*y_p[k];
4040                              EM_F[INDEX2(k,2,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,2,numEq)]+=tmp0;
4041                              EM_F[INDEX2(k,3,numEq)]+=tmp0_1;                              EM_F[INDEX2(k,3,numEq)]+=tmp0;
4042                          }                          }
4043                      }                      }
4044                      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.4368  
changed lines
  Added in v.4370

  ViewVC Help
Powered by ViewVC 1.1.26