/[escript]/branches/domexper/dudley/src/Assemble_PDE_System2_3D.c
ViewVC logotype

Diff of /branches/domexper/dudley/src/Assemble_PDE_System2_3D.c

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

revision 3144 by jfenwick, Fri Sep 3 00:49:02 2010 UTC revision 3184 by jfenwick, Wed Sep 15 00:23:42 2010 UTC
# Line 53  void  Dudley_Assemble_PDE_System2_3D(Ass Line 53  void  Dudley_Assemble_PDE_System2_3D(Ass
53      index_t color;      index_t color;
54      dim_t e;      dim_t e;
55      __const double  *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p, *A_q, *B_q, *C_q, *D_q, *X_q, *Y_q;      __const double  *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p, *A_q, *B_q, *C_q, *D_q, *X_q, *Y_q;
56      double *EM_S, *EM_F, *Vol, *DSDX;      double *EM_S, *EM_F, *DSDX;
57      index_t *row_index;      index_t *row_index;
58      register dim_t q, s,r,k,m;      register dim_t q, s,r,k,m;
59      register double rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22;      register double rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22;
# Line 90  void  Dudley_Assemble_PDE_System2_3D(Ass Line 90  void  Dudley_Assemble_PDE_System2_3D(Ass
90                     D_p=getSampleDataRO(D,e);                     D_p=getSampleDataRO(D,e);
91                     X_p=getSampleDataRO(X,e);                     X_p=getSampleDataRO(X,e);
92                     Y_p=getSampleDataRO(Y,e);                     Y_p=getSampleDataRO(Y,e);
93              double vol=p.row_jac->absD[e]*p.row_jac->quadweight;
94    
                       Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);  
95                        DSDX=&(p.row_jac->DSDX[INDEX5(0,0,0,0,e,p.row_numShapesTotal,DIM,p.numQuadTotal,1)]);                        DSDX=&(p.row_jac->DSDX[INDEX5(0,0,0,0,e,p.row_numShapesTotal,DIM,p.numQuadTotal,1)]);
96                        for (q=0;q<len_EM_S;++q) EM_S[q]=0;                        for (q=0;q<len_EM_S;++q) EM_S[q]=0;
97                        for (q=0;q<len_EM_F;++q) EM_F[q]=0;                        for (q=0;q<len_EM_F;++q) EM_F[q]=0;
# Line 111  void  Dudley_Assemble_PDE_System2_3D(Ass Line 111  void  Dudley_Assemble_PDE_System2_3D(Ass
111                                    for (m=0;m<p.numComp;m++) {                                    for (m=0;m<p.numComp;m++) {
112                                       rtmp=0;                                       rtmp=0;
113                                       for (q=0;q<p.numQuadTotal;q++) {                                       for (q=0;q<p.numQuadTotal;q++) {
114                                          rtmp+=Vol[q]*(                                          rtmp+=vol*(
115                                            DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,0,q, p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]                                            DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,0,q, p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
116                                           +DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,1,q, p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]                                           +DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,1,q, p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]
117                                           +DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,2,q, p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]                                           +DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,2,q, p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]
# Line 141  void  Dudley_Assemble_PDE_System2_3D(Ass Line 141  void  Dudley_Assemble_PDE_System2_3D(Ass
141                                    rtmp21=0;                                    rtmp21=0;
142                                    rtmp22=0;                                    rtmp22=0;
143                                    for (q=0;q<p.numQuadTotal;q++) {                                    for (q=0;q<p.numQuadTotal;q++) {
144                                          rtmp0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                                          rtmp0=vol*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
145                                          rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                                          rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
146                                          rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                                          rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
147                                          rtmp02+=rtmp0*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];                                          rtmp02+=rtmp0*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];
148                
149                                          rtmp1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];                                          rtmp1=vol*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
150                                          rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                                          rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
151                                          rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                                          rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
152                                          rtmp12+=rtmp1*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];                                          rtmp12+=rtmp1*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];
153                
154                                          rtmp2=Vol[q]*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];                                          rtmp2=vol*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];
155                                          rtmp20+=rtmp2*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                                          rtmp20+=rtmp2*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
156                                          rtmp21+=rtmp2*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                                          rtmp21+=rtmp2*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
157                                          rtmp22+=rtmp2*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];                                          rtmp22+=rtmp2*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];
# Line 186  void  Dudley_Assemble_PDE_System2_3D(Ass Line 186  void  Dudley_Assemble_PDE_System2_3D(Ass
186                                    for (m=0;m<p.numComp;m++) {                                    for (m=0;m<p.numComp;m++) {
187                                       rtmp=0;                                       rtmp=0;
188                                       for (q=0;q<p.numQuadTotal;q++) {                                       for (q=0;q<p.numQuadTotal;q++) {
189                                           rtmp+=Vol[q]*S[INDEX2(r,q,p.row_numShapes)]*(                                           rtmp+=vol*S[INDEX2(r,q,p.row_numShapes)]*(
190                                                   DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]                                                   DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]
191                                                 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]                                                 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]
192                                                 + DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,2,m,q,p.numEqu,DIM,p.numComp)] );                                                 + DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,2,m,q,p.numEqu,DIM,p.numComp)] );
# Line 203  void  Dudley_Assemble_PDE_System2_3D(Ass Line 203  void  Dudley_Assemble_PDE_System2_3D(Ass
203                                       rtmp1=0;                                       rtmp1=0;
204                                       rtmp2=0;                                       rtmp2=0;
205                                       for (q=0;q<p.numQuadTotal;q++) {                                       for (q=0;q<p.numQuadTotal;q++) {
206                                            rtmp=Vol[q]*S[INDEX2(r,q,p.row_numShapes)];                                            rtmp=vol*S[INDEX2(r,q,p.row_numShapes)];
207                                            rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                                            rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
208                                            rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];                                            rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
209                                            rtmp2+=rtmp*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];                                            rtmp2+=rtmp*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];
# Line 232  void  Dudley_Assemble_PDE_System2_3D(Ass Line 232  void  Dudley_Assemble_PDE_System2_3D(Ass
232                                    for (m=0;m<p.numComp;m++) {                                    for (m=0;m<p.numComp;m++) {
233                                      rtmp=0;                                      rtmp=0;
234                                      for (q=0;q<p.numQuadTotal;q++) {                                      for (q=0;q<p.numQuadTotal;q++) {
235                                           rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*(                                           rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*(
236                                                      C_q[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]                                                      C_q[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
237                                                     +C_q[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]                                                     +C_q[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]
238                                                     +C_q[INDEX4(k,m,2,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]);                                                     +C_q[INDEX4(k,m,2,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]);
# Line 249  void  Dudley_Assemble_PDE_System2_3D(Ass Line 249  void  Dudley_Assemble_PDE_System2_3D(Ass
249                                       rtmp1=0;                                       rtmp1=0;
250                                       rtmp2=0;                                       rtmp2=0;
251                                       for (q=0;q<p.numQuadTotal;q++) {                                       for (q=0;q<p.numQuadTotal;q++) {
252                                          rtmp=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];                                          rtmp=vol*S[INDEX2(s,q,p.row_numShapes)];
253                                          rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                                          rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
254                                          rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                                          rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
255                                          rtmp2+=rtmp*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];                                          rtmp2+=rtmp*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];
# Line 278  void  Dudley_Assemble_PDE_System2_3D(Ass Line 278  void  Dudley_Assemble_PDE_System2_3D(Ass
278                                    for (m=0;m<p.numComp;m++) {                                    for (m=0;m<p.numComp;m++) {
279                                      rtmp=0;                                      rtmp=0;
280                                      for (q=0;q<p.numQuadTotal;q++) {                                      for (q=0;q<p.numQuadTotal;q++) {
281                                          rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_q[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_numShapes)];                                          rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*D_q[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_numShapes)];
282                                      }                                      }
283                                      EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                                      EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;
284                                    }                                    }
# Line 289  void  Dudley_Assemble_PDE_System2_3D(Ass Line 289  void  Dudley_Assemble_PDE_System2_3D(Ass
289                              for (s=0;s<p.row_numShapes;s++) {                              for (s=0;s<p.row_numShapes;s++) {
290                                for (r=0;r<p.col_numShapes;r++) {                                for (r=0;r<p.col_numShapes;r++) {
291                                    rtmp=0;                                    rtmp=0;
292                                    for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(r,q,p.row_numShapes)];                                    for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(r,q,p.row_numShapes)];
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_numShapesTotal)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];
# Line 310  void  Dudley_Assemble_PDE_System2_3D(Ass Line 310  void  Dudley_Assemble_PDE_System2_3D(Ass
310                                for (k=0;k<p.numEqu;k++) {                                for (k=0;k<p.numEqu;k++) {
311                                   rtmp=0;                                   rtmp=0;
312                                   for (q=0;q<p.numQuadTotal;q++) {                                   for (q=0;q<p.numQuadTotal;q++) {
313                                         rtmp+=Vol[q]* ( DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,0,q,p.numEqu,DIM)]                                         rtmp+=vol* ( DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,0,q,p.numEqu,DIM)]
314                                                       + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,1,q,p.numEqu,DIM)]                                                       + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,1,q,p.numEqu,DIM)]
315                                                       + DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,2,q,p.numEqu,DIM)]);                                                       + DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,2,q,p.numEqu,DIM)]);
316                                   }                                   }
# Line 323  void  Dudley_Assemble_PDE_System2_3D(Ass Line 323  void  Dudley_Assemble_PDE_System2_3D(Ass
323                               rtmp1=0;                               rtmp1=0;
324                               rtmp2=0;                               rtmp2=0;
325                               for (q=0;q<p.numQuadTotal;q++) {                               for (q=0;q<p.numQuadTotal;q++) {
326                                  rtmp0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                                  rtmp0+=vol*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
327                                  rtmp1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];                                  rtmp1+=vol*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
328                                  rtmp2+=Vol[q]*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];                                  rtmp2+=vol*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];
329                               }                               }
330                               for (k=0;k<p.numEqu;k++) {                               for (k=0;k<p.numEqu;k++) {
331                                         EM_F[INDEX2(k,s,p.numEqu)]+=rtmp0*X_p[INDEX2(k,0,p.numEqu)]                                         EM_F[INDEX2(k,s,p.numEqu)]+=rtmp0*X_p[INDEX2(k,0,p.numEqu)]
# Line 345  void  Dudley_Assemble_PDE_System2_3D(Ass Line 345  void  Dudley_Assemble_PDE_System2_3D(Ass
345                             for (s=0;s<p.row_numShapes;s++) {                             for (s=0;s<p.row_numShapes;s++) {
346                                for (k=0;k<p.numEqu;k++) {                                for (k=0;k<p.numEqu;k++) {
347                                   rtmp=0.;                                   rtmp=0.;
348                                   for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*Y_q[INDEX2(k,q,p.numEqu)];                                   for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*Y_q[INDEX2(k,q,p.numEqu)];
349                                   EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;                                   EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
350                                }                                }
351                             }                             }
352                           } else {                           } else {
353                             for (s=0;s<p.row_numShapes;s++) {                             for (s=0;s<p.row_numShapes;s++) {
354                                 rtmp=0;                                 rtmp=0;
355                                 for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];                                 for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)];
356                                 for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];                                 for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];
357                             }                             }
358                           }                           }

Legend:
Removed from v.3144  
changed lines
  Added in v.3184

  ViewVC Help
Powered by ViewVC 1.1.26