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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26