/[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 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_2D(Ass Line 52  void  Finley_Assemble_PDE_System2_2D(Ass
52      #define DIM 2      #define DIM 2
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, rtmp00, rtmp10, rtmp01, rtmp11;
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_2D(Ass Line 70  void  Finley_Assemble_PDE_System2_2D(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, rtmp00, rtmp10, rtmp01, rtmp11,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, rtmp00, rtmp10, rtmp01, rtmp11;                                                                                                                                                                                                      
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++) {
102                  A_p=getSampleData(A,e);                                 for (m=0;m<p.numComp;m++) {
103                  if (NULL!=A_p) {                                   rtmp=0;
104                     add_EM_S=TRUE;                                   for (q=0;q<p.numQuad;q++) {
105                     if (extendedA) {                                      rtmp+=Vol[q]* (
106                        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)]
107                          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)]
108                            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)]
109                              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)]);
110                                rtmp=0;                                   }
111                                     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
112                                   }
113                                 }
114                               }
115                             }
116                          } else {
117                             for (s=0;s<p.row_NS;s++) {
118                               for (r=0;r<p.col_NS;r++) {
119                                  rtmp00=0;
120                                  rtmp01=0;
121                                  rtmp10=0;
122                                  rtmp11=0;
123                                for (q=0;q<p.numQuad;q++) {                                for (q=0;q<p.numQuad;q++) {
124                                   rtmp+=Vol[q]* (                                      rtmp0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
125                                       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)];
126                                      +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)];
127                                      +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)];
128                                      +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)];
129                                        rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
130                                  }
131                                  for (k=0;k<p.numEqu;k++) {
132                                      for (m=0;m<p.numComp;m++) {
133                                         EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=
134                                                    rtmp00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]
135                                                   +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]
136                                                   +rtmp10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]
137                                                   +rtmp11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)];
138                                      }
139                                }                                }
                               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)];  
140                             }                             }
141                             for (k=0;k<p.numEqu;k++) {                           }
142                         }
143                       }
144                       /**************************************************************/
145                       /*   process B: */
146                       /**************************************************************/
147                       B_p=getSampleData(B,e);
148                       if (NULL!=B_p) {
149                          add_EM_S=TRUE;
150                          if (extendedB) {
151                             for (s=0;s<p.row_NS;s++) {
152                               for (r=0;r<p.col_NS;r++) {
153                                 for (k=0;k<p.numEqu;k++) {
154                                 for (m=0;m<p.numComp;m++) {                                 for (m=0;m<p.numComp;m++) {
155                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=                                    rtmp=0;
156                                               rtmp00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]                                    for (q=0;q<p.numQuad;q++) {
157                                              +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]                                       rtmp+=Vol[q]*S[INDEX2(r,q,p.row_NS)]*
158                                              +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)]
159                                              +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)]);
160                                      }
161                                      EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
162                                 }                                 }
163                                 }
164                             }                             }
165                          }                           }
166                        }                        } else {
167                    }                           for (s=0;s<p.row_NS;s++) {
168                  }                             for (r=0;r<p.col_NS;r++) {
169                  /**************************************************************/                                 rtmp0=0;
170                  /*   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;  
171                                 for (q=0;q<p.numQuad;q++) {                                 for (q=0;q<p.numQuad;q++) {
172                                    rtmp+=Vol[q]*S[INDEX2(r,q,p.row_NS)]*                                     rtmp=Vol[q]*S[INDEX2(r,q,p.row_NS)];
173                                         ( 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)];
174                                         + 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)];
175                                 }                                 }
176                                 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                                 for (k=0;k<p.numEqu;k++) {
177                              }                                    for (m=0;m<p.numComp;m++) {
178                            }                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+= rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]
179                          }                                                                                         + rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)];
180                                      }
181                                   }
182                               }
183                             }
184                        }                        }
185                     } else {                     }
186                        for (s=0;s<p.row_NS;s++) {                     /**************************************************************/
187                          for (r=0;r<p.col_NS;r++) {                     /*   process C: */
188                              rtmp0=0;                     /**************************************************************/
189                              rtmp1=0;                     C_p=getSampleData(C,e);
190                              for (q=0;q<p.numQuad;q++) {                     if (NULL!=C_p) {
191                                  rtmp=Vol[q]*S[INDEX2(r,q,p.row_NS)];                       add_EM_S=TRUE;
192                                  rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                       if (extendedC) {
193                                  rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                           for (s=0;s<p.row_NS;s++) {
194                              }                             for (r=0;r<p.col_NS;r++) {
195                              for (k=0;k<p.numEqu;k++) {                               for (k=0;k<p.numEqu;k++) {
196                                 for (m=0;m<p.numComp;m++) {                                 for (m=0;m<p.numComp;m++) {
197                                    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;
198                                                                                      + rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)];                                    for (q=0;q<p.numQuad;q++) {
199                                          rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*
200                                               ( C_p[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
201                                               + C_p[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]);
202                                      }
203                                      EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
204                                 }                                 }
205                              }                               }
206                          }                             }
207                        }                           }
208                         } else {
209                             for (s=0;s<p.row_NS;s++) {
210                               for (r=0;r<p.col_NS;r++) {
211                                  rtmp0=0;
212                                  rtmp1=0;
213                                  for (q=0;q<p.numQuad;q++) {
214                                          rtmp=Vol[q]*S[INDEX2(s,q,p.row_NS)];
215                                          rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
216                                          rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
217                                  }
218                                  for (k=0;k<p.numEqu;k++) {
219                                     for (m=0;m<p.numComp;m++) {
220                                            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)]
221                                                                                              +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)];
222                                     }
223                                  }
224                               }
225                             }
226                         }
227                     }                     }
228                  }                     /************************************************************* */
229                  /**************************************************************/                     /* process D */
230                  /*   process C: */                     /**************************************************************/
231                  /**************************************************************/                     D_p=getSampleData(D,e);
232                  C_p=getSampleData(C,e);                     if (NULL!=D_p) {
233                  if (NULL!=C_p) {                       add_EM_S=TRUE;
234                    add_EM_S=TRUE;                       if (extendedD) {
235                    if (extendedC) {                           for (s=0;s<p.row_NS;s++) {
236                        for (s=0;s<p.row_NS;s++) {                             for (r=0;r<p.col_NS;r++) {
237                          for (r=0;r<p.col_NS;r++) {                               for (k=0;k<p.numEqu;k++) {
238                            for (k=0;k<p.numEqu;k++) {                                 for (m=0;m<p.numComp;m++) {
239                              for (m=0;m<p.numComp;m++) {                                   rtmp=0;
240                                     for (q=0;q<p.numQuad;q++) {
241                                         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)];
242                                     }
243                                     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
244                                   }
245                                 }
246                               }
247                             }
248                         } else {
249                             for (s=0;s<p.row_NS;s++) {
250                               for (r=0;r<p.col_NS;r++) {
251                                 rtmp=0;                                 rtmp=0;
252                                 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)];
253                                     rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*                                 for (k=0;k<p.numEqu;k++) {
254                                          ( 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++) {
255                                          + 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)];
256                                      }
257                                 }                                 }
                                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)];  
258                             }                             }
259                             }
260                         }
261                       }
262                       /**************************************************************/
263                       /*   process X: */
264                       /**************************************************************/
265                       X_p=getSampleData(X,e);
266                       if (NULL!=X_p) {
267                         add_EM_F=TRUE;
268                         if (extendedX) {
269                            for (s=0;s<p.row_NS;s++) {
270                             for (k=0;k<p.numEqu;k++) {                             for (k=0;k<p.numEqu;k++) {
271                                for (m=0;m<p.numComp;m++) {                               rtmp=0;
272                                       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++) {
273                                                                                         +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)]
274                                }                                                 +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*X_p[INDEX3(k,1,q,p.numEqu,DIM)]);
275                                 }
276                                 EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
277                             }                             }
278                          }                          }
279                        }                       } else {
280                    }                          for (s=0;s<p.row_NS;s++) {
281                  }                                rtmp0=0;
282                  /************************************************************* */                                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;  
283                                for (q=0;q<p.numQuad;q++) {                                for (q=0;q<p.numQuad;q++) {
284                                    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)];
285                                      rtmp1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
286                                }                                }
287                                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;  
288                          }                          }
289                       }                       }
                   } 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)];  
                      }  
290                    }                    }
291                 }                    /**************************************************************/
292                 /**************************************************************/                    /*   process Y: */
293                 /*   process Y: */                    /**************************************************************/
294                 /**************************************************************/                     Y_p=getSampleData(Y,e);
295                  Y_p=getSampleData(Y,e);                     if (NULL!=Y_p) {
296                  if (NULL!=Y_p) {                       add_EM_F=TRUE;
297                    add_EM_F=TRUE;                       if (extendedY) {
298                    if (extendedY) {                          for (s=0;s<p.row_NS;s++) {
299                       for (s=0;s<p.row_NS;s++) {                             for (k=0;k<p.numEqu;k++) {
300                          for (k=0;k<p.numEqu;k++) {                                rtmp=0;
301                             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)];
302                             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;
303                             EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;                             }
304                          }                          }
305                       }                        } else {
306                     } else {                          for (s=0;s<p.row_NS;s++) {
307                       for (s=0;s<p.row_NS;s++) {                              rtmp=0;
308                           rtmp=0;                              for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
309                           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];
310                           for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];                          }
311                       }                        }
312                     }                      }
313                   }                      /***********************************************************************************************/
314                   /***********************************************************************************************/                      /* add the element matrices onto the matrix and right hand side                                */
315                   /* add the element matrices onto the matrix and right hand side                                */                      /***********************************************************************************************/
316                   /***********************************************************************************************/                      for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
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)]];                      if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);
318                   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);
319                   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                    } /* end color check */
321               } /* end color check */               } /* end element loop */
322            } /* end element loop */           } /* end color loop */
323        } /* end color loop */            
324             THREAD_MEMFREE(EM_S);
325             THREAD_MEMFREE(EM_F);
326             THREAD_MEMFREE(row_index);
327    
328          } /* end of pointer check */
329     } /* end parallel region */     } /* end parallel region */
330  }  }
331  /*  /*

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

  ViewVC Help
Powered by ViewVC 1.1.26