/[escript]/branches/doubleplusgood/dudley/src/Assemble_PDE_System2_3D.cpp
ViewVC logotype

Diff of /branches/doubleplusgood/dudley/src/Assemble_PDE_System2_3D.cpp

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

revision 852 by gross, Fri Aug 4 01:05:36 2006 UTC revision 853 by gross, Wed Sep 20 05:56:36 2006 UTC
# Line 41  Line 41 
41    
42  #include "Assemble.h"  #include "Assemble.h"
43  #include "Util.h"  #include "Util.h"
44    #ifdef _OPENMP
45    #include <omp.h>
46    #endif
47    
48    
49  /**************************************************************/  /**************************************************************/
50    
# Line 51  void  Finley_Assemble_PDE_System2_3D(Ass Line 55  void  Finley_Assemble_PDE_System2_3D(Ass
55      #define DIM 3      #define DIM 3
56      index_t color;      index_t color;
57      dim_t e;      dim_t e;
58        double *EM_S, *EM_F, *Vol, *DSDX, *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p;
59        index_t *row_index;
60        register dim_t q, s,r,k,m;
61        register double rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22;
62        bool_t add_EM_F, add_EM_S;
63    
64      bool_t extendedA=isExpanded(A);      bool_t extendedA=isExpanded(A);
65      bool_t extendedB=isExpanded(B);      bool_t extendedB=isExpanded(B);
66      bool_t extendedC=isExpanded(C);      bool_t extendedC=isExpanded(C);
# Line 63  void  Finley_Assemble_PDE_System2_3D(Ass Line 73  void  Finley_Assemble_PDE_System2_3D(Ass
73      dim_t len_EM_F=p.row_NN*p.numEqu;      dim_t len_EM_F=p.row_NN*p.numEqu;
74    
75    
76      #pragma omp parallel private(color)      #pragma omp parallel private(color,EM_S, EM_F, Vol, DSDX, A_p, B_p, C_p, D_p, X_p, Y_p,row_index,q, s,r,k,m,rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22,add_EM_F, add_EM_S)
77      {      {
78         double EM_S[len_EM_S], EM_F[len_EM_F];         EM_S=THREAD_MEMALLOC(len_EM_S,double);
79         index_t row_index[p.row_NN];         EM_F=THREAD_MEMALLOC(len_EM_F,double);
80         register dim_t q, s,r,k,m;         row_index=THREAD_MEMALLOC(p.row_NN,index_t);
81         register double rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22;                                                                                                                                                                                                      
82         double *Vol, *DSDX, *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p;         if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) {
        bool_t add_EM_F, add_EM_S;  
   
        #ifndef PASO_MPI  
        for (color=elements->minColor;color<=elements->maxColor;color++) {  
           /*  open loop over all elements: */  
           #pragma omp for private(e) schedule(static)  
           for(e=0;e<elements->numElements;e++){  
              if (elements->Color[e]==color) {  
        #else  
        {  
           for(e=0;e<elements->numElements;e++) {  
              {  
        #endif  
                 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);  
                 DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);  
                 for (q=0;q<len_EM_S;++q) EM_S[q]=0;  
                 for (q=0;q<len_EM_F;++q) EM_F[q]=0;  
                 add_EM_F=FALSE;  
                 add_EM_S=FALSE;  
                 /**************************************************************/  
                 /*   process A: */  
                 /**************************************************************/  
                 A_p=getSampleData(A,e);  
                 if (NULL!=A_p) {  
                    add_EM_S=TRUE;  
                    if (extendedA) {  
                       for (s=0;s<p.row_NS;s++) {  
                         for (r=0;r<p.col_NS;r++) {  
                           for (k=0;k<p.numEqu;k++) {  
                             for (m=0;m<p.numComp;m++) {  
                                rtmp=0;  
                                for (q=0;q<p.numQuad;q++) {  
                                   rtmp+=Vol[q]*(  
                                     DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX5(k,0,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]  
                                    +DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX5(k,0,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]  
                                    +DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX5(k,0,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]  
                                    +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX5(k,1,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]  
                                    +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX5(k,1,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]  
                                    +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX5(k,1,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]  
                                    +DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*A_p[INDEX5(k,2,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]  
                                    +DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*A_p[INDEX5(k,2,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]  
                                    +DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*A_p[INDEX5(k,2,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]);  
                         
                                }  
                                EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;  
                             }  
                           }  
                         }  
                       }  
                    } else {  
                       for (s=0;s<p.row_NS;s++) {  
                         for (r=0;r<p.col_NS;r++) {  
                             rtmp00=0;  
                             rtmp01=0;  
                             rtmp02=0;  
                             rtmp10=0;  
                             rtmp11=0;  
                             rtmp12=0;  
                             rtmp20=0;  
                             rtmp21=0;  
                             rtmp22=0;  
                             for (q=0;q<p.numQuad;q++) {  
                                   rtmp0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];  
                                   rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];  
                                   rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];  
                                   rtmp02+=rtmp0*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];  
   
                                   rtmp1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];  
                                   rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];  
                                   rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];  
                                   rtmp12+=rtmp1*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];  
83    
84                                    rtmp2=Vol[q]*DSDX[INDEX3(s,2,q,p.row_NN,DIM)];            #ifndef PASO_MPI
85                                    rtmp20+=rtmp2*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];            for (color=elements->minColor;color<=elements->maxColor;color++) {
86                                    rtmp21+=rtmp2*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];               /*  open loop over all elements: */
87                                    rtmp22+=rtmp2*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];               #pragma omp for private(e) schedule(static)
88                              }               for(e=0;e<elements->numElements;e++){
89                              for (k=0;k<p.numEqu;k++) {                  if (elements->Color[e]==color) {
90              #else
91              {
92                 for(e=0;e<elements->numElements;e++) {
93                    {
94              #endif
95                       Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
96                       DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);
97                       for (q=0;q<len_EM_S;++q) EM_S[q]=0;
98                       for (q=0;q<len_EM_F;++q) EM_F[q]=0;
99                       add_EM_F=FALSE;
100                       add_EM_S=FALSE;
101                       /**************************************************************/
102                       /*   process A: */
103                       /**************************************************************/
104                       A_p=getSampleData(A,e);
105                       if (NULL!=A_p) {
106                          add_EM_S=TRUE;
107                          if (extendedA) {
108                             for (s=0;s<p.row_NS;s++) {
109                               for (r=0;r<p.col_NS;r++) {
110                                 for (k=0;k<p.numEqu;k++) {
111                                 for (m=0;m<p.numComp;m++) {                                 for (m=0;m<p.numComp;m++) {
112                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+= rtmp00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]                                    rtmp=0;
113                                                                                       +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]                                    for (q=0;q<p.numQuad;q++) {
114                                                                                       +rtmp02*A_p[INDEX4(k,0,m,2,p.numEqu,DIM,p.numComp)]                                       rtmp+=Vol[q]*(
115                                                                                       +rtmp10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]                                         DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX5(k,0,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
116                                                                                       +rtmp11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)]                                        +DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX5(k,0,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]
117                                                                                       +rtmp12*A_p[INDEX4(k,1,m,2,p.numEqu,DIM,p.numComp)]                                        +DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX5(k,0,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]
118                                                                                       +rtmp20*A_p[INDEX4(k,2,m,0,p.numEqu,DIM,p.numComp)]                                        +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX5(k,1,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
119                                                                                       +rtmp21*A_p[INDEX4(k,2,m,1,p.numEqu,DIM,p.numComp)]                                        +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX5(k,1,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]
120                                                                                       +rtmp22*A_p[INDEX4(k,2,m,2,p.numEqu,DIM,p.numComp)];                                        +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX5(k,1,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]
121                                 }                                        +DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*A_p[INDEX5(k,2,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
122                              }                                        +DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*A_p[INDEX5(k,2,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]
123                          }                                        +DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*A_p[INDEX5(k,2,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]);
124                        }                            
125                      }                                    }
126                  }                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
                 /**************************************************************/  
                 /*   process B: */  
                 /**************************************************************/  
                 B_p=getSampleData(B,e);  
                 if (NULL!=B_p) {  
                    add_EM_S=TRUE;  
                    if (extendedB) {  
                       for (s=0;s<p.row_NS;s++) {  
                         for (r=0;r<p.col_NS;r++) {  
                           for (k=0;k<p.numEqu;k++) {  
                             for (m=0;m<p.numComp;m++) {  
                                rtmp=0;  
                                for (q=0;q<p.numQuad;q++) {  
                                    rtmp+=Vol[q]*S[INDEX2(r,q,p.row_NS)]*(  
                                            DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*B_p[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]  
                                          + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*B_p[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]  
                                          + DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*B_p[INDEX4(k,2,m,q,p.numEqu,DIM,p.numComp)] );  
127                                 }                                 }
128                                 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                               }
129                              }                             }
130                            }                           }
131                          }                        } else {
132                        }                           for (s=0;s<p.row_NS;s++) {
133                     } else {                             for (r=0;r<p.col_NS;r++) {
134                        for (s=0;s<p.row_NS;s++) {                                 rtmp00=0;
135                          for (r=0;r<p.col_NS;r++) {                                 rtmp01=0;
136                                 rtmp0=0;                                 rtmp02=0;
137                                 rtmp1=0;                                 rtmp10=0;
138                                 rtmp2=0;                                 rtmp11=0;
139                                   rtmp12=0;
140                                   rtmp20=0;
141                                   rtmp21=0;
142                                   rtmp22=0;
143                                 for (q=0;q<p.numQuad;q++) {                                 for (q=0;q<p.numQuad;q++) {
144                                      rtmp=Vol[q]*S[INDEX2(r,q,p.row_NS)];                                       rtmp0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
145                                      rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                                       rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
146                                      rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                                       rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
147                                      rtmp2+=rtmp*DSDX[INDEX3(s,2,q,p.row_NN,DIM)];                                       rtmp02+=rtmp0*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];
148      
149                                         rtmp1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
150                                         rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
151                                         rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
152                                         rtmp12+=rtmp1*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];
153      
154                                         rtmp2=Vol[q]*DSDX[INDEX3(s,2,q,p.row_NN,DIM)];
155                                         rtmp20+=rtmp2*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
156                                         rtmp21+=rtmp2*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
157                                         rtmp22+=rtmp2*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];
158                                 }                                 }
159                                 for (k=0;k<p.numEqu;k++) {                                 for (k=0;k<p.numEqu;k++) {
160                                    for (m=0;m<p.numComp;m++) {                                    for (m=0;m<p.numComp;m++) {
161                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+= rtmp00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]
162                                                                                         +rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)]                                                                                          +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]
163                                                                                         +rtmp2*B_p[INDEX3(k,2,m,p.numEqu,DIM)];                                                                                          +rtmp02*A_p[INDEX4(k,0,m,2,p.numEqu,DIM,p.numComp)]
164                                                                                            +rtmp10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]
165                                                                                            +rtmp11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)]
166                                                                                            +rtmp12*A_p[INDEX4(k,1,m,2,p.numEqu,DIM,p.numComp)]
167                                                                                            +rtmp20*A_p[INDEX4(k,2,m,0,p.numEqu,DIM,p.numComp)]
168                                                                                            +rtmp21*A_p[INDEX4(k,2,m,1,p.numEqu,DIM,p.numComp)]
169                                                                                            +rtmp22*A_p[INDEX4(k,2,m,2,p.numEqu,DIM,p.numComp)];
170                                    }                                    }
171                                 }                                 }
172                          }                             }
173                        }                           }
174                           }
175                     }                     }
176                  }                     /**************************************************************/
177                  /**************************************************************/                     /*   process B: */
178                  /*   process C: */                     /**************************************************************/
179                  /**************************************************************/                     B_p=getSampleData(B,e);
180                  C_p=getSampleData(C,e);                     if (NULL!=B_p) {
181                  if (NULL!=C_p) {                        add_EM_S=TRUE;
182                    add_EM_S=TRUE;                        if (extendedB) {
183                    if (extendedC) {                           for (s=0;s<p.row_NS;s++) {
184                        for (s=0;s<p.row_NS;s++) {                             for (r=0;r<p.col_NS;r++) {
185                          for (r=0;r<p.col_NS;r++) {                               for (k=0;k<p.numEqu;k++) {
186                            for (k=0;k<p.numEqu;k++) {                                 for (m=0;m<p.numComp;m++) {
187                              for (m=0;m<p.numComp;m++) {                                    rtmp=0;
188                                rtmp=0;                                    for (q=0;q<p.numQuad;q++) {
189                                for (q=0;q<p.numQuad;q++) {                                        rtmp+=Vol[q]*S[INDEX2(r,q,p.row_NS)]*(
190                                     rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*(                                                DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*B_p[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]
191                                                C_p[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]                                              + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*B_p[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]
192                                               +C_p[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]                                              + DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*B_p[INDEX4(k,2,m,q,p.numEqu,DIM,p.numComp)] );
193                                               +C_p[INDEX4(k,m,2,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]);                                    }
194                                }                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
195                                EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                                 }
196                              }                               }
197                            }                             }
198                          }                           }
199                          } else {
200                             for (s=0;s<p.row_NS;s++) {
201                               for (r=0;r<p.col_NS;r++) {
202                                      rtmp0=0;
203                                      rtmp1=0;
204                                      rtmp2=0;
205                                      for (q=0;q<p.numQuad;q++) {
206                                           rtmp=Vol[q]*S[INDEX2(r,q,p.row_NS)];
207                                           rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
208                                           rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
209                                           rtmp2+=rtmp*DSDX[INDEX3(s,2,q,p.row_NN,DIM)];
210                                      }
211                                      for (k=0;k<p.numEqu;k++) {
212                                         for (m=0;m<p.numComp;m++) {
213                                            EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]
214                                                                                              +rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)]
215                                                                                              +rtmp2*B_p[INDEX3(k,2,m,p.numEqu,DIM)];
216                                         }
217                                      }
218                               }
219                             }
220                        }                        }
221                    } else {                     }
222                        for (s=0;s<p.row_NS;s++) {                     /**************************************************************/
223                          for (r=0;r<p.col_NS;r++) {                     /*   process C: */
224                                 rtmp0=0;                     /**************************************************************/
225                                 rtmp1=0;                     C_p=getSampleData(C,e);
226                                 rtmp2=0;                     if (NULL!=C_p) {
227                                 for (q=0;q<p.numQuad;q++) {                       add_EM_S=TRUE;
228                                    rtmp=Vol[q]*S[INDEX2(s,q,p.row_NS)];                       if (extendedC) {
229                                    rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                           for (s=0;s<p.row_NS;s++) {
230                                    rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                             for (r=0;r<p.col_NS;r++) {
231                                    rtmp2+=rtmp*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];                               for (k=0;k<p.numEqu;k++) {
232                                   for (m=0;m<p.numComp;m++) {
233                                     rtmp=0;
234                                     for (q=0;q<p.numQuad;q++) {
235                                          rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*(
236                                                     C_p[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
237                                                    +C_p[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]
238                                                    +C_p[INDEX4(k,m,2,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]);
239                                     }
240                                     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
241                                 }                                 }
242                                 }
243                               }
244                             }
245                         } else {
246                             for (s=0;s<p.row_NS;s++) {
247                               for (r=0;r<p.col_NS;r++) {
248                                      rtmp0=0;
249                                      rtmp1=0;
250                                      rtmp2=0;
251                                      for (q=0;q<p.numQuad;q++) {
252                                         rtmp=Vol[q]*S[INDEX2(s,q,p.row_NS)];
253                                         rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
254                                         rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
255                                         rtmp2+=rtmp*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];
256                                      }
257                                      for (k=0;k<p.numEqu;k++) {
258                                         for (m=0;m<p.numComp;m++) {
259                                               EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp0*C_p[INDEX3(k,m,0,p.numEqu,p.numComp)]
260                                                                                                 +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)]
261                                                                                                 +rtmp2*C_p[INDEX3(k,m,2,p.numEqu,p.numComp)];
262                                          }
263                                      }
264                               }
265                             }
266                         }
267                       }
268                       /************************************************************* */
269                       /* process D */
270                       /**************************************************************/
271                       D_p=getSampleData(D,e);
272                       if (NULL!=D_p) {
273                         add_EM_S=TRUE;
274                         if (extendedD) {
275                             for (s=0;s<p.row_NS;s++) {
276                               for (r=0;r<p.col_NS;r++) {
277                                 for (k=0;k<p.numEqu;k++) {
278                                   for (m=0;m<p.numComp;m++) {
279                                     rtmp=0;
280                                     for (q=0;q<p.numQuad;q++) {
281                                         rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_NS)];
282                                     }
283                                     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
284                                   }
285                                 }
286                               }
287                             }
288                         } else {
289                             for (s=0;s<p.row_NS;s++) {
290                               for (r=0;r<p.col_NS;r++) {
291                                   rtmp=0;
292                                   for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];
293                                 for (k=0;k<p.numEqu;k++) {                                 for (k=0;k<p.numEqu;k++) {
294                                    for (m=0;m<p.numComp;m++) {                                     for (m=0;m<p.numComp;m++) {
295                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp0*C_p[INDEX3(k,m,0,p.numEqu,p.numComp)]                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];
296                                                                                            +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)]                                    }
                                                                                           +rtmp2*C_p[INDEX3(k,m,2,p.numEqu,p.numComp)];  
                                    }  
297                                 }                                 }
298                          }                             }
299                        }                           }
300                    }                       }
301                  }                     }
302                  /************************************************************* */                     /**************************************************************/
303                  /* process D */                     /*   process X: */
304                  /**************************************************************/                     /**************************************************************/
305                  D_p=getSampleData(D,e);                     X_p=getSampleData(X,e);
306                  if (NULL!=D_p) {                     if (NULL!=X_p) {
307                    add_EM_S=TRUE;                       add_EM_F=TRUE;
308                    if (extendedD) {                       if (extendedX) {
309                        for (s=0;s<p.row_NS;s++) {                          for (s=0;s<p.row_NS;s++) {
310                          for (r=0;r<p.col_NS;r++) {                             for (k=0;k<p.numEqu;k++) {
                           for (k=0;k<p.numEqu;k++) {  
                             for (m=0;m<p.numComp;m++) {  
311                                rtmp=0;                                rtmp=0;
312                                for (q=0;q<p.numQuad;q++) {                                for (q=0;q<p.numQuad;q++) {
313                                    rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_NS)];                                      rtmp+=Vol[q]* ( DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX3(k,0,q,p.numEqu,DIM)]
314                                                      + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*X_p[INDEX3(k,1,q,p.numEqu,DIM)]
315                                                      + DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*X_p[INDEX3(k,2,q,p.numEqu,DIM)]);
316                                }                                }
317                                EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                                EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
318                              }                             }
                           }  
319                          }                          }
320                        }                       } else {
321                    } else {                          for (s=0;s<p.row_NS;s++) {
322                        for (s=0;s<p.row_NS;s++) {                            rtmp0=0;
323                          for (r=0;r<p.col_NS;r++) {                            rtmp1=0;
324                              rtmp=0;                            rtmp2=0;
325                              for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];                            for (q=0;q<p.numQuad;q++) {
326                              for (k=0;k<p.numEqu;k++) {                               rtmp0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
327                                  for (m=0;m<p.numComp;m++) {                               rtmp1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
328                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];                               rtmp2+=Vol[q]*DSDX[INDEX3(s,2,q,p.row_NN,DIM)];
329                                 }                            }
330                              }                            for (k=0;k<p.numEqu;k++) {
331                                        EM_F[INDEX2(k,s,p.numEqu)]+=rtmp0*X_p[INDEX2(k,0,p.numEqu)]
332                                                                   +rtmp1*X_p[INDEX2(k,1,p.numEqu)]
333                                                                   +rtmp2*X_p[INDEX2(k,2,p.numEqu)];
334                              }
335                          }                          }
336                        }                       }
337                    }                    }
338                  }                    /**************************************************************/
339                  /**************************************************************/                    /*   process Y: */
340                  /*   process X: */                    /**************************************************************/
341                  /**************************************************************/                     Y_p=getSampleData(Y,e);
342                  X_p=getSampleData(X,e);                     if (NULL!=Y_p) {
343                  if (NULL!=X_p) {                       add_EM_F=TRUE;
344                    add_EM_F=TRUE;                       if (extendedY) {
345                    if (extendedX) {                          for (s=0;s<p.row_NS;s++) {
346                       for (s=0;s<p.row_NS;s++) {                             for (k=0;k<p.numEqu;k++) {
347                          for (k=0;k<p.numEqu;k++) {                                rtmp=0.;
348                             rtmp=0;                                for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[INDEX2(k,q,p.numEqu)];
349                             for (q=0;q<p.numQuad;q++) {                                EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
                                  rtmp+=Vol[q]* ( DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX3(k,0,q,p.numEqu,DIM)]  
                                                + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*X_p[INDEX3(k,1,q,p.numEqu,DIM)]  
                                                + DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*X_p[INDEX3(k,2,q,p.numEqu,DIM)]);  
350                             }                             }
                            EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;  
351                          }                          }
352                       }                        } else {
353                    } else {                          for (s=0;s<p.row_NS;s++) {
354                       for (s=0;s<p.row_NS;s++) {                              rtmp=0;
355                         rtmp0=0;                              for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
356                         rtmp1=0;                              for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];
                        rtmp2=0;  
                        for (q=0;q<p.numQuad;q++) {  
                           rtmp0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];  
                           rtmp1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];  
                           rtmp2+=Vol[q]*DSDX[INDEX3(s,2,q,p.row_NN,DIM)];  
                        }  
                        for (k=0;k<p.numEqu;k++) {  
                                  EM_F[INDEX2(k,s,p.numEqu)]+=rtmp0*X_p[INDEX2(k,0,p.numEqu)]  
                                                             +rtmp1*X_p[INDEX2(k,1,p.numEqu)]  
                                                             +rtmp2*X_p[INDEX2(k,2,p.numEqu)];  
                        }  
                      }  
                   }  
                }  
                /**************************************************************/  
                /*   process Y: */  
                /**************************************************************/  
                 Y_p=getSampleData(Y,e);  
                 if (NULL!=Y_p) {  
                   add_EM_F=TRUE;  
                   if (extendedY) {  
                      for (s=0;s<p.row_NS;s++) {  
                         for (k=0;k<p.numEqu;k++) {  
                            rtmp=0.;  
                            for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[INDEX2(k,q,p.numEqu)];  
                            EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;  
357                          }                          }
358                       }                        }
359                     } else {                      }
360                       for (s=0;s<p.row_NS;s++) {                      /***********************************************************************************************/
361                           rtmp=0;                      /* add the element matrices onto the matrix and right hand side                                */
362                           for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];                      /***********************************************************************************************/
363                           for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];                      for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
364                       }                      if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);
365                     }                      if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);
366                   }    
367                   /***********************************************************************************************/                  } /* end color check */
368                   /* add the element matrices onto the matrix and right hand side                                */               } /* end element loop */
369                   /***********************************************************************************************/           } /* end color loop */
370                   for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];            
371                   if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);           THREAD_MEMFREE(EM_S);
372                   if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);           THREAD_MEMFREE(EM_F);
373             THREAD_MEMFREE(row_index);
374    
375               } /* end color check */        } /* end of pointer check */
           } /* end element loop */  
       } /* end color loop */  
376     } /* end parallel region */     } /* end parallel region */
377  }  }
378  /*  /*

Legend:
Removed from v.852  
changed lines
  Added in v.853

  ViewVC Help
Powered by ViewVC 1.1.26