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

Diff of /branches/trilinos_from_5897/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 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC branches/domexper/dudley/src/Assemble_PDE_System2_2D.c revision 3203 by jfenwick, Thu Sep 23 23:51:26 2010 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2009 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 11  Line 11 
11  *  *
12  *******************************************************/  *******************************************************/
13    
   
14  /**************************************************************/  /**************************************************************/
15    
16  /*    assembles the system of numEq PDEs into the stiffness matrix S right hand side F  */  /*    assembles the system of numEq PDEs into the stiffness matrix S right hand side F  */
17  /*    the shape functions for test and solution must be identical */  /*    the shape functions for test and solution must be identical */
18    
   
19  /*      -(A_{k,i,m,j} u_m,j)_i-(B_{k,i,m} u_m)_i+C_{k,m,j} u_m,j-D_{k,m} u_m  and -(X_{k,i})_i + Y_k */  /*      -(A_{k,i,m,j} u_m,j)_i-(B_{k,i,m} u_m)_i+C_{k,m,j} u_m,j-D_{k,m} u_m  and -(X_{k,i})_i + Y_k */
20    
21  /*    u has p.numComp components in a 2D domain. The shape functions for test and solution must be identical  */  /*    u has p.numComp components in a 2D domain. The shape functions for test and solution must be identical  */
# Line 32  Line 30 
30  /*      X = p.numEqu x 2  */  /*      X = p.numEqu x 2  */
31  /*      Y = p.numEqu   */  /*      Y = p.numEqu   */
32    
   
33  /**************************************************************/  /**************************************************************/
34    
   
35  #include "Assemble.h"  #include "Assemble.h"
36  #include "Util.h"  #include "Util.h"
37  #ifdef _OPENMP  #ifdef _OPENMP
38  #include <omp.h>  #include <omp.h>
39  #endif  #endif
40    
   
41  /**************************************************************/  /**************************************************************/
42    
43  void  Finley_Assemble_PDE_System2_2D(Assemble_Parameters p, Finley_ElementFile* elements,  void Dudley_Assemble_PDE_System2_2D(Assemble_Parameters p, Dudley_ElementFile * elements,
44                                       Paso_SystemMatrix* Mat, escriptDataC* F,                      Paso_SystemMatrix * Mat, escriptDataC * F,
45                                       escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y) {                      escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,
46                        escriptDataC * X, escriptDataC * Y)
47    {
48    
49      #define DIM 2  #define DIM 2
50      index_t color;      index_t color;
51      dim_t e;      dim_t e;
52      __const double *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p;      __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;
53      double *EM_S, *EM_F, *Vol, *DSDX;      double *EM_S, *EM_F, *DSDX;
54      index_t *row_index;      index_t *row_index;
55      register dim_t q, s,r,k,m;      register dim_t q, s, r, k, m;
56      register double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11;      register double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11;
57      bool_t add_EM_F, add_EM_S;      bool_t add_EM_F, add_EM_S;
58    
59      bool_t extendedA=isExpanded(A);      bool_t extendedA = isExpanded(A);
60      bool_t extendedB=isExpanded(B);      bool_t extendedB = isExpanded(B);
61      bool_t extendedC=isExpanded(C);      bool_t extendedC = isExpanded(C);
62      bool_t extendedD=isExpanded(D);      bool_t extendedD = isExpanded(D);
63      bool_t extendedX=isExpanded(X);      bool_t extendedX = isExpanded(X);
64      bool_t extendedY=isExpanded(Y);      bool_t extendedY = isExpanded(Y);
65      double *F_p=(requireWrite(F), getSampleDataRW(F,0));    /* use comma, to get around the mixed code and declarations thing */      double *F_p = (requireWrite(F), getSampleDataRW(F, 0)); /* use comma, to get around the mixed code and declarations thing */
66      double *S=p.row_jac->ReferenceElement->S;  //    double *S = p.row_jac->BasisFunctions->S;
67      dim_t len_EM_S=p.row_NN*p.col_NN*p.numEqu*p.numComp;      const double* S = p.shapeFns;
68      dim_t len_EM_F=p.row_NN*p.numEqu;      dim_t len_EM_S = p.row_numShapesTotal * p.row_numShapesTotal * p.numEqu * p.numComp;
69        dim_t len_EM_F = p.row_numShapesTotal * p.numEqu;
70      void* ABuff=allocSampleBuffer(A);  
71      void* BBuff=allocSampleBuffer(B);  #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)
     void* CBuff=allocSampleBuffer(C);  
     void* DBuff=allocSampleBuffer(D);  
     void* XBuff=allocSampleBuffer(X);  
     void* YBuff=allocSampleBuffer(Y);  
     #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)  
72      {      {
73    
74         EM_S=THREAD_MEMALLOC(len_EM_S,double);      EM_S = THREAD_MEMALLOC(len_EM_S, double);
75         EM_F=THREAD_MEMALLOC(len_EM_F,double);      EM_F = THREAD_MEMALLOC(len_EM_F, double);
76         row_index=THREAD_MEMALLOC(p.row_NN,index_t);      row_index = THREAD_MEMALLOC(p.row_numShapesTotal, index_t);
77                                                                                                                                                                                                        
78         if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) {      if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
79        {
80            for (color=elements->minColor;color<=elements->maxColor;color++) {  
81               /*  open loop over all elements: */          for (color = elements->minColor; color <= elements->maxColor; color++)
82               #pragma omp for private(e) schedule(static)          {
83               for(e=0;e<elements->numElements;e++){          /*  open loop over all elements: */
84                  if (elements->Color[e]==color) {  #pragma omp for private(e) schedule(static)
85                     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);          for (e = 0; e < elements->numElements; e++)
86                     DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);          {
87                     for (q=0;q<len_EM_S;++q) EM_S[q]=0;              if (elements->Color[e] == color)
88                     for (q=0;q<len_EM_F;++q) EM_F[q]=0;              {
89                     add_EM_F=FALSE;  
90                     add_EM_S=FALSE;              A_p = getSampleDataRO(A, e);
91                     /**************************************************************/              B_p = getSampleDataRO(B, e);
92                     /*   process A: */              C_p = getSampleDataRO(C, e);
93                     /**************************************************************/              D_p = getSampleDataRO(D, e);
94                     A_p=getSampleDataRO(A,e,ABuff);              X_p = getSampleDataRO(X, e);
95                     if (NULL!=A_p) {              Y_p = getSampleDataRO(Y, e);
96                        add_EM_S=TRUE;              double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
97                        if (extendedA) {              DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapesTotal, DIM, p.numQuadTotal, 1)]);
98                           for (s=0;s<p.row_NS;s++) {              for (q = 0; q < len_EM_S; ++q)
99                             for (r=0;r<p.col_NS;r++) {                  EM_S[q] = 0;
100                               for (k=0;k<p.numEqu;k++) {              for (q = 0; q < len_EM_F; ++q)
101                                 for (m=0;m<p.numComp;m++) {                  EM_F[q] = 0;
102                                   rtmp=0;              add_EM_F = FALSE;
103                                   for (q=0;q<p.numQuad;q++) {              add_EM_S = FALSE;
104                                      rtmp+=Vol[q]* (  
105                                          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)]                /**************************************************************/
106                                         +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)]              /*   process A: */
107                                         +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)]                /**************************************************************/
108                                         +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)]);              A_p = getSampleDataRO(A, e);
109                                   }              if (NULL != A_p)
110                                   EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;              {
111                                 }                  add_EM_S = TRUE;
112                               }                  if (extendedA)
113                             }                  {
114                           }                  A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuadTotal)]);
115                        } else {                  for (s = 0; s < p.row_numShapes; s++)
116                           for (s=0;s<p.row_NS;s++) {                  {
117                             for (r=0;r<p.col_NS;r++) {                      for (r = 0; r < p.row_numShapes; r++)
118                                rtmp00=0;                      {
119                                rtmp01=0;                      for (k = 0; k < p.numEqu; k++)
120                                rtmp10=0;                      {
121                                rtmp11=0;                          for (m = 0; m < p.numComp; m++)
122                                for (q=0;q<p.numQuad;q++) {                          {
123                                      rtmp0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                          rtmp = 0;
124                                      rtmp1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                          for (q = 0; q < p.numQuadTotal; q++)
125                                      rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                          {
126                                      rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                              rtmp +=
127                                      rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                              vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
128                                      rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                                     A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
129                                }                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
130                                for (k=0;k<p.numEqu;k++) {                                     DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
131                                    for (m=0;m<p.numComp;m++) {                                     A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
132                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
133                                                  rtmp00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
134                                                 +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]                                     A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
135                                                 +rtmp10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
136                                                 +rtmp11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)];                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
137                                    }                                     A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
138                                }                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);
139                             }                          }
140                           }                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
141                       }                              rtmp;
142                     }                          }
143                     /**************************************************************/                      }
144                     /*   process B: */                      }
145                     /**************************************************************/                  }
146                     B_p=getSampleDataRO(B,e,BBuff);                  } else
147                     if (NULL!=B_p) {                  {
148                        add_EM_S=TRUE;                  for (s = 0; s < p.row_numShapes; s++)
149                        if (extendedB) {                  {
150                           for (s=0;s<p.row_NS;s++) {                      for (r = 0; r < p.row_numShapes; r++)
151                             for (r=0;r<p.col_NS;r++) {                      {
152                               for (k=0;k<p.numEqu;k++) {                      rtmp00 = 0;
153                                 for (m=0;m<p.numComp;m++) {                      rtmp01 = 0;
154                                    rtmp=0;                      rtmp10 = 0;
155                                    for (q=0;q<p.numQuad;q++) {                      rtmp11 = 0;
156                                       rtmp+=Vol[q]*S[INDEX2(r,q,p.row_NS)]*                      for (q = 0; q < p.numQuadTotal; q++)
157                                            ( DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*B_p[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]                      {
158                                            + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*B_p[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]);                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
159                                    }                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
160                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
161                                 }                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
162                               }                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
163                             }                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
164                           }                      }
165                        } else {                      for (k = 0; k < p.numEqu; k++)
166                           for (s=0;s<p.row_NS;s++) {                      {
167                             for (r=0;r<p.col_NS;r++) {                          for (m = 0; m < p.numComp; m++)
168                                 rtmp0=0;                          {
169                                 rtmp1=0;                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
170                                 for (q=0;q<p.numQuad;q++) {                              rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)]
171                                     rtmp=Vol[q]*S[INDEX2(r,q,p.row_NS)];                              + rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)]
172                                     rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                              + rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)]
173                                     rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                              + rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)];
174                                 }                          }
175                                 for (k=0;k<p.numEqu;k++) {                      }
176                                    for (m=0;m<p.numComp;m++) {                      }
177                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+= rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]                  }
178                                                                                         + rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)];                  }
179                                    }              }
180                                 }                /**************************************************************/
181                             }              /*   process B: */
182                           }                /**************************************************************/
183                        }              B_p = getSampleDataRO(B, e);
184                     }              if (NULL != B_p)
185                     /**************************************************************/              {
186                     /*   process C: */                  add_EM_S = TRUE;
187                     /**************************************************************/                  if (extendedB)
188                     C_p=getSampleDataRO(C,e,CBuff);                  {
189                     if (NULL!=C_p) {                  B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);
190                       add_EM_S=TRUE;                  for (s = 0; s < p.row_numShapes; s++)
191                       if (extendedC) {                  {
192                           for (s=0;s<p.row_NS;s++) {                      for (r = 0; r < p.row_numShapes; r++)
193                             for (r=0;r<p.col_NS;r++) {                      {
194                               for (k=0;k<p.numEqu;k++) {                      for (k = 0; k < p.numEqu; k++)
195                                 for (m=0;m<p.numComp;m++) {                      {
196                                    rtmp=0;                          for (m = 0; m < p.numComp; m++)
197                                    for (q=0;q<p.numQuad;q++) {                          {
198                                        rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*                          rtmp = 0;
199                                             ( C_p[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]                          for (q = 0; q < p.numQuadTotal; q++)
200                                             + C_p[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]);                          {
201                                    }                              rtmp += vol * S[INDEX2(r, q, p.row_numShapes)] *
202                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                              (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
203                                 }                               B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
204                               }                               DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
205                             }                               B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);
206                           }                          }
207                       } else {                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
208                           for (s=0;s<p.row_NS;s++) {                              rtmp;
209                             for (r=0;r<p.col_NS;r++) {                          }
210                                rtmp0=0;                      }
211                                rtmp1=0;                      }
212                                for (q=0;q<p.numQuad;q++) {                  }
213                                        rtmp=Vol[q]*S[INDEX2(s,q,p.row_NS)];                  } else
214                                        rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                  {
215                                        rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                  for (s = 0; s < p.row_numShapes; s++)
216                                }                  {
217                                for (k=0;k<p.numEqu;k++) {                      for (r = 0; r < p.row_numShapes; r++)
218                                   for (m=0;m<p.numComp;m++) {                      {
219                                          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)]                      rtmp0 = 0;
220                                                                                            +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)];                      rtmp1 = 0;
221                                   }                      for (q = 0; q < p.numQuadTotal; q++)
222                                }                      {
223                             }                          rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];
224                           }                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
225                       }                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
226                     }                      }
227                     /************************************************************* */                      for (k = 0; k < p.numEqu; k++)
228                     /* process D */                      {
229                     /**************************************************************/                          for (m = 0; m < p.numComp; m++)
230                     D_p=getSampleDataRO(D,e,DBuff);                          {
231                     if (NULL!=D_p) {                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
232                       add_EM_S=TRUE;                              rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
233                       if (extendedD) {                              rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];
234                           for (s=0;s<p.row_NS;s++) {                          }
235                             for (r=0;r<p.col_NS;r++) {                      }
236                               for (k=0;k<p.numEqu;k++) {                      }
237                                 for (m=0;m<p.numComp;m++) {                  }
238                                   rtmp=0;                  }
239                                   for (q=0;q<p.numQuad;q++) {              }
240                                       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)];                /**************************************************************/
241                                   }              /*   process C: */
242                                   EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                /**************************************************************/
243                                 }              C_p = getSampleDataRO(C, e);
244                               }              if (NULL != C_p)
245                             }              {
246                           }                  add_EM_S = TRUE;
247                       } else {                  if (extendedC)
248                           for (s=0;s<p.row_NS;s++) {                  {
249                             for (r=0;r<p.col_NS;r++) {                  C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);
250                                 rtmp=0;                  for (s = 0; s < p.row_numShapes; s++)
251                                 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];                  {
252                                 for (k=0;k<p.numEqu;k++) {                      for (r = 0; r < p.row_numShapes; r++)
253                                     for (m=0;m<p.numComp;m++) {                      {
254                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];                      for (k = 0; k < p.numEqu; k++)
255                                    }                      {
256                                 }                          for (m = 0; m < p.numComp; m++)
257                             }                          {
258                           }                          rtmp = 0;
259                       }                          for (q = 0; q < p.numQuadTotal; q++)
260                     }                          {
261                     /**************************************************************/                              rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] *
262                     /*   process X: */                              (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
263                     /**************************************************************/                               DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
264                     X_p=getSampleDataRO(X,e,XBuff);                               C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
265                     if (NULL!=X_p) {                               DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);
266                       add_EM_F=TRUE;                          }
267                       if (extendedX) {                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
268                          for (s=0;s<p.row_NS;s++) {                              rtmp;
269                             for (k=0;k<p.numEqu;k++) {                          }
270                               rtmp=0;                      }
271                               for (q=0;q<p.numQuad;q++) {                      }
272                                    rtmp+=Vol[q]*(DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX3(k,0,q,p.numEqu,DIM)]                  }
273                                                 +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*X_p[INDEX3(k,1,q,p.numEqu,DIM)]);                  } else
274                               }                  {
275                               EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;                  for (s = 0; s < p.row_numShapes; s++)
276                             }                  {
277                          }                      for (r = 0; r < p.row_numShapes; r++)
278                       } else {                      {
279                          for (s=0;s<p.row_NS;s++) {                      rtmp0 = 0;
280                                rtmp0=0;                      rtmp1 = 0;
281                                rtmp1=0;                      for (q = 0; q < p.numQuadTotal; q++)
282                                for (q=0;q<p.numQuad;q++) {                      {
283                                    rtmp0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                          rtmp = vol * S[INDEX2(s, q, p.row_numShapes)];
284                                    rtmp1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
285                                }                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
286                                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)];                      }
287                          }                      for (k = 0; k < p.numEqu; k++)
288                       }                      {
289                    }                          for (m = 0; m < p.numComp; m++)
290                    /**************************************************************/                          {
291                    /*   process Y: */                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
292                    /**************************************************************/                              rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
293                     Y_p=getSampleDataRO(Y,e,YBuff);                              rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];
294                     if (NULL!=Y_p) {                          }
295                       add_EM_F=TRUE;                      }
296                       if (extendedY) {                      }
297                          for (s=0;s<p.row_NS;s++) {                  }
298                             for (k=0;k<p.numEqu;k++) {                  }
299                                rtmp=0;              }
300                                for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[INDEX2(k,q,p.numEqu)];                /************************************************************* */
301                                EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;              /* process D */
302                             }                /**************************************************************/
303                          }              D_p = getSampleDataRO(D, e);
304                        } else {              if (NULL != D_p)
305                          for (s=0;s<p.row_NS;s++) {              {
306                              rtmp=0;                  add_EM_S = TRUE;
307                              for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];                  if (extendedD)
308                              for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];                  {
309                          }                  D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);
310                        }                  for (s = 0; s < p.row_numShapes; s++)
311                      }                  {
312                      /***********************************************************************************************/                      for (r = 0; r < p.row_numShapes; r++)
313                      /* add the element matrices onto the matrix and right hand side                                */                      {
314                      /***********************************************************************************************/                      for (k = 0; k < p.numEqu; k++)
315                      for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];                      {
316                      if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);                          for (m = 0; m < p.numComp; m++)
317                      if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);                          {
318                              rtmp = 0;
319                  } /* end color check */                          for (q = 0; q < p.numQuadTotal; q++)
320               } /* end element loop */                          {
321           } /* end color loop */                              rtmp +=
322                                          vol * S[INDEX2(s, q, p.row_numShapes)] *
323           THREAD_MEMFREE(EM_S);                              D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
324           THREAD_MEMFREE(EM_F);                              S[INDEX2(r, q, p.row_numShapes)];
325           THREAD_MEMFREE(row_index);                          }
326                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
327        } /* end of pointer check */                              rtmp;
328     } /* end parallel region */                          }
329     freeSampleBuffer(ABuff);                      }
330     freeSampleBuffer(BBuff);                      }
331     freeSampleBuffer(CBuff);                  }
332     freeSampleBuffer(DBuff);                  } else
333     freeSampleBuffer(XBuff);                  {
334     freeSampleBuffer(YBuff);                  for (s = 0; s < p.row_numShapes; s++)
335                    {
336                        for (r = 0; r < p.row_numShapes; r++)
337                        {
338                        rtmp = 0;
339                        for (q = 0; q < p.numQuadTotal; q++)
340                            rtmp +=
341                            vol * S[INDEX2(s, q, p.row_numShapes)] *
342                            S[INDEX2(r, q, p.row_numShapes)];
343                        for (k = 0; k < p.numEqu; k++)
344                        {
345                            for (m = 0; m < p.numComp; m++)
346                            {
347                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
348                                rtmp * D_p[INDEX2(k, m, p.numEqu)];
349                            }
350                        }
351                        }
352                    }
353                    }
354                }
355                  /**************************************************************/
356                /*   process X: */
357                  /**************************************************************/
358                X_p = getSampleDataRO(X, e);
359                if (NULL != X_p)
360                {
361                    add_EM_F = TRUE;
362                    if (extendedX)
363                    {
364                    X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuadTotal)]);
365                    for (s = 0; s < p.row_numShapes; s++)
366                    {
367                        for (k = 0; k < p.numEqu; k++)
368                        {
369                        rtmp = 0;
370                        for (q = 0; q < p.numQuadTotal; q++)
371                        {
372                            rtmp +=
373                            vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
374                                   X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
375                                   DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
376                                   X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);
377                        }
378                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
379                        }
380                    }
381                    } else
382                    {
383                    for (s = 0; s < p.row_numShapes; s++)
384                    {
385                        rtmp0 = 0;
386                        rtmp1 = 0;
387                        for (q = 0; q < p.numQuadTotal; q++)
388                        {
389                        rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
390                        rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
391                        }
392                        for (k = 0; k < p.numEqu; k++)
393                        EM_F[INDEX2(k, s, p.numEqu)] +=
394                            rtmp0 * X_p[INDEX2(k, 0, p.numEqu)] + rtmp1 * X_p[INDEX2(k, 1, p.numEqu)];
395                    }
396                    }
397                }
398                 /**************************************************************/
399                /*   process Y: */
400                 /**************************************************************/
401                Y_p = getSampleDataRO(Y, e);
402                if (NULL != Y_p)
403                {
404                    add_EM_F = TRUE;
405                    if (extendedY)
406                    {
407                    Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);
408                    for (s = 0; s < p.row_numShapes; s++)
409                    {
410                        for (k = 0; k < p.numEqu; k++)
411                        {
412                        rtmp = 0;
413                        for (q = 0; q < p.numQuadTotal; q++)
414                            rtmp +=
415                            vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
416                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
417                        }
418                    }
419                    } else
420                    {
421                    for (s = 0; s < p.row_numShapes; s++)
422                    {
423                        rtmp = 0;
424                        for (q = 0; q < p.numQuadTotal; q++)
425                        rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
426                        for (k = 0; k < p.numEqu; k++)
427                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
428                    }
429                    }
430                }
431                   /***********************************************************************************************/
432                /* add the element matrices onto the matrix and right hand side                                */
433                   /***********************************************************************************************/
434                for (q = 0; q < p.row_numShapesTotal; q++)
435                    row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
436    
437                if (add_EM_F)
438                    Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
439                               p.row_DOF_UpperBound);
440                if (add_EM_S)
441                    Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
442                                      p.row_numShapesTotal, row_index, p.numComp, EM_S);
443    
444                }       /* end color check */
445            }       /* end element loop */
446            }           /* end color loop */
447    
448            THREAD_MEMFREE(EM_S);
449            THREAD_MEMFREE(EM_F);
450            THREAD_MEMFREE(row_index);
451    
452        }           /* end of pointer check */
453        }               /* end parallel region */
454  }  }
455    
456  /*  /*
457   * $Log$   * $Log$
458   */   */

Legend:
Removed from v.2548  
changed lines
  Added in v.3203

  ViewVC Help
Powered by ViewVC 1.1.26