/[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

trunk/finley/src/Assemble_PDE_System2_2D.c revision 798 by gross, Fri Aug 4 01:05:36 2006 UTC branches/domexper/dudley/src/Assemble_PDE_System2_2D.c revision 3144 by jfenwick, Fri Sep 3 00:49:02 2010 UTC
# Line 1  Line 1 
1  /*  
2   ************************************************************  /*******************************************************
3   *          Copyright 2006 by ACcESS MNRF                   *  *
4   *                                                          *  * Copyright (c) 2003-2010 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    
48  void  Finley_Assemble_PDE_System2_2D(Assemble_Parameters p, Finley_ElementFile* elements,  void  Dudley_Assemble_PDE_System2_2D(Assemble_Parameters p, Dudley_ElementFile* elements,
49                                       Paso_SystemMatrix* Mat, escriptDataC* F,                                       Paso_SystemMatrix* Mat, escriptDataC* F,
50                                       escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y) {                                       escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y) {
51    
52      #define DIM 2      #define DIM 2
53      index_t color;      index_t color;
54      dim_t e;      dim_t e;
55        __const double *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p, *A_q, *B_q, *C_q, *D_q, *X_q, *Y_q;
56        double *EM_S, *EM_F, *Vol, *DSDX;
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);
65      bool_t extendedD=isExpanded(D);      bool_t extendedD=isExpanded(D);
66      bool_t extendedX=isExpanded(X);      bool_t extendedX=isExpanded(X);
67      bool_t extendedY=isExpanded(Y);      bool_t extendedY=isExpanded(Y);
68      double *F_p=getSampleData(F,0);      double *F_p=(requireWrite(F), getSampleDataRW(F,0));    /* use comma, to get around the mixed code and declarations thing */
69      double *S=p.row_jac->ReferenceElement->S;      double *S=p.row_jac->BasisFunctions->S;
70      dim_t len_EM_S=p.row_NN*p.col_NN*p.numEqu*p.numComp;      dim_t len_EM_S=p.row_numShapesTotal*p.col_numShapesTotal*p.numEqu*p.numComp;
71      dim_t len_EM_F=p.row_NN*p.numEqu;      dim_t len_EM_F=p.row_numShapesTotal*p.numEqu;
   
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, A_q, B_q, C_q, D_q, X_q, Y_q,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];  
76         index_t row_index[p.row_NN];         EM_S=THREAD_MEMALLOC(len_EM_S,double);
77         register dim_t q, s,r,k,m;         EM_F=THREAD_MEMALLOC(len_EM_F,double);
78         register double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11;         row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);
79         double *Vol, *DSDX, *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p;                                                                                                                                                                                                      
80         bool_t add_EM_F, add_EM_S;         if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index) ) {
81    
82         #ifndef PASO_MPI            for (color=elements->minColor;color<=elements->maxColor;color++) {
83         for (color=elements->minColor;color<=elements->maxColor;color++) {               /*  open loop over all elements: */
84            /*  open loop over all elements: */               #pragma omp for private(e) schedule(static)
85            #pragma omp for private(e) schedule(static)               for(e=0;e<elements->numElements;e++){
86            for(e=0;e<elements->numElements;e++){                  if (elements->Color[e]==color) {
87               if (elements->Color[e]==color) {                        
88         #else             A_p=getSampleDataRO(A,e);
89         {             B_p=getSampleDataRO(B,e);
90            for(e=0;e<elements->numElements;e++) {                     C_p=getSampleDataRO(C,e);    
91               {                     D_p=getSampleDataRO(D,e);
92         #endif                     X_p=getSampleDataRO(X,e);
93                  Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);                     Y_p=getSampleDataRO(Y,e);
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;                        Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);
96                  for (q=0;q<len_EM_F;++q) EM_F[q]=0;                        DSDX=&(p.row_jac->DSDX[INDEX5(0,0,0,0,e,p.row_numShapesTotal,DIM,p.numQuadTotal,1)]);
97                  add_EM_F=FALSE;                        for (q=0;q<len_EM_S;++q) EM_S[q]=0;
98                  add_EM_S=FALSE;                        for (q=0;q<len_EM_F;++q) EM_F[q]=0;
99                  /**************************************************************/                        add_EM_F=FALSE;
100                  /*   process A: */                        add_EM_S=FALSE;
101                  /**************************************************************/  
102                  A_p=getSampleData(A,e);                        /**************************************************************/
103                  if (NULL!=A_p) {                        /*   process A: */
104                     add_EM_S=TRUE;                        /**************************************************************/
105                     if (extendedA) {                        A_p=getSampleDataRO(A,e);
106                        for (s=0;s<p.row_NS;s++) {                        if (NULL!=A_p) {
107                          for (r=0;r<p.col_NS;r++) {                           add_EM_S=TRUE;
108                            for (k=0;k<p.numEqu;k++) {                           if (extendedA) {
109                              for (m=0;m<p.numComp;m++) {                  A_q=&(A_p[INDEX6(0,0,0,0,0,0, p.numEqu,DIM,p.numComp,DIM,p.numQuadTotal)]);
110                                rtmp=0;                              for (s=0;s<p.row_numShapes;s++) {
111                                for (q=0;q<p.numQuad;q++) {                                for (r=0;r<p.col_numShapes;r++) {
112                                   rtmp+=Vol[q]* (                                  for (k=0;k<p.numEqu;k++) {
113                                       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)]                                    for (m=0;m<p.numComp;m++) {
114                                      +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)]                                      rtmp=0;
115                                      +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)]                                      for (q=0;q<p.numQuadTotal;q++) {
116                                      +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)]);                                         rtmp+=Vol[q]* (
117                                               DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
118                                              +DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]
119                                              +DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,1,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
120                                              +DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,1,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]);
121                                        }
122                                        EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;
123                                      }
124                                    }
125                                }                                }
                               EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;  
126                              }                              }
127                            }                           } else {
128                          }                              for (s=0;s<p.row_numShapes;s++) {
129                        }                                for (r=0;r<p.col_numShapes;r++) {
130                     } else {                                   rtmp00=0;
131                        for (s=0;s<p.row_NS;s++) {                                   rtmp01=0;
132                          for (r=0;r<p.col_NS;r++) {                                   rtmp10=0;
133                             rtmp00=0;                                   rtmp11=0;
134                             rtmp01=0;                                   for (q=0;q<p.numQuadTotal;q++) {
135                             rtmp10=0;                                         rtmp0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
136                             rtmp11=0;                                         rtmp1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
137                             for (q=0;q<p.numQuad;q++) {                                         rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
138                                   rtmp0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                                         rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
139                                   rtmp1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                                         rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
140                                   rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                                         rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
141                                   rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                                   }
142                                   rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                                   for (k=0;k<p.numEqu;k++) {
143                                   rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                                       for (m=0;m<p.numComp;m++) {
144                             }                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=
145                             for (k=0;k<p.numEqu;k++) {                                                     rtmp00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]
146                                 for (m=0;m<p.numComp;m++) {                                                    +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]
147                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=                                                    +rtmp10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]
148                                               rtmp00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]                                                    +rtmp11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)];
149                                              +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]                                       }
150                                              +rtmp10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]                                   }
151                                              +rtmp11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)];                                }
                                }  
                            }  
                         }  
                       }  
                   }  
                 }  
                 /**************************************************************/  
                 /*   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)]);  
                                }  
                                EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;  
152                              }                              }
                           }  
153                          }                          }
154                        }                        }
155                     } else {                        /**************************************************************/
156                        for (s=0;s<p.row_NS;s++) {                        /*   process B: */
157                          for (r=0;r<p.col_NS;r++) {                        /**************************************************************/
158                              rtmp0=0;                        B_p=getSampleDataRO(B,e);
159                              rtmp1=0;                        if (NULL!=B_p) {
160                              for (q=0;q<p.numQuad;q++) {                           add_EM_S=TRUE;
161                                  rtmp=Vol[q]*S[INDEX2(r,q,p.row_NS)];                           if (extendedB) {
162                                  rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                      B_q=&(B_p[INDEX5(0,0,0,0,0 ,p.numEqu,DIM,p.numComp,p.numQuadTotal)]);
163                                  rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                              for (s=0;s<p.row_numShapes;s++) {
164                                  for (r=0;r<p.col_numShapes;r++) {
165                                    for (k=0;k<p.numEqu;k++) {
166                                      for (m=0;m<p.numComp;m++) {
167                                         rtmp=0;
168                                         for (q=0;q<p.numQuadTotal;q++) {
169                                            rtmp+=Vol[q]*S[INDEX2(r,q,p.row_numShapes)]*
170                                                 ( DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]
171                                                 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]);
172                                         }
173                                         EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;
174                                      }
175                                    }
176                                  }
177                              }                              }
178                              for (k=0;k<p.numEqu;k++) {                           } else {
179                                 for (m=0;m<p.numComp;m++) {                              for (s=0;s<p.row_numShapes;s++) {
180                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+= rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]                                for (r=0;r<p.col_numShapes;r++) {
181                                                                                      + rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)];                                    rtmp0=0;
182                                 }                                    rtmp1=0;
183                                      for (q=0;q<p.numQuadTotal;q++) {
184                                          rtmp=Vol[q]*S[INDEX2(r,q,p.row_numShapes)];
185                                          rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
186                                          rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
187                                      }
188                                      for (k=0;k<p.numEqu;k++) {
189                                         for (m=0;m<p.numComp;m++) {
190                                            EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+= rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]
191                                                                                              + rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)];
192                                         }
193                                      }
194                                  }
195                              }                              }
196                          }                           }
197                        }                        }
198                     }                        /**************************************************************/
199                  }                        /*   process C: */
200                  /**************************************************************/                        /**************************************************************/
201                  /*   process C: */                        C_p=getSampleDataRO(C,e);
202                  /**************************************************************/                        if (NULL!=C_p) {
203                  C_p=getSampleData(C,e);                          add_EM_S=TRUE;
204                  if (NULL!=C_p) {                          if (extendedC) {
205                    add_EM_S=TRUE;                  C_q=&(C_p[INDEX5(0,0,0,0,0, p.numEqu,p.numComp,DIM,p.numQuadTotal)]);
206                    if (extendedC) {                              for (s=0;s<p.row_numShapes;s++) {
207                        for (s=0;s<p.row_NS;s++) {                                for (r=0;r<p.col_numShapes;r++) {
208                          for (r=0;r<p.col_NS;r++) {                                  for (k=0;k<p.numEqu;k++) {
209                            for (k=0;k<p.numEqu;k++) {                                    for (m=0;m<p.numComp;m++) {
210                              for (m=0;m<p.numComp;m++) {                                       rtmp=0;
211                                 rtmp=0;                                       for (q=0;q<p.numQuadTotal;q++) {
212                                 for (q=0;q<p.numQuad;q++) {                                           rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*
213                                     rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*                                                ( C_q[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
214                                          ( C_p[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]                                                + C_q[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]);
215                                          + C_p[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]);                                       }
216                                 }                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;
217                                 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                                    }
218                                    }
219                                  }
220                              }                              }
221                            }                          } else {
222                          }                              for (s=0;s<p.row_numShapes;s++) {
223                        }                                for (r=0;r<p.col_numShapes;r++) {
224                    } else {                                   rtmp0=0;
225                        for (s=0;s<p.row_NS;s++) {                                   rtmp1=0;
226                          for (r=0;r<p.col_NS;r++) {                                   for (q=0;q<p.numQuadTotal;q++) {
227                             rtmp0=0;                                           rtmp=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
228                             rtmp1=0;                                           rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
229                             for (q=0;q<p.numQuad;q++) {                                           rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
230                                     rtmp=Vol[q]*S[INDEX2(s,q,p.row_NS)];                                   }
231                                     rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                                   for (k=0;k<p.numEqu;k++) {
232                                     rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                                      for (m=0;m<p.numComp;m++) {
233                             }                                             EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp0*C_p[INDEX3(k,m,0,p.numEqu,p.numComp)]
234                             for (k=0;k<p.numEqu;k++) {                                                                                               +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)];
235                                for (m=0;m<p.numComp;m++) {                                      }
236                                       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)]                                   }
                                                                                        +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)];  
237                                }                                }
238                             }                              }
239                          }                          }
240                        }                        }
241                    }                        /************************************************************* */
242                  }                        /* process D */
243                  /************************************************************* */                        /**************************************************************/
244                  /* process D */                        D_p=getSampleDataRO(D,e);
245                  /**************************************************************/                        if (NULL!=D_p) {
246                  D_p=getSampleData(D,e);                          add_EM_S=TRUE;
247                  if (NULL!=D_p) {                          if (extendedD) {
248                    add_EM_S=TRUE;                  D_q=&(D_p[INDEX4(0,0,0,0, p.numEqu,p.numComp,p.numQuadTotal)]);
249                    if (extendedD) {                              for (s=0;s<p.row_numShapes;s++) {
250                        for (s=0;s<p.row_NS;s++) {                                for (r=0;r<p.col_numShapes;r++) {
251                          for (r=0;r<p.col_NS;r++) {                                  for (k=0;k<p.numEqu;k++) {
252                            for (k=0;k<p.numEqu;k++) {                                    for (m=0;m<p.numComp;m++) {
253                              for (m=0;m<p.numComp;m++) {                                      rtmp=0;
254                                rtmp=0;                                      for (q=0;q<p.numQuadTotal;q++) {
255                                for (q=0;q<p.numQuad;q++) {                                          rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_q[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_numShapes)];
256                                    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)];                                      }
257                                        EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;
258                                      }
259                                    }
260                                }                                }
                               EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;  
261                              }                              }
262                            }                          } else {
263                          }                              for (s=0;s<p.row_numShapes;s++) {
264                        }                                for (r=0;r<p.col_numShapes;r++) {
265                    } else {                                    rtmp=0;
266                        for (s=0;s<p.row_NS;s++) {                                    for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(r,q,p.row_numShapes)];
267                          for (r=0;r<p.col_NS;r++) {                                    for (k=0;k<p.numEqu;k++) {
268                              rtmp=0;                                        for (m=0;m<p.numComp;m++) {
269                              for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];
270                              for (k=0;k<p.numEqu;k++) {                                       }
271                                  for (m=0;m<p.numComp;m++) {                                    }
272                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];                                }
                                }  
273                              }                              }
274                          }                          }
275                        }                        }
276                    }                        /**************************************************************/
277                  }                        /*   process X: */
278                  /**************************************************************/                        /**************************************************************/
279                  /*   process X: */                        X_p=getSampleDataRO(X,e);
280                  /**************************************************************/                        if (NULL!=X_p) {
281                  X_p=getSampleData(X,e);                          add_EM_F=TRUE;
282                  if (NULL!=X_p) {                          if (extendedX) {
283                    add_EM_F=TRUE;                 X_q=&(X_p[INDEX4(0,0,0,0, p.numEqu,DIM,p.numQuadTotal)]);
284                    if (extendedX) {                             for (s=0;s<p.row_numShapes;s++) {
285                       for (s=0;s<p.row_NS;s++) {                                for (k=0;k<p.numEqu;k++) {
286                          for (k=0;k<p.numEqu;k++) {                                  rtmp=0;
287                            rtmp=0;                                  for (q=0;q<p.numQuadTotal;q++) {
288                            for (q=0;q<p.numQuad;q++) {                                       rtmp+=Vol[q]*(DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,0,q,p.numEqu,DIM)]
289                                 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_numShapesTotal,DIM)]*X_q[INDEX3(k,1,q,p.numEqu,DIM)]);
290                                              +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*X_p[INDEX3(k,1,q,p.numEqu,DIM)]);                                  }
291                            }                                  EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
292                            EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;                                }
293                          }                             }
294                       }                          } else {
295                    } else {                             for (s=0;s<p.row_numShapes;s++) {
296                       for (s=0;s<p.row_NS;s++) {                                   rtmp0=0;
297                             rtmp0=0;                                   rtmp1=0;
298                             rtmp1=0;                                   for (q=0;q<p.numQuadTotal;q++) {
299                             for (q=0;q<p.numQuad;q++) {                                       rtmp0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
300                                 rtmp0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                                       rtmp1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
301                                 rtmp1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                                   }
302                                     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)];
303                             }                             }
                            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)];  
                      }  
                   }  
                }  
                /**************************************************************/  
                /*   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;  
304                          }                          }
305                       }                       }
306                     } else {                       /**************************************************************/
307                       for (s=0;s<p.row_NS;s++) {                       /*   process Y: */
308                           rtmp=0;                       /**************************************************************/
309                           for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];                        Y_p=getSampleDataRO(Y,e);
310                           for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];                        if (NULL!=Y_p) {
311                       }                          add_EM_F=TRUE;
312                     }                          if (extendedY) {
313                   }                 Y_q=&(Y_p[INDEX3(0,0,0, p.numEqu,p.numQuadTotal)]);
314                   /***********************************************************************************************/                             for (s=0;s<p.row_numShapes;s++) {
315                   /* add the element matrices onto the matrix and right hand side                                */                                for (k=0;k<p.numEqu;k++) {
316                   /***********************************************************************************************/                                   rtmp=0;
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)]];                                   for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*Y_q[INDEX2(k,q,p.numEqu)];
318                   if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);                                   EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
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                               }
321               } /* end color check */                           } else {
322            } /* end element loop */                             for (s=0;s<p.row_numShapes;s++) {
323        } /* end color loop */                                 rtmp=0;
324                                   for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
325                                   for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];
326                               }
327                             }
328                           }
329                           /***********************************************************************************************/
330                           /* add the element matrices onto the matrix and right hand side                                */
331                           /***********************************************************************************************/
332                           for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
333    
334                           if (add_EM_F) Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);
335                           if (add_EM_S) Dudley_Assemble_addToSystemMatrix(Mat,p.row_numShapesTotal,row_index,p.numEqu,p.col_numShapesTotal,row_index,p.numComp,EM_S);
336    
337                    } /* end color check */
338                 } /* end element loop */
339             } /* end color loop */
340              
341             THREAD_MEMFREE(EM_S);
342             THREAD_MEMFREE(EM_F);
343             THREAD_MEMFREE(row_index);
344    
345          } /* end of pointer check */
346     } /* end parallel region */     } /* end parallel region */
347  }  }
348  /*  /*

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

  ViewVC Help
Powered by ViewVC 1.1.26