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

Legend:
Removed from v.798  
changed lines
  Added in v.1811

  ViewVC Help
Powered by ViewVC 1.1.26