/[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 3761 by caltinay, Mon Jan 9 06:50:33 2012 UTC revision 3762 by caltinay, Tue Jan 10 00:01:45 2012 UTC
# Line 2178  void Rectangle::assemblePDESingle(Paso_S Line 2178  void Rectangle::assemblePDESingle(Paso_S
2178  }  }
2179    
2180  //protected  //protected
2181    void Rectangle::assemblePDESingleReduced(Paso_SystemMatrix* mat,
2182            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2183            const escript::Data& C, const escript::Data& D,
2184            const escript::Data& X, const escript::Data& Y,
2185            const escript::Data& d, const escript::Data& y) const
2186    {
2187        const double h0 = m_l0/m_gNE0;
2188        const double h1 = m_l1/m_gNE1;
2189        /*** GENERATOR SNIP_PDE_SINGLE_REDUCED_PRE TOP */
2190        const double w0 = -0.25*h1/h0;
2191        const double w1 = 0.25;
2192        const double w2 = -0.25;
2193        const double w3 = 0.25*h0/h1;
2194        const double w4 = -0.25*h0/h1;
2195        const double w5 = 0.25*h1/h0;
2196        const double w6 = -0.125*h1;
2197        const double w7 = -0.125*h0;
2198        const double w8 = 0.125*h1;
2199        const double w9 = 0.125*h0;
2200        const double w10 = 0.0625*h0*h1;
2201        const double w11 = -0.5*h1;
2202        const double w12 = -0.5*h0;
2203        const double w13 = 0.5*h1;
2204        const double w14 = 0.5*h0;
2205        const double w15 = 0.25*h0*h1;
2206        /* GENERATOR SNIP_PDE_SINGLE_REDUCED_PRE BOTTOM */
2207    
2208        rhs.requireWrite();
2209    #pragma omp parallel
2210        {
2211            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
2212    #pragma omp for
2213                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
2214                    for (index_t k0=0; k0<m_NE0; ++k0)  {
2215                        bool add_EM_S=false;
2216                        bool add_EM_F=false;
2217                        vector<double> EM_S(4*4, 0);
2218                        vector<double> EM_F(4, 0);
2219                        const index_t e = k0 + m_NE0*k1;
2220                        /*** GENERATOR SNIP_PDE_SINGLE_REDUCED TOP */
2221                        ///////////////
2222                        // process A //
2223                        ///////////////
2224                        if (!A.isEmpty()) {
2225                            add_EM_S=true;
2226                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
2227                            const register double A_00 = A_p[INDEX2(0,0,2)];
2228                            const register double A_01 = A_p[INDEX2(0,1,2)];
2229                            const register double A_10 = A_p[INDEX2(1,0,2)];
2230                            const register double A_11 = A_p[INDEX2(1,1,2)];
2231                            const register double tmp0_0 = A_01 + A_10;
2232                            const register double tmp2_1 = A_01*w1;
2233                            const register double tmp7_1 = tmp0_0*w2;
2234                            const register double tmp3_1 = A_10*w2;
2235                            const register double tmp6_1 = A_00*w5;
2236                            const register double tmp1_1 = A_00*w0;
2237                            const register double tmp4_1 = tmp0_0*w1;
2238                            const register double tmp9_1 = A_01*w2;
2239                            const register double tmp5_1 = A_11*w4;
2240                            const register double tmp8_1 = A_10*w1;
2241                            const register double tmp0_1 = A_11*w3;
2242                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2243                            EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2244                            EM_S[INDEX2(3,2,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
2245                            EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2246                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
2247                            EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2248                            EM_S[INDEX2(3,1,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2249                            EM_S[INDEX2(2,1,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
2250                            EM_S[INDEX2(0,2,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
2251                            EM_S[INDEX2(2,0,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2252                            EM_S[INDEX2(1,3,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
2253                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2254                            EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2255                            EM_S[INDEX2(1,0,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
2256                            EM_S[INDEX2(0,3,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
2257                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
2258                        }
2259                        ///////////////
2260                        // process B //
2261                        ///////////////
2262                        if (!B.isEmpty()) {
2263                            add_EM_S=true;
2264                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
2265                            const register double B_0 = B_p[0];
2266                            const register double B_1 = B_p[1];
2267                            const register double tmp2_1 = B_0*w8;
2268                            const register double tmp0_1 = B_0*w6;
2269                            const register double tmp3_1 = B_1*w9;
2270                            const register double tmp1_1 = B_1*w7;
2271                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2272                            EM_S[INDEX2(1,2,4)]+=tmp1_1 + tmp2_1;
2273                            EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2274                            EM_S[INDEX2(0,0,4)]+=tmp0_1 + tmp1_1;
2275                            EM_S[INDEX2(3,3,4)]+=tmp2_1 + tmp3_1;
2276                            EM_S[INDEX2(3,0,4)]+=tmp2_1 + tmp3_1;
2277                            EM_S[INDEX2(3,1,4)]+=tmp2_1 + tmp3_1;
2278                            EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp3_1;
2279                            EM_S[INDEX2(0,2,4)]+=tmp0_1 + tmp1_1;
2280                            EM_S[INDEX2(2,0,4)]+=tmp0_1 + tmp3_1;
2281                            EM_S[INDEX2(1,3,4)]+=tmp1_1 + tmp2_1;
2282                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp3_1;
2283                            EM_S[INDEX2(2,2,4)]+=tmp0_1 + tmp3_1;
2284                            EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp2_1;
2285                            EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp1_1;
2286                            EM_S[INDEX2(1,1,4)]+=tmp1_1 + tmp2_1;
2287                        }
2288                        ///////////////
2289                        // process C //
2290                        ///////////////
2291                        if (!C.isEmpty()) {
2292                            add_EM_S=true;
2293                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
2294                            const register double C_0 = C_p[0];
2295                            const register double C_1 = C_p[1];
2296                            const register double tmp1_1 = C_1*w7;
2297                            const register double tmp0_1 = C_0*w8;
2298                            const register double tmp3_1 = C_0*w6;
2299                            const register double tmp2_1 = C_1*w9;
2300                            EM_S[INDEX2(0,1,4)]+=tmp0_1 + tmp1_1;
2301                            EM_S[INDEX2(1,2,4)]+=tmp2_1 + tmp3_1;
2302                            EM_S[INDEX2(3,2,4)]+=tmp2_1 + tmp3_1;
2303                            EM_S[INDEX2(0,0,4)]+=tmp1_1 + tmp3_1;
2304                            EM_S[INDEX2(3,3,4)]+=tmp0_1 + tmp2_1;
2305                            EM_S[INDEX2(3,0,4)]+=tmp1_1 + tmp3_1;
2306                            EM_S[INDEX2(3,1,4)]+=tmp0_1 + tmp1_1;
2307                            EM_S[INDEX2(2,1,4)]+=tmp0_1 + tmp1_1;
2308                            EM_S[INDEX2(0,2,4)]+=tmp2_1 + tmp3_1;
2309                            EM_S[INDEX2(2,0,4)]+=tmp1_1 + tmp3_1;
2310                            EM_S[INDEX2(1,3,4)]+=tmp0_1 + tmp2_1;
2311                            EM_S[INDEX2(2,3,4)]+=tmp0_1 + tmp2_1;
2312                            EM_S[INDEX2(2,2,4)]+=tmp2_1 + tmp3_1;
2313                            EM_S[INDEX2(1,0,4)]+=tmp1_1 + tmp3_1;
2314                            EM_S[INDEX2(0,3,4)]+=tmp0_1 + tmp2_1;
2315                            EM_S[INDEX2(1,1,4)]+=tmp0_1 + tmp1_1;
2316                        }
2317                        ///////////////
2318                        // process D //
2319                        ///////////////
2320                        if (!D.isEmpty()) {
2321                            add_EM_S=true;
2322                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
2323                            const register double D_0 = D_p[0];
2324                            const register double tmp0_1 = D_0*w10;
2325                            EM_S[INDEX2(0,1,4)]+=tmp0_1;
2326                            EM_S[INDEX2(1,2,4)]+=tmp0_1;
2327                            EM_S[INDEX2(3,2,4)]+=tmp0_1;
2328                            EM_S[INDEX2(0,0,4)]+=tmp0_1;
2329                            EM_S[INDEX2(3,3,4)]+=tmp0_1;
2330                            EM_S[INDEX2(3,0,4)]+=tmp0_1;
2331                            EM_S[INDEX2(3,1,4)]+=tmp0_1;
2332                            EM_S[INDEX2(2,1,4)]+=tmp0_1;
2333                            EM_S[INDEX2(0,2,4)]+=tmp0_1;
2334                            EM_S[INDEX2(2,0,4)]+=tmp0_1;
2335                            EM_S[INDEX2(1,3,4)]+=tmp0_1;
2336                            EM_S[INDEX2(2,3,4)]+=tmp0_1;
2337                            EM_S[INDEX2(2,2,4)]+=tmp0_1;
2338                            EM_S[INDEX2(1,0,4)]+=tmp0_1;
2339                            EM_S[INDEX2(0,3,4)]+=tmp0_1;
2340                            EM_S[INDEX2(1,1,4)]+=tmp0_1;
2341                        }
2342                        ///////////////
2343                        // process X //
2344                        ///////////////
2345                        if (!X.isEmpty()) {
2346                            add_EM_F=true;
2347                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
2348                            const register double X_0 = X_p[0];
2349                            const register double X_1 = X_p[1];
2350                            const register double tmp0_1 = X_0*w11;
2351                            const register double tmp2_1 = X_0*w13;
2352                            const register double tmp1_1 = X_1*w12;
2353                            const register double tmp3_1 = X_1*w14;
2354                            EM_F[0]+=tmp0_1 + tmp1_1;
2355                            EM_F[1]+=tmp1_1 + tmp2_1;
2356                            EM_F[2]+=tmp0_1 + tmp3_1;
2357                            EM_F[3]+=tmp2_1 + tmp3_1;
2358                        }
2359                        ///////////////
2360                        // process Y //
2361                        ///////////////
2362                        if (!Y.isEmpty()) {
2363                            add_EM_F=true;
2364                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
2365                            const register double Y_0 = Y_p[0];
2366                            const register double tmp0_1 = Y_0*w15;
2367                            EM_F[0]+=tmp0_1;
2368                            EM_F[1]+=tmp0_1;
2369                            EM_F[2]+=tmp0_1;
2370                            EM_F[3]+=tmp0_1;
2371                        }
2372                        /* GENERATOR SNIP_PDE_SINGLE_REDUCED BOTTOM */
2373    
2374                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2375                        const index_t firstNode=m_N0*k1+k0;
2376                        IndexVector rowIndex;
2377                        rowIndex.push_back(m_dofMap[firstNode]);
2378                        rowIndex.push_back(m_dofMap[firstNode+1]);
2379                        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
2380                        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
2381                        if (add_EM_F) {
2382                            //cout << "-- AddtoRHS -- " << endl;
2383                            double *F_p=rhs.getSampleDataRW(0);
2384                            for (index_t i=0; i<rowIndex.size(); i++) {
2385                                if (rowIndex[i]<getNumDOF()) {
2386                                    F_p[rowIndex[i]]+=EM_F[i];
2387                                    //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;
2388                                }
2389                            }
2390                            //cout << "---"<<endl;
2391                        }
2392                        if (add_EM_S) {
2393                            //cout << "-- AddtoSystem -- " << endl;
2394                            addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
2395                        }
2396    
2397                    } // end k0 loop
2398                } // end k1 loop
2399            } // end of colouring
2400        } // end of parallel region
2401    }
2402    
2403    //protected
2404  void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESystem(Paso_SystemMatrix* mat,
2405          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
2406          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
# Line 2937  void Rectangle::assemblePDESystem(Paso_S Line 3160  void Rectangle::assemblePDESystem(Paso_S
3160    
3161                      // 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)
3162                      const index_t firstNode=m_N0*k1+k0;                      const index_t firstNode=m_N0*k1+k0;
3163                        IndexVector rowIndex;
3164                        rowIndex.push_back(m_dofMap[firstNode]);
3165                        rowIndex.push_back(m_dofMap[firstNode+1]);
3166                        rowIndex.push_back(m_dofMap[firstNode+m_N0]);
3167                        rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
3168                        if (add_EM_F) {
3169                            //cout << "-- AddtoRHS -- " << endl;
3170                            double *F_p=rhs.getSampleDataRW(0);
3171                            for (index_t i=0; i<rowIndex.size(); i++) {
3172                                if (rowIndex[i]<getNumDOF()) {
3173                                    for (index_t eq=0; eq<numEq; eq++) {
3174                                        F_p[INDEX2(eq,rowIndex[i],numEq)]+=EM_F[INDEX2(eq,i,numEq)];
3175                                        //cout << "F[" << INDEX2(eq,rowIndex[i],numEq) << "]=" << F_p[INDEX2(eq,rowIndex[i],numEq)] << endl;
3176                                    }
3177                                }
3178                            }
3179                            //cout << "---"<<endl;
3180                        }
3181                        if (add_EM_S) {
3182                            //cout << "-- AddtoSystem -- " << endl;
3183                            addToSystemMatrix(mat, rowIndex, numEq, rowIndex, numComp, EM_S);
3184                        }
3185    
3186                    } // end k0 loop
3187                } // end k1 loop
3188            } // end of colouring
3189        } // end of parallel region
3190    }
3191    
3192    //protected
3193    void Rectangle::assemblePDESystemReduced(Paso_SystemMatrix* mat,
3194            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
3195            const escript::Data& C, const escript::Data& D,
3196            const escript::Data& X, const escript::Data& Y,
3197            const escript::Data& d, const escript::Data& y) const
3198    {
3199        const double h0 = m_l0/m_gNE0;
3200        const double h1 = m_l1/m_gNE1;
3201        dim_t numEq, numComp;
3202        if (!mat)
3203            numEq=numComp=(rhs.isEmpty() ? 1 : rhs.getDataPointSize());
3204        else {
3205            numEq=mat->logical_row_block_size;
3206            numComp=mat->logical_col_block_size;
3207        }
3208    
3209        /*** GENERATOR SNIP_PDE_SYSTEM_REDUCED_PRE TOP */
3210        const double w0 = -0.25*h1/h0;
3211        const double w1 = 0.25;
3212        const double w2 = -0.25;
3213        const double w3 = 0.25*h0/h1;
3214        const double w4 = -0.25*h0/h1;
3215        const double w5 = 0.25*h1/h0;
3216        const double w6 = -0.125*h1;
3217        const double w7 = -0.125*h0;
3218        const double w8 = 0.125*h1;
3219        const double w9 = 0.125*h0;
3220        const double w10 = 0.0625*h0*h1;
3221        const double w11 = -0.5*h1;
3222        const double w12 = -0.5*h0;
3223        const double w13 = 0.5*h1;
3224        const double w14 = 0.5*h0;
3225        const double w15 = 0.25*h0*h1;
3226        /* GENERATOR SNIP_PDE_SYSTEM_REDUCED_PRE BOTTOM */
3227    
3228        rhs.requireWrite();
3229    #pragma omp parallel
3230        {
3231            for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3232    #pragma omp for
3233                for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3234                    for (index_t k0=0; k0<m_NE0; ++k0)  {
3235                        bool add_EM_S=false;
3236                        bool add_EM_F=false;
3237                        vector<double> EM_S(4*4*numEq*numComp, 0);
3238                        vector<double> EM_F(4*numEq, 0);
3239                        const index_t e = k0 + m_NE0*k1;
3240                        /*** GENERATOR SNIP_PDE_SYSTEM_REDUCED TOP */
3241                        ///////////////
3242                        // process A //
3243                        ///////////////
3244                        if (!A.isEmpty()) {
3245                            add_EM_S=true;
3246                            const double* A_p=const_cast<escript::Data*>(&A)->getSampleDataRO(e);
3247                            for (index_t k=0; k<numEq; k++) {
3248                                for (index_t m=0; m<numComp; m++) {
3249                                    const register double A_00 = A_p[INDEX4(k,0,m,0, numEq,2, numComp)];
3250                                    const register double A_01 = A_p[INDEX4(k,0,m,1, numEq,2, numComp)];
3251                                    const register double A_10 = A_p[INDEX4(k,1,m,0, numEq,2, numComp)];
3252                                    const register double A_11 = A_p[INDEX4(k,1,m,1, numEq,2, numComp)];
3253                                    const register double tmp0_0 = A_01 + A_10;
3254                                    const register double tmp2_1 = A_01*w1;
3255                                    const register double tmp7_1 = tmp0_0*w2;
3256                                    const register double tmp3_1 = A_10*w2;
3257                                    const register double tmp6_1 = A_00*w5;
3258                                    const register double tmp1_1 = A_00*w0;
3259                                    const register double tmp4_1 = tmp0_0*w1;
3260                                    const register double tmp9_1 = A_01*w2;
3261                                    const register double tmp5_1 = A_11*w4;
3262                                    const register double tmp8_1 = A_10*w1;
3263                                    const register double tmp0_1 = A_11*w3;
3264                                    EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
3265                                    EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
3266                                    EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
3267                                    EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
3268                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
3269                                    EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
3270                                    EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
3271                                    EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
3272                                    EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
3273                                    EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
3274                                    EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
3275                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
3276                                    EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
3277                                    EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
3278                                    EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
3279                                    EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
3280                                }
3281                            }
3282                        }
3283                        ///////////////
3284                        // process B //
3285                        ///////////////
3286                        if (!B.isEmpty()) {
3287                            add_EM_S=true;
3288                            const double* B_p=const_cast<escript::Data*>(&B)->getSampleDataRO(e);
3289                            for (index_t k=0; k<numEq; k++) {
3290                                for (index_t m=0; m<numComp; m++) {
3291                                    const register double B_0 = B_p[INDEX3(k,0,m, numEq, 2)];
3292                                    const register double B_1 = B_p[INDEX3(k,1,m, numEq, 2)];
3293                                    const register double tmp2_1 = B_0*w8;
3294                                    const register double tmp0_1 = B_0*w6;
3295                                    const register double tmp3_1 = B_1*w9;
3296                                    const register double tmp1_1 = B_1*w7;
3297                                    EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3298                                    EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3299                                    EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3300                                    EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3301                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3302                                    EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3303                                    EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3304                                    EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
3305                                    EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3306                                    EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
3307                                    EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3308                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
3309                                    EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1 + tmp3_1;
3310                                    EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3311                                    EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3312                                    EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp1_1 + tmp2_1;
3313                                }
3314                            }
3315                        }
3316                        ///////////////
3317                        // process C //
3318                        ///////////////
3319                        if (!C.isEmpty()) {
3320                            add_EM_S=true;
3321                            const double* C_p=const_cast<escript::Data*>(&C)->getSampleDataRO(e);
3322                            for (index_t k=0; k<numEq; k++) {
3323                                for (index_t m=0; m<numComp; m++) {
3324                                    const register double C_0 = C_p[INDEX3(k, m, 0, numEq, numComp)];
3325                                    const register double C_1 = C_p[INDEX3(k, m, 1, numEq, numComp)];
3326                                    const register double tmp1_1 = C_1*w7;
3327                                    const register double tmp0_1 = C_0*w8;
3328                                    const register double tmp3_1 = C_0*w6;
3329                                    const register double tmp2_1 = C_1*w9;
3330                                    EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3331                                    EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3332                                    EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3333                                    EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3334                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
3335                                    EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3336                                    EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3337                                    EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3338                                    EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3339                                    EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3340                                    EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
3341                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
3342                                    EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp2_1 + tmp3_1;
3343                                    EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp1_1 + tmp3_1;
3344                                    EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1 + tmp2_1;
3345                                    EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1 + tmp1_1;
3346                                }
3347                            }
3348                        }
3349                        ///////////////
3350                        // process D //
3351                        ///////////////
3352                        if (!D.isEmpty()) {
3353                            add_EM_S=true;
3354                            const double* D_p=const_cast<escript::Data*>(&D)->getSampleDataRO(e);
3355                            for (index_t k=0; k<numEq; k++) {
3356                                for (index_t m=0; m<numComp; m++) {
3357                                    const register double D_0 = D_p[INDEX2(k, m, numEq)];
3358                                    const register double tmp0_1 = D_0*w10;
3359                                    EM_S[INDEX4(k,m,0,1,numEq,numComp,4)]+=tmp0_1;
3360                                    EM_S[INDEX4(k,m,1,2,numEq,numComp,4)]+=tmp0_1;
3361                                    EM_S[INDEX4(k,m,3,2,numEq,numComp,4)]+=tmp0_1;
3362                                    EM_S[INDEX4(k,m,0,0,numEq,numComp,4)]+=tmp0_1;
3363                                    EM_S[INDEX4(k,m,3,3,numEq,numComp,4)]+=tmp0_1;
3364                                    EM_S[INDEX4(k,m,3,0,numEq,numComp,4)]+=tmp0_1;
3365                                    EM_S[INDEX4(k,m,3,1,numEq,numComp,4)]+=tmp0_1;
3366                                    EM_S[INDEX4(k,m,2,1,numEq,numComp,4)]+=tmp0_1;
3367                                    EM_S[INDEX4(k,m,0,2,numEq,numComp,4)]+=tmp0_1;
3368                                    EM_S[INDEX4(k,m,2,0,numEq,numComp,4)]+=tmp0_1;
3369                                    EM_S[INDEX4(k,m,1,3,numEq,numComp,4)]+=tmp0_1;
3370                                    EM_S[INDEX4(k,m,2,3,numEq,numComp,4)]+=tmp0_1;
3371                                    EM_S[INDEX4(k,m,2,2,numEq,numComp,4)]+=tmp0_1;
3372                                    EM_S[INDEX4(k,m,1,0,numEq,numComp,4)]+=tmp0_1;
3373                                    EM_S[INDEX4(k,m,0,3,numEq,numComp,4)]+=tmp0_1;
3374                                    EM_S[INDEX4(k,m,1,1,numEq,numComp,4)]+=tmp0_1;
3375                                }
3376                            }
3377                        }
3378                        ///////////////
3379                        // process X //
3380                        ///////////////
3381                        if (!X.isEmpty()) {
3382                            add_EM_F=true;
3383                            const double* X_p=const_cast<escript::Data*>(&X)->getSampleDataRO(e);
3384                            for (index_t k=0; k<numEq; k++) {
3385                                const register double X_0 = X_p[INDEX2(k, 0, numEq)];
3386                                const register double X_1 = X_p[INDEX2(k, 1, numEq)];
3387                                const register double tmp0_1 = X_0*w11;
3388                                const register double tmp2_1 = X_0*w13;
3389                                const register double tmp1_1 = X_1*w12;
3390                                const register double tmp3_1 = X_1*w14;
3391                                EM_F[INDEX2(k,0,numEq)]+=tmp0_1 + tmp1_1;
3392                                EM_F[INDEX2(k,1,numEq)]+=tmp1_1 + tmp2_1;
3393                                EM_F[INDEX2(k,2,numEq)]+=tmp0_1 + tmp3_1;
3394                                EM_F[INDEX2(k,3,numEq)]+=tmp2_1 + tmp3_1;
3395                            }
3396                        }
3397                        ///////////////
3398                        // process Y //
3399                        ///////////////
3400                        if (!Y.isEmpty()) {
3401                            add_EM_F=true;
3402                            const double* Y_p=const_cast<escript::Data*>(&Y)->getSampleDataRO(e);
3403                            for (index_t k=0; k<numEq; k++) {
3404                                const register double Y_0 = Y_p[k];
3405                                const register double tmp0_1 = Y_0*w15;
3406                                EM_F[INDEX2(k,0,numEq)]+=tmp0_1;
3407                                EM_F[INDEX2(k,1,numEq)]+=tmp0_1;
3408                                EM_F[INDEX2(k,2,numEq)]+=tmp0_1;
3409                                EM_F[INDEX2(k,3,numEq)]+=tmp0_1;
3410                            }
3411                        }
3412                        /* GENERATOR SNIP_PDE_SYSTEM_REDUCED BOTTOM */
3413    
3414                        // add to matrix (if add_EM_S) and RHS (if add_EM_F)
3415                        const index_t firstNode=m_N0*k1+k0;
3416                      IndexVector rowIndex;                      IndexVector rowIndex;
3417                      rowIndex.push_back(m_dofMap[firstNode]);                      rowIndex.push_back(m_dofMap[firstNode]);
3418                      rowIndex.push_back(m_dofMap[firstNode+1]);                      rowIndex.push_back(m_dofMap[firstNode+1]);

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

  ViewVC Help
Powered by ViewVC 1.1.26