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

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

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

revision 4644 by sshaw, Fri Jan 24 03:29:25 2014 UTC revision 4645 by sshaw, Mon Feb 3 00:06:24 2014 UTC
# Line 1978  void WaveAssembler3D::assemblePDESystem( Line 1978  void WaveAssembler3D::assemblePDESystem(
1978                                           *c33_p = c33.getSampleDataRO(e),                                           *c33_p = c33.getSampleDataRO(e),
1979                                           *c44_p = c44.getSampleDataRO(e),                                           *c44_p = c44.getSampleDataRO(e),
1980                                           *c66_p = c66.getSampleDataRO(e);                                           *c66_p = c66.getSampleDataRO(e);
 #define NUM_QUADPOINTS 8  
                             double *X_00 = new double[NUM_QUADPOINTS];  
                             double *X_01 = new double[NUM_QUADPOINTS];  
                             double *X_02 = new double[NUM_QUADPOINTS];  
                             double *X_11 = new double[NUM_QUADPOINTS];  
                             double *X_12 = new double[NUM_QUADPOINTS];  
                             double *X_22 = new double[NUM_QUADPOINTS];  
1981                              if (du.actsExpanded()) {                              if (du.actsExpanded()) {
1982    const double X_00_0 = -(c11_p[0]*du_p[INDEX3(0,0,0,numEq,3)]+c12_p[0]*du_p[INDEX3(1,1,0,numEq,3)]+c13_p[0]*du_p[INDEX3(2,2,0,numEq,3)]);
1983    const double X_00_1 = -(c11_p[1]*du_p[INDEX3(0,0,1,numEq,3)]+c12_p[1]*du_p[INDEX3(1,1,1,numEq,3)]+c13_p[1]*du_p[INDEX3(2,2,1,numEq,3)]);
1984    const double X_00_2 = -(c11_p[2]*du_p[INDEX3(0,0,2,numEq,3)]+c12_p[2]*du_p[INDEX3(1,1,2,numEq,3)]+c13_p[2]*du_p[INDEX3(2,2,2,numEq,3)]);
1985    const double X_00_3 = -(c11_p[3]*du_p[INDEX3(0,0,3,numEq,3)]+c12_p[3]*du_p[INDEX3(1,1,3,numEq,3)]+c13_p[3]*du_p[INDEX3(2,2,3,numEq,3)]);
1986    const double X_00_4 = -(c11_p[4]*du_p[INDEX3(0,0,4,numEq,3)]+c12_p[4]*du_p[INDEX3(1,1,4,numEq,3)]+c13_p[4]*du_p[INDEX3(2,2,4,numEq,3)]);
1987    const double X_00_5 = -(c11_p[5]*du_p[INDEX3(0,0,5,numEq,3)]+c12_p[5]*du_p[INDEX3(1,1,5,numEq,3)]+c13_p[5]*du_p[INDEX3(2,2,5,numEq,3)]);
1988    const double X_00_6 = -(c11_p[6]*du_p[INDEX3(0,0,6,numEq,3)]+c12_p[6]*du_p[INDEX3(1,1,6,numEq,3)]+c13_p[6]*du_p[INDEX3(2,2,6,numEq,3)]);
1989    const double X_00_7 = -(c11_p[7]*du_p[INDEX3(0,0,7,numEq,3)]+c12_p[7]*du_p[INDEX3(1,1,7,numEq,3)]+c13_p[7]*du_p[INDEX3(2,2,7,numEq,3)]);
1990    const double X_01_0 = -(c66_p[0]*(du_p[INDEX3(0,1,0,numEq,3)]+du_p[INDEX3(1,0,0,numEq,3)]));
1991    const double X_01_1 = -(c66_p[1]*(du_p[INDEX3(0,1,1,numEq,3)]+du_p[INDEX3(1,0,1,numEq,3)]));
1992    const double X_01_2 = -(c66_p[2]*(du_p[INDEX3(0,1,2,numEq,3)]+du_p[INDEX3(1,0,2,numEq,3)]));
1993    const double X_01_3 = -(c66_p[3]*(du_p[INDEX3(0,1,3,numEq,3)]+du_p[INDEX3(1,0,3,numEq,3)]));
1994    const double X_01_4 = -(c66_p[4]*(du_p[INDEX3(0,1,4,numEq,3)]+du_p[INDEX3(1,0,4,numEq,3)]));
1995    const double X_01_5 = -(c66_p[5]*(du_p[INDEX3(0,1,5,numEq,3)]+du_p[INDEX3(1,0,5,numEq,3)]));
1996    const double X_01_6 = -(c66_p[6]*(du_p[INDEX3(0,1,6,numEq,3)]+du_p[INDEX3(1,0,6,numEq,3)]));
1997    const double X_01_7 = -(c66_p[7]*(du_p[INDEX3(0,1,7,numEq,3)]+du_p[INDEX3(1,0,7,numEq,3)]));
1998    const double X_02_0 = -(c44_p[0]*(du_p[INDEX3(2,0,0,numEq,3)]+du_p[INDEX3(0,2,0,numEq,3)]));
1999    const double X_02_1 = -(c44_p[1]*(du_p[INDEX3(2,0,1,numEq,3)]+du_p[INDEX3(0,2,1,numEq,3)]));
2000    const double X_02_2 = -(c44_p[2]*(du_p[INDEX3(2,0,2,numEq,3)]+du_p[INDEX3(0,2,2,numEq,3)]));
2001    const double X_02_3 = -(c44_p[3]*(du_p[INDEX3(2,0,3,numEq,3)]+du_p[INDEX3(0,2,3,numEq,3)]));
2002    const double X_02_4 = -(c44_p[4]*(du_p[INDEX3(2,0,4,numEq,3)]+du_p[INDEX3(0,2,4,numEq,3)]));
2003    const double X_02_5 = -(c44_p[5]*(du_p[INDEX3(2,0,5,numEq,3)]+du_p[INDEX3(0,2,5,numEq,3)]));
2004    const double X_02_6 = -(c44_p[6]*(du_p[INDEX3(2,0,6,numEq,3)]+du_p[INDEX3(0,2,6,numEq,3)]));
2005    const double X_02_7 = -(c44_p[7]*(du_p[INDEX3(2,0,7,numEq,3)]+du_p[INDEX3(0,2,7,numEq,3)]));
2006    const double X_11_0 = -(c12_p[0]*du_p[INDEX3(0,0,0,numEq,3)]+c11_p[0]*du_p[INDEX3(1,1,0,numEq,3)]+c13_p[0]*du_p[INDEX3(2,2,0,numEq,3)]);
2007    const double X_11_1 = -(c12_p[1]*du_p[INDEX3(0,0,1,numEq,3)]+c11_p[1]*du_p[INDEX3(1,1,1,numEq,3)]+c13_p[1]*du_p[INDEX3(2,2,1,numEq,3)]);
2008    const double X_11_2 = -(c12_p[2]*du_p[INDEX3(0,0,2,numEq,3)]+c11_p[2]*du_p[INDEX3(1,1,2,numEq,3)]+c13_p[2]*du_p[INDEX3(2,2,2,numEq,3)]);
2009    const double X_11_3 = -(c12_p[3]*du_p[INDEX3(0,0,3,numEq,3)]+c11_p[3]*du_p[INDEX3(1,1,3,numEq,3)]+c13_p[3]*du_p[INDEX3(2,2,3,numEq,3)]);
2010    const double X_11_4 = -(c12_p[4]*du_p[INDEX3(0,0,4,numEq,3)]+c11_p[4]*du_p[INDEX3(1,1,4,numEq,3)]+c13_p[4]*du_p[INDEX3(2,2,4,numEq,3)]);
2011    const double X_11_5 = -(c12_p[5]*du_p[INDEX3(0,0,5,numEq,3)]+c11_p[5]*du_p[INDEX3(1,1,5,numEq,3)]+c13_p[5]*du_p[INDEX3(2,2,5,numEq,3)]);
2012    const double X_11_6 = -(c12_p[6]*du_p[INDEX3(0,0,6,numEq,3)]+c11_p[6]*du_p[INDEX3(1,1,6,numEq,3)]+c13_p[6]*du_p[INDEX3(2,2,6,numEq,3)]);
2013    const double X_11_7 = -(c12_p[7]*du_p[INDEX3(0,0,7,numEq,3)]+c11_p[7]*du_p[INDEX3(1,1,7,numEq,3)]+c13_p[7]*du_p[INDEX3(2,2,7,numEq,3)]);
2014    const double X_12_0 = -(c44_p[0]*(du_p[INDEX3(2,1,0,numEq,3)]+du_p[INDEX3(1,2,0,numEq,3)]));
2015    const double X_12_1 = -(c44_p[1]*(du_p[INDEX3(2,1,1,numEq,3)]+du_p[INDEX3(1,2,1,numEq,3)]));
2016    const double X_12_2 = -(c44_p[2]*(du_p[INDEX3(2,1,2,numEq,3)]+du_p[INDEX3(1,2,2,numEq,3)]));
2017    const double X_12_3 = -(c44_p[3]*(du_p[INDEX3(2,1,3,numEq,3)]+du_p[INDEX3(1,2,3,numEq,3)]));
2018    const double X_12_4 = -(c44_p[4]*(du_p[INDEX3(2,1,4,numEq,3)]+du_p[INDEX3(1,2,4,numEq,3)]));
2019    const double X_12_5 = -(c44_p[5]*(du_p[INDEX3(2,1,5,numEq,3)]+du_p[INDEX3(1,2,5,numEq,3)]));
2020    const double X_12_6 = -(c44_p[6]*(du_p[INDEX3(2,1,6,numEq,3)]+du_p[INDEX3(1,2,6,numEq,3)]));
2021    const double X_12_7 = -(c44_p[7]*(du_p[INDEX3(2,1,7,numEq,3)]+du_p[INDEX3(1,2,7,numEq,3)]));
2022    const double X_22_0 = -(c13_p[0]*(du_p[INDEX3(0,0,0,numEq,3)] + du_p[INDEX3(1,1,0,numEq,3)]) + c33_p[0]*du_p[INDEX3(2,2,0,numEq,3)]);
2023    const double X_22_1 = -(c13_p[1]*(du_p[INDEX3(0,0,1,numEq,3)] + du_p[INDEX3(1,1,1,numEq,3)]) + c33_p[1]*du_p[INDEX3(2,2,1,numEq,3)]);
2024    const double X_22_2 = -(c13_p[2]*(du_p[INDEX3(0,0,2,numEq,3)] + du_p[INDEX3(1,1,2,numEq,3)]) + c33_p[2]*du_p[INDEX3(2,2,2,numEq,3)]);
2025    const double X_22_3 = -(c13_p[3]*(du_p[INDEX3(0,0,3,numEq,3)] + du_p[INDEX3(1,1,3,numEq,3)]) + c33_p[3]*du_p[INDEX3(2,2,3,numEq,3)]);
2026    const double X_22_4 = -(c13_p[4]*(du_p[INDEX3(0,0,4,numEq,3)] + du_p[INDEX3(1,1,4,numEq,3)]) + c33_p[4]*du_p[INDEX3(2,2,4,numEq,3)]);
2027    const double X_22_5 = -(c13_p[5]*(du_p[INDEX3(0,0,5,numEq,3)] + du_p[INDEX3(1,1,5,numEq,3)]) + c33_p[5]*du_p[INDEX3(2,2,5,numEq,3)]);
2028    const double X_22_6 = -(c13_p[6]*(du_p[INDEX3(0,0,6,numEq,3)] + du_p[INDEX3(1,1,6,numEq,3)]) + c33_p[6]*du_p[INDEX3(2,2,6,numEq,3)]);
2029    const double X_22_7 = -(c13_p[7]*(du_p[INDEX3(0,0,7,numEq,3)] + du_p[INDEX3(1,1,7,numEq,3)]) + c33_p[7]*du_p[INDEX3(2,2,7,numEq,3)]);
2030    
2031    const double Atmp0 = w72*(X_00_6 + X_00_7);
2032    const double Atmp1 = w66*(X_02_0 + X_02_4);
2033    const double Atmp2 = w64*(X_00_0 + X_00_1);
2034    const double Atmp3 = w68*(X_02_1 + X_02_2 + X_02_5 + X_02_6);
2035    const double Atmp4 = w65*(X_01_0 + X_01_2);
2036    const double Atmp5 = w70*(X_02_3 + X_02_7);
2037    const double Atmp6 = w67*(X_01_1 + X_01_3 + X_01_4 + X_01_6);
2038    const double Atmp7 = w71*(X_01_5 + X_01_7);
2039    const double Atmp8 = w69*(X_00_2 + X_00_3 + X_00_4 + X_00_5);
2040    const double Atmp9 = w72*(-X_00_6 - X_00_7);
2041    const double Atmp10 = w66*(X_02_1 + X_02_5);
2042    const double Atmp11 = w64*(-X_00_0 - X_00_1);
2043    const double Atmp12 = w68*(X_02_0 + X_02_3 + X_02_4 + X_02_7);
2044    const double Atmp13 = w65*(X_01_1 + X_01_3);
2045    const double Atmp14 = w70*(X_02_2 + X_02_6);
2046    const double Atmp15 = w67*(X_01_0 + X_01_2 + X_01_5 + X_01_7);
2047    const double Atmp16 = w71*(X_01_4 + X_01_6);
2048    const double Atmp17 = w69*(-X_00_2 - X_00_3 - X_00_4 - X_00_5);
2049    const double Atmp18 = w72*(X_00_4 + X_00_5);
2050    const double Atmp19 = w66*(X_02_2 + X_02_6);
2051    const double Atmp20 = w64*(X_00_2 + X_00_3);
2052    const double Atmp21 = w65*(-X_01_0 - X_01_2);
2053    const double Atmp22 = w70*(X_02_1 + X_02_5);
2054    const double Atmp23 = w67*(-X_01_1 - X_01_3 - X_01_4 - X_01_6);
2055    const double Atmp24 = w71*(-X_01_5 - X_01_7);
2056    const double Atmp25 = w69*(X_00_0 + X_00_1 + X_00_6 + X_00_7);
2057    const double Atmp26 = w72*(-X_00_4 - X_00_5);
2058    const double Atmp27 = w66*(X_02_3 + X_02_7);
2059    const double Atmp28 = w64*(-X_00_2 - X_00_3);
2060    const double Atmp29 = w65*(-X_01_1 - X_01_3);
2061    const double Atmp30 = w70*(X_02_0 + X_02_4);
2062    const double Atmp31 = w67*(-X_01_0 - X_01_2 - X_01_5 - X_01_7);
2063    const double Atmp32 = w71*(-X_01_4 - X_01_6);
2064    const double Atmp33 = w69*(-X_00_0 - X_00_1 - X_00_6 - X_00_7);
2065    const double Atmp34 = w72*(X_00_2 + X_00_3);
2066    const double Atmp35 = w66*(-X_02_0 - X_02_4);
2067    const double Atmp36 = w64*(X_00_4 + X_00_5);
2068    const double Atmp37 = w68*(-X_02_1 - X_02_2 - X_02_5 - X_02_6);
2069    const double Atmp38 = w65*(X_01_4 + X_01_6);
2070    const double Atmp39 = w70*(-X_02_3 - X_02_7);
2071    const double Atmp40 = w71*(X_01_1 + X_01_3);
2072    const double Atmp41 = w72*(-X_00_2 - X_00_3);
2073    const double Atmp42 = w66*(-X_02_1 - X_02_5);
2074    const double Atmp43 = w64*(-X_00_4 - X_00_5);
2075    const double Atmp44 = w68*(-X_02_0 - X_02_3 - X_02_4 - X_02_7);
2076    const double Atmp45 = w65*(X_01_5 + X_01_7);
2077    const double Atmp46 = w70*(-X_02_2 - X_02_6);
2078    const double Atmp47 = w71*(X_01_0 + X_01_2);
2079    const double Atmp48 = w72*(X_00_0 + X_00_1);
2080    const double Atmp49 = w66*(-X_02_2 - X_02_6);
2081    const double Atmp50 = w64*(X_00_6 + X_00_7);
2082    const double Atmp51 = w65*(-X_01_4 - X_01_6);
2083    const double Atmp52 = w70*(-X_02_1 - X_02_5);
2084    const double Atmp53 = w71*(-X_01_1 - X_01_3);
2085    const double Atmp54 = w72*(-X_00_0 - X_00_1);
2086    const double Atmp55 = w66*(-X_02_3 - X_02_7);
2087    const double Atmp56 = w64*(-X_00_6 - X_00_7);
2088    const double Atmp57 = w65*(-X_01_5 - X_01_7);
2089    const double Atmp58 = w70*(-X_02_0 - X_02_4);
2090    const double Atmp59 = w71*(-X_01_0 - X_01_2);
2091    EM_F[INDEX2(0,0,numEq)]+=Atmp0 + Atmp1 + Atmp2 + Atmp3 + Atmp4 + Atmp5 + Atmp6 + Atmp7 + Atmp8;
2092    EM_F[INDEX2(0,1,numEq)]+=Atmp10 + Atmp11 + Atmp12 + Atmp13 + Atmp14 + Atmp15 + Atmp16 + Atmp17 + Atmp9;
2093    EM_F[INDEX2(0,2,numEq)]+=Atmp12 + Atmp18 + Atmp19 + Atmp20 + Atmp21 + Atmp22 + Atmp23 + Atmp24 + Atmp25;
2094    EM_F[INDEX2(0,3,numEq)]+=Atmp26 + Atmp27 + Atmp28 + Atmp29 + Atmp3 + Atmp30 + Atmp31 + Atmp32 + Atmp33;
2095    EM_F[INDEX2(0,4,numEq)]+=Atmp15 + Atmp25 + Atmp34 + Atmp35 + Atmp36 + Atmp37 + Atmp38 + Atmp39 + Atmp40;
2096    EM_F[INDEX2(0,5,numEq)]+=Atmp33 + Atmp41 + Atmp42 + Atmp43 + Atmp44 + Atmp45 + Atmp46 + Atmp47 + Atmp6;
2097    EM_F[INDEX2(0,6,numEq)]+=Atmp31 + Atmp44 + Atmp48 + Atmp49 + Atmp50 + Atmp51 + Atmp52 + Atmp53 + Atmp8;
2098    EM_F[INDEX2(0,7,numEq)]+=Atmp17 + Atmp23 + Atmp37 + Atmp54 + Atmp55 + Atmp56 + Atmp57 + Atmp58 + Atmp59;
2099    const double Btmp0 = w72*(X_01_6 + X_01_7);
2100    const double Btmp1 = w66*(X_12_0 + X_12_4);
2101    const double Btmp2 = w64*(X_01_0 + X_01_1);
2102    const double Btmp3 = w68*(X_12_1 + X_12_2 + X_12_5 + X_12_6);
2103    const double Btmp4 = w65*(X_11_0 + X_11_2);
2104    const double Btmp5 = w70*(X_12_3 + X_12_7);
2105    const double Btmp6 = w67*(X_11_1 + X_11_3 + X_11_4 + X_11_6);
2106    const double Btmp7 = w71*(X_11_5 + X_11_7);
2107    const double Btmp8 = w69*(X_01_2 + X_01_3 + X_01_4 + X_01_5);
2108    const double Btmp9 = w72*(-X_01_6 - X_01_7);
2109    const double Btmp10 = w66*(X_12_1 + X_12_5);
2110    const double Btmp11 = w64*(-X_01_0 - X_01_1);
2111    const double Btmp12 = w68*(X_12_0 + X_12_3 + X_12_4 + X_12_7);
2112    const double Btmp13 = w65*(X_11_1 + X_11_3);
2113    const double Btmp14 = w70*(X_12_2 + X_12_6);
2114    const double Btmp15 = w67*(X_11_0 + X_11_2 + X_11_5 + X_11_7);
2115    const double Btmp16 = w71*(X_11_4 + X_11_6);
2116    const double Btmp17 = w69*(-X_01_2 - X_01_3 - X_01_4 - X_01_5);
2117    const double Btmp18 = w72*(X_01_4 + X_01_5);
2118    const double Btmp19 = w66*(X_12_2 + X_12_6);
2119    const double Btmp20 = w64*(X_01_2 + X_01_3);
2120    const double Btmp21 = w65*(-X_11_0 - X_11_2);
2121    const double Btmp22 = w70*(X_12_1 + X_12_5);
2122    const double Btmp23 = w67*(-X_11_1 - X_11_3 - X_11_4 - X_11_6);
2123    const double Btmp24 = w71*(-X_11_5 - X_11_7);
2124    const double Btmp25 = w69*(X_01_0 + X_01_1 + X_01_6 + X_01_7);
2125    const double Btmp26 = w72*(-X_01_4 - X_01_5);
2126    const double Btmp27 = w66*(X_12_3 + X_12_7);
2127    const double Btmp28 = w64*(-X_01_2 - X_01_3);
2128    const double Btmp29 = w65*(-X_11_1 - X_11_3);
2129    const double Btmp30 = w70*(X_12_0 + X_12_4);
2130    const double Btmp31 = w67*(-X_11_0 - X_11_2 - X_11_5 - X_11_7);
2131    const double Btmp32 = w71*(-X_11_4 - X_11_6);
2132    const double Btmp33 = w69*(-X_01_0 - X_01_1 - X_01_6 - X_01_7);
2133    const double Btmp34 = w72*(X_01_2 + X_01_3);
2134    const double Btmp35 = w66*(-X_12_0 - X_12_4);
2135    const double Btmp36 = w64*(X_01_4 + X_01_5);
2136    const double Btmp37 = w68*(-X_12_1 - X_12_2 - X_12_5 - X_12_6);
2137    const double Btmp38 = w65*(X_11_4 + X_11_6);
2138    const double Btmp39 = w70*(-X_12_3 - X_12_7);
2139    const double Btmp40 = w71*(X_11_1 + X_11_3);
2140    const double Btmp41 = w72*(-X_01_2 - X_01_3);
2141    const double Btmp42 = w66*(-X_12_1 - X_12_5);
2142    const double Btmp43 = w64*(-X_01_4 - X_01_5);
2143    const double Btmp44 = w68*(-X_12_0 - X_12_3 - X_12_4 - X_12_7);
2144    const double Btmp45 = w65*(X_11_5 + X_11_7);
2145    const double Btmp46 = w70*(-X_12_2 - X_12_6);
2146    const double Btmp47 = w71*(X_11_0 + X_11_2);
2147    const double Btmp48 = w72*(X_01_0 + X_01_1);
2148    const double Btmp49 = w66*(-X_12_2 - X_12_6);
2149    const double Btmp50 = w64*(X_01_6 + X_01_7);
2150    const double Btmp51 = w65*(-X_11_4 - X_11_6);
2151    const double Btmp52 = w70*(-X_12_1 - X_12_5);
2152    const double Btmp53 = w71*(-X_11_1 - X_11_3);
2153    const double Btmp54 = w72*(-X_01_0 - X_01_1);
2154    const double Btmp55 = w66*(-X_12_3 - X_12_7);
2155    const double Btmp56 = w64*(-X_01_6 - X_01_7);
2156    const double Btmp57 = w65*(-X_11_5 - X_11_7);
2157    const double Btmp58 = w70*(-X_12_0 - X_12_4);
2158    const double Btmp59 = w71*(-X_11_0 - X_11_2);
2159    EM_F[INDEX2(1,0,numEq)]+=Btmp0 + Btmp1 + Btmp2 + Btmp3 + Btmp4 + Btmp5 + Btmp6 + Btmp7 + Btmp8;
2160    EM_F[INDEX2(1,1,numEq)]+=Btmp10 + Btmp11 + Btmp12 + Btmp13 + Btmp14 + Btmp15 + Btmp16 + Btmp17 + Btmp9;
2161    EM_F[INDEX2(1,2,numEq)]+=Btmp12 + Btmp18 + Btmp19 + Btmp20 + Btmp21 + Btmp22 + Btmp23 + Btmp24 + Btmp25;
2162    EM_F[INDEX2(1,3,numEq)]+=Btmp26 + Btmp27 + Btmp28 + Btmp29 + Btmp3 + Btmp30 + Btmp31 + Btmp32 + Btmp33;
2163    EM_F[INDEX2(1,4,numEq)]+=Btmp15 + Btmp25 + Btmp34 + Btmp35 + Btmp36 + Btmp37 + Btmp38 + Btmp39 + Btmp40;
2164    EM_F[INDEX2(1,5,numEq)]+=Btmp33 + Btmp41 + Btmp42 + Btmp43 + Btmp44 + Btmp45 + Btmp46 + Btmp47 + Btmp6;
2165    EM_F[INDEX2(1,6,numEq)]+=Btmp31 + Btmp44 + Btmp48 + Btmp49 + Btmp50 + Btmp51 + Btmp52 + Btmp53 + Btmp8;
2166    EM_F[INDEX2(1,7,numEq)]+=Btmp17 + Btmp23 + Btmp37 + Btmp54 + Btmp55 + Btmp56 + Btmp57 + Btmp58 + Btmp59;
2167    const double Ctmp0 = w72*(X_02_6 + X_02_7);
2168    const double Ctmp1 = w66*(X_22_0 + X_22_4);
2169    const double Ctmp2 = w64*(X_02_0 + X_02_1);
2170    const double Ctmp3 = w68*(X_22_1 + X_22_2 + X_22_5 + X_22_6);
2171    const double Ctmp4 = w65*(X_12_0 + X_12_2);
2172    const double Ctmp5 = w70*(X_22_3 + X_22_7);
2173    const double Ctmp6 = w67*(X_12_1 + X_12_3 + X_12_4 + X_12_6);
2174    const double Ctmp7 = w71*(X_12_5 + X_12_7);
2175    const double Ctmp8 = w69*(X_02_2 + X_02_3 + X_02_4 + X_02_5);
2176    const double Ctmp9 = w72*(-X_02_6 - X_02_7);
2177    const double Ctmp10 = w66*(X_22_1 + X_22_5);
2178    const double Ctmp11 = w64*(-X_02_0 - X_02_1);
2179    const double Ctmp12 = w68*(X_22_0 + X_22_3 + X_22_4 + X_22_7);
2180    const double Ctmp13 = w65*(X_12_1 + X_12_3);
2181    const double Ctmp14 = w70*(X_22_2 + X_22_6);
2182    const double Ctmp15 = w67*(X_12_0 + X_12_2 + X_12_5 + X_12_7);
2183    const double Ctmp16 = w71*(X_12_4 + X_12_6);
2184    const double Ctmp17 = w69*(-X_02_2 - X_02_3 - X_02_4 - X_02_5);
2185    const double Ctmp18 = w72*(X_02_4 + X_02_5);
2186    const double Ctmp19 = w66*(X_22_2 + X_22_6);
2187    const double Ctmp20 = w64*(X_02_2 + X_02_3);
2188    const double Ctmp21 = w65*(-X_12_0 - X_12_2);
2189    const double Ctmp22 = w70*(X_22_1 + X_22_5);
2190    const double Ctmp23 = w67*(-X_12_1 - X_12_3 - X_12_4 - X_12_6);
2191    const double Ctmp24 = w71*(-X_12_5 - X_12_7);
2192    const double Ctmp25 = w69*(X_02_0 + X_02_1 + X_02_6 + X_02_7);
2193    const double Ctmp26 = w72*(-X_02_4 - X_02_5);
2194    const double Ctmp27 = w66*(X_22_3 + X_22_7);
2195    const double Ctmp28 = w64*(-X_02_2 - X_02_3);
2196    const double Ctmp29 = w65*(-X_12_1 - X_12_3);
2197    const double Ctmp30 = w70*(X_22_0 + X_22_4);
2198    const double Ctmp31 = w67*(-X_12_0 - X_12_2 - X_12_5 - X_12_7);
2199    const double Ctmp32 = w71*(-X_12_4 - X_12_6);
2200    const double Ctmp33 = w69*(-X_02_0 - X_02_1 - X_02_6 - X_02_7);
2201    const double Ctmp34 = w72*(X_02_2 + X_02_3);
2202    const double Ctmp35 = w66*(-X_22_0 - X_22_4);
2203    const double Ctmp36 = w64*(X_02_4 + X_02_5);
2204    const double Ctmp37 = w68*(-X_22_1 - X_22_2 - X_22_5 - X_22_6);
2205    const double Ctmp38 = w65*(X_12_4 + X_12_6);
2206    const double Ctmp39 = w70*(-X_22_3 - X_22_7);
2207    const double Ctmp40 = w71*(X_12_1 + X_12_3);
2208    const double Ctmp41 = w72*(-X_02_2 - X_02_3);
2209    const double Ctmp42 = w66*(-X_22_1 - X_22_5);
2210    const double Ctmp43 = w64*(-X_02_4 - X_02_5);
2211    const double Ctmp44 = w68*(-X_22_0 - X_22_3 - X_22_4 - X_22_7);
2212    const double Ctmp45 = w65*(X_12_5 + X_12_7);
2213    const double Ctmp46 = w70*(-X_22_2 - X_22_6);
2214    const double Ctmp47 = w71*(X_12_0 + X_12_2);
2215    const double Ctmp48 = w72*(X_02_0 + X_02_1);
2216    const double Ctmp49 = w66*(-X_22_2 - X_22_6);
2217    const double Ctmp50 = w64*(X_02_6 + X_02_7);
2218    const double Ctmp51 = w65*(-X_12_4 - X_12_6);
2219    const double Ctmp52 = w70*(-X_22_1 - X_22_5);
2220    const double Ctmp53 = w71*(-X_12_1 - X_12_3);
2221    const double Ctmp54 = w72*(-X_02_0 - X_02_1);
2222    const double Ctmp55 = w66*(-X_22_3 - X_22_7);
2223    const double Ctmp56 = w64*(-X_02_6 - X_02_7);
2224    const double Ctmp57 = w65*(-X_12_5 - X_12_7);
2225    const double Ctmp58 = w70*(-X_22_0 - X_22_4);
2226    const double Ctmp59 = w71*(-X_12_0 - X_12_2);
2227    EM_F[INDEX2(2,0,numEq)]+=Ctmp0 + Ctmp1 + Ctmp2 + Ctmp3 + Ctmp4 + Ctmp5 + Ctmp6 + Ctmp7 + Ctmp8;
2228    EM_F[INDEX2(2,1,numEq)]+=Ctmp10 + Ctmp11 + Ctmp12 + Ctmp13 + Ctmp14 + Ctmp15 + Ctmp16 + Ctmp17 + Ctmp9;
2229    EM_F[INDEX2(2,2,numEq)]+=Ctmp12 + Ctmp18 + Ctmp19 + Ctmp20 + Ctmp21 + Ctmp22 + Ctmp23 + Ctmp24 + Ctmp25;
2230    EM_F[INDEX2(2,3,numEq)]+=Ctmp26 + Ctmp27 + Ctmp28 + Ctmp29 + Ctmp3 + Ctmp30 + Ctmp31 + Ctmp32 + Ctmp33;
2231    EM_F[INDEX2(2,4,numEq)]+=Ctmp15 + Ctmp25 + Ctmp34 + Ctmp35 + Ctmp36 + Ctmp37 + Ctmp38 + Ctmp39 + Ctmp40;
2232    EM_F[INDEX2(2,5,numEq)]+=Ctmp33 + Ctmp41 + Ctmp42 + Ctmp43 + Ctmp44 + Ctmp45 + Ctmp46 + Ctmp47 + Ctmp6;
2233    EM_F[INDEX2(2,6,numEq)]+=Ctmp31 + Ctmp44 + Ctmp48 + Ctmp49 + Ctmp50 + Ctmp51 + Ctmp52 + Ctmp53 + Ctmp8;
2234    EM_F[INDEX2(2,7,numEq)]+=Ctmp17 + Ctmp23 + Ctmp37 + Ctmp54 + Ctmp55 + Ctmp56 + Ctmp57 + Ctmp58 + Ctmp59;
2235    #if 0
2236    #define NUM_QUADPOINTS 8
2237                                double *X_00[NUM_QUADPOINTS];
2238                                double *X_01[NUM_QUADPOINTS];
2239                                double *X_02[NUM_QUADPOINTS];
2240                                double *X_11[NUM_QUADPOINTS];
2241                                double *X_12[NUM_QUADPOINTS];
2242                                double *X_22[NUM_QUADPOINTS];
2243    
2244  #define GET(a,b) [INDEX3((a),(b),qp,numEq,3)]  #define GET(a,b) [INDEX3((a),(b),qp,numEq,3)]
2245                                  for (int qp = 0; qp < NUM_QUADPOINTS; qp++) {                                  for (int qp = 0; qp < NUM_QUADPOINTS; qp++) {
# Line 2101  X_12[qp] = -(c44_p[qp]*(du_p GET(2,1) + Line 2355  X_12[qp] = -(c44_p[qp]*(du_p GET(2,1) +
2355                                      EM_F[INDEX2(k,5,numEq)]+=tmp33 + tmp41 + tmp42 + tmp43 + tmp44 + tmp45 + tmp46 + tmp47 + tmp6;                                      EM_F[INDEX2(k,5,numEq)]+=tmp33 + tmp41 + tmp42 + tmp43 + tmp44 + tmp45 + tmp46 + tmp47 + tmp6;
2356                                      EM_F[INDEX2(k,6,numEq)]+=tmp31 + tmp44 + tmp48 + tmp49 + tmp50 + tmp51 + tmp52 + tmp53 + tmp8;                                      EM_F[INDEX2(k,6,numEq)]+=tmp31 + tmp44 + tmp48 + tmp49 + tmp50 + tmp51 + tmp52 + tmp53 + tmp8;
2357                                      EM_F[INDEX2(k,7,numEq)]+=tmp17 + tmp23 + tmp37 + tmp54 + tmp55 + tmp56 + tmp57 + tmp58 + tmp59;                                      EM_F[INDEX2(k,7,numEq)]+=tmp17 + tmp23 + tmp37 + tmp54 + tmp55 + tmp56 + tmp57 + tmp58 + tmp59;
                                 }                                
                             } else { // constant data  
                                 throw RipleyException("WaveAssembler3D support"  
                                     " for constant data not yet implemented");  
                                 /*  
                                 for (index_t k=0; k<numEq; k++) {  
                                     const double wX0 = 18*X_p[INDEX2(k, 0, numEq)]*w55;  
                                     const double wX1 = 18*X_p[INDEX2(k, 1, numEq)]*w56;  
                                     const double wX2 = 18*X_p[INDEX2(k, 2, numEq)]*w54;  
                                     EM_F[INDEX2(k,0,numEq)]+= wX0 + wX1 + wX2;  
                                     EM_F[INDEX2(k,1,numEq)]+=-wX0 + wX1 + wX2;  
                                     EM_F[INDEX2(k,2,numEq)]+= wX0 - wX1 + wX2;  
                                     EM_F[INDEX2(k,3,numEq)]+=-wX0 - wX1 + wX2;  
                                     EM_F[INDEX2(k,4,numEq)]+= wX0 + wX1 - wX2;  
                                     EM_F[INDEX2(k,5,numEq)]+=-wX0 + wX1 - wX2;  
                                     EM_F[INDEX2(k,6,numEq)]+= wX0 - wX1 - wX2;  
                                     EM_F[INDEX2(k,7,numEq)]+=-wX0 - wX1 - wX2;  
2358                                  }                                  }
2359                                  */  #endif
2360                                } else { // constant data
2361                                    const double wX00 = 18*w55*-(c11_p[0]* du_p[INDEX2(0,0,numEq)] + c12_p[0] * du_p[INDEX2(1,1,numEq)] + c13_p[0]*du_p[INDEX2(2,2,numEq)]);
2362                                    const double wX11 = 18*w56*-(c12_p[0]* du_p[INDEX2(0,0,numEq)] + c11_p[0] * du_p[INDEX2(1,1,numEq)] + c13_p[0]*du_p[INDEX2(2,2,numEq)]);
2363                                    const double wX22 = 18*w54*-(c13_p[0]*(du_p[INDEX2(0,0,numEq)] + du_p[INDEX2(1,1,numEq)]) + c33_p[0]* du_p[INDEX2(2,2,numEq)]);
2364                                    
2365                                    const double wX01 = 18*w56*-(c66_p[0]*(du_p[INDEX2(0,1,numEq)] + du_p[INDEX2(1,0,numEq)]));
2366                                    const double wX10 = 18*w55*-(c66_p[0]*(du_p[INDEX2(0,1,numEq)] + du_p[INDEX2(1,0,numEq)]));
2367    
2368                                    const double wX02 = 18*w54*-(c44_p[0]*(du_p[INDEX2(2,0,numEq)] + du_p[INDEX2(0,2,numEq)]));
2369                                    const double wX20 = 18*w55*-(c44_p[0]*(du_p[INDEX2(2,0,numEq)] + du_p[INDEX2(0,2,numEq)]));
2370    
2371                                    const double wX12 = 18*w54*-(c44_p[0]*(du_p[INDEX2(2,1,numEq)] + du_p[INDEX2(1,2,numEq)]));
2372                                    const double wX21 = 18*w56*-(c44_p[0]*(du_p[INDEX2(2,1,numEq)] + du_p[INDEX2(1,2,numEq)]));
2373                                    EM_F[INDEX2(0,0,numEq)]+= wX00 + wX01 + wX02;
2374                                    EM_F[INDEX2(0,1,numEq)]+=-wX00 + wX01 + wX02;
2375                                    EM_F[INDEX2(0,2,numEq)]+= wX00 - wX01 + wX02;
2376                                    EM_F[INDEX2(0,3,numEq)]+=-wX00 - wX01 + wX02;
2377                                    EM_F[INDEX2(0,4,numEq)]+= wX00 + wX01 - wX02;
2378                                    EM_F[INDEX2(0,5,numEq)]+=-wX00 + wX01 - wX02;
2379                                    EM_F[INDEX2(0,6,numEq)]+= wX00 - wX01 - wX02;
2380                                    EM_F[INDEX2(0,7,numEq)]+=-wX00 - wX01 - wX02;
2381                                    
2382                                    EM_F[INDEX2(1,0,numEq)]+= wX10 + wX11 + wX12;
2383                                    EM_F[INDEX2(1,1,numEq)]+=-wX10 + wX11 + wX12;
2384                                    EM_F[INDEX2(1,2,numEq)]+= wX10 - wX11 + wX12;
2385                                    EM_F[INDEX2(1,3,numEq)]+=-wX10 - wX11 + wX12;
2386                                    EM_F[INDEX2(1,4,numEq)]+= wX10 + wX11 - wX12;
2387                                    EM_F[INDEX2(1,5,numEq)]+=-wX10 + wX11 - wX12;
2388                                    EM_F[INDEX2(1,6,numEq)]+= wX10 - wX11 - wX12;
2389                                    EM_F[INDEX2(1,7,numEq)]+=-wX10 - wX11 - wX12;
2390                                    
2391                                    EM_F[INDEX2(2,0,numEq)]+= wX20 + wX21 + wX22;
2392                                    EM_F[INDEX2(2,1,numEq)]+=-wX20 + wX21 + wX22;
2393                                    EM_F[INDEX2(2,2,numEq)]+= wX20 - wX21 + wX22;
2394                                    EM_F[INDEX2(2,3,numEq)]+=-wX20 - wX21 + wX22;
2395                                    EM_F[INDEX2(2,4,numEq)]+= wX20 + wX21 - wX22;
2396                                    EM_F[INDEX2(2,5,numEq)]+=-wX20 + wX21 - wX22;
2397                                    EM_F[INDEX2(2,6,numEq)]+= wX20 - wX21 - wX22;
2398                                    EM_F[INDEX2(2,7,numEq)]+=-wX20 - wX21 - wX22;
2399                              }                              }
                                 delete[] X_00;  
                                 delete[] X_01;  
                                 delete[] X_02;  
                                 delete[] X_11;  
                                 delete[] X_12;  
                                 delete[] X_22;  
2400                          }                          }
2401                          ///////////////                          ///////////////
2402                          // process Y //                          // process Y //

Legend:
Removed from v.4644  
changed lines
  Added in v.4645

  ViewVC Help
Powered by ViewVC 1.1.26