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

revision 3184 by jfenwick, Wed Sep 15 00:23:42 2010 UTC revision 3187 by jfenwick, Thu Sep 16 02:57:17 2010 UTC
# 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  Dudley_Assemble_PDE_System2_2D(Assemble_Parameters p, Dudley_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, *A_q, *B_q, *C_q, *D_q, *X_q, *Y_q;      __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, *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;
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->BasisFunctions->S;      double *S = p.row_jac->BasisFunctions->S;
67      dim_t len_EM_S=p.row_numShapesTotal*p.col_numShapesTotal*p.numEqu*p.numComp;      dim_t len_EM_S = p.row_numShapesTotal * p.col_numShapesTotal * p.numEqu * p.numComp;
68      dim_t len_EM_F=p.row_numShapesTotal*p.numEqu;      dim_t len_EM_F = p.row_numShapesTotal * p.numEqu;
69
70      #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)  #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)
71      {      {
72
76
77         if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index) ) {      if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
78        {
79            for (color=elements->minColor;color<=elements->maxColor;color++) {
80               /*  open loop over all elements: */          for (color = elements->minColor; color <= elements->maxColor; color++)
81               #pragma omp for private(e) schedule(static)          {
82               for(e=0;e<elements->numElements;e++){          /*  open loop over all elements: */
83                  if (elements->Color[e]==color) {  #pragma omp for private(e) schedule(static)
84                                  for (e = 0; e < elements->numElements; e++)
85             A_p=getSampleDataRO(A,e);          {
86             B_p=getSampleDataRO(B,e);              if (elements->Color[e] == color)
87                     C_p=getSampleDataRO(C,e);                  {
88                     D_p=getSampleDataRO(D,e);
89                     X_p=getSampleDataRO(X,e);              A_p = getSampleDataRO(A, e);
90                     Y_p=getSampleDataRO(Y,e);              B_p = getSampleDataRO(B, e);
91            double vol=p.row_jac->absD[e]*p.row_jac->quadweight;                      C_p = getSampleDataRO(C, e);
92                        DSDX=&(p.row_jac->DSDX[INDEX5(0,0,0,0,e,p.row_numShapesTotal,DIM,p.numQuadTotal,1)]);              D_p = getSampleDataRO(D, e);
93                        for (q=0;q<len_EM_S;++q) EM_S[q]=0;              X_p = getSampleDataRO(X, e);
94                        for (q=0;q<len_EM_F;++q) EM_F[q]=0;              Y_p = getSampleDataRO(Y, e);
96                        add_EM_S=FALSE;              DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapesTotal, DIM, p.numQuadTotal, 1)]);
97                for (q = 0; q < len_EM_S; ++q)
98                        /**************************************************************/                  EM_S[q] = 0;
99                        /*   process A: */              for (q = 0; q < len_EM_F; ++q)
100                        /**************************************************************/                  EM_F[q] = 0;
102                        if (NULL!=A_p) {              add_EM_S = FALSE;
104                           if (extendedA) {                /**************************************************************/
105                  A_q=&(A_p[INDEX6(0,0,0,0,0,0, p.numEqu,DIM,p.numComp,DIM,p.numQuadTotal)]);              /*   process A: */
106                              for (s=0;s<p.row_numShapes;s++) {                /**************************************************************/
107                                for (r=0;r<p.col_numShapes;r++) {              A_p = getSampleDataRO(A, e);
108                                  for (k=0;k<p.numEqu;k++) {              if (NULL != A_p)
109                                    for (m=0;m<p.numComp;m++) {              {
111                                      for (q=0;q<p.numQuadTotal;q++) {                  if (extendedA)
112                                         rtmp+=vol* (                  {
113                                             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)]                  A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuadTotal)]);
114                                            +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)]                  for (s = 0; s < p.row_numShapes; s++)
115                                            +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)]                  {
116                                            +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)]);                      for (r = 0; r < p.col_numShapes; r++)
117                                      }                      {
118                                      EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                      for (k = 0; k < p.numEqu; k++)
119                                    }                      {
120                                  }                          for (m = 0; m < p.numComp; m++)
121                                }                          {
122                              }                          rtmp = 0;
123                           } else {                          for (q = 0; q < p.numQuadTotal; q++)
124                              for (s=0;s<p.row_numShapes;s++) {                          {
125                                for (r=0;r<p.col_numShapes;r++) {                              rtmp +=
126                                   rtmp00=0;                              vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
127                                   rtmp01=0;                                     A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
128                                   rtmp10=0;                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
129                                   rtmp11=0;                                     DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
130                                   for (q=0;q<p.numQuadTotal;q++) {                                     A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
131                                         rtmp0=vol*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
132                                         rtmp1=vol*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
133                                         rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                                     A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
134                                         rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
135                                         rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
136                                         rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                                     A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
137                                   }                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);
138                                   for (k=0;k<p.numEqu;k++) {                          }
139                                       for (m=0;m<p.numComp;m++) {                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
140                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=                              rtmp;
141                                                     rtmp00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]                          }
142                                                    +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]                      }
143                                                    +rtmp10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]                      }
144                                                    +rtmp11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)];                  }
145                                       }                  } else
146                                   }                  {
147                                }                  for (s = 0; s < p.row_numShapes; s++)
148                              }                  {
149                          }                      for (r = 0; r < p.col_numShapes; r++)
150                        }                      {
151                        /**************************************************************/                      rtmp00 = 0;
152                        /*   process B: */                      rtmp01 = 0;
153                        /**************************************************************/                      rtmp10 = 0;
154                        B_p=getSampleDataRO(B,e);                      rtmp11 = 0;
155                        if (NULL!=B_p) {                      for (q = 0; q < p.numQuadTotal; q++)
157                           if (extendedB) {                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
158                      B_q=&(B_p[INDEX5(0,0,0,0,0 ,p.numEqu,DIM,p.numComp,p.numQuadTotal)]);                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
159                              for (s=0;s<p.row_numShapes;s++) {                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
160                                for (r=0;r<p.col_numShapes;r++) {                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
161                                  for (k=0;k<p.numEqu;k++) {                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
162                                    for (m=0;m<p.numComp;m++) {                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
163                                       rtmp=0;                      }
164                                       for (q=0;q<p.numQuadTotal;q++) {                      for (k = 0; k < p.numEqu; k++)
165                                          rtmp+=vol*S[INDEX2(r,q,p.row_numShapes)]*                      {
166                                               ( DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]                          for (m = 0; m < p.numComp; m++)
167                                               + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]);                          {
168                                       }                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
169                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                              rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)]
170                                    }                              + rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)]
171                                  }                              + rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)]
172                                }                              + rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)];
173                              }                          }
174                           } else {                      }
175                              for (s=0;s<p.row_numShapes;s++) {                      }
176                                for (r=0;r<p.col_numShapes;r++) {                  }
177                                    rtmp0=0;                  }
178                                    rtmp1=0;              }
180                                        rtmp=vol*S[INDEX2(r,q,p.row_numShapes)];              /*   process B: */
181                                        rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                /**************************************************************/
182                                        rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];              B_p = getSampleDataRO(B, e);
183                                    }              if (NULL != B_p)
184                                    for (k=0;k<p.numEqu;k++) {              {
185                                       for (m=0;m<p.numComp;m++) {                  add_EM_S = TRUE;
186                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+= rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]                  if (extendedB)
187                                                                                            + rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)];                  {
188                                       }                  B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);
189                                    }                  for (s = 0; s < p.row_numShapes; s++)
190                                }                  {
191                              }                      for (r = 0; r < p.col_numShapes; r++)
192                           }                      {
193                        }                      for (k = 0; k < p.numEqu; k++)
194                        /**************************************************************/                      {
195                        /*   process C: */                          for (m = 0; m < p.numComp; m++)
196                        /**************************************************************/                          {
197                        C_p=getSampleDataRO(C,e);                          rtmp = 0;
198                        if (NULL!=C_p) {                          for (q = 0; q < p.numQuadTotal; q++)
200                          if (extendedC) {                              rtmp += vol * S[INDEX2(r, q, p.row_numShapes)] *
201                  C_q=&(C_p[INDEX5(0,0,0,0,0, p.numEqu,p.numComp,DIM,p.numQuadTotal)]);                              (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
202                              for (s=0;s<p.row_numShapes;s++) {                               B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
203                                for (r=0;r<p.col_numShapes;r++) {                               DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
204                                  for (k=0;k<p.numEqu;k++) {                               B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);
205                                    for (m=0;m<p.numComp;m++) {                          }
206                                       rtmp=0;                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
208                                           rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*                          }
209                                                ( C_q[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]                      }
210                                                + C_q[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]);                      }
211                                       }                  }
212                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                  } else
213                                    }                  {
214                                  }                  for (s = 0; s < p.row_numShapes; s++)
215                                }                  {
216                              }                      for (r = 0; r < p.col_numShapes; r++)
217                          } else {                      {
218                              for (s=0;s<p.row_numShapes;s++) {                      rtmp0 = 0;
219                                for (r=0;r<p.col_numShapes;r++) {                      rtmp1 = 0;
220                                   rtmp0=0;                      for (q = 0; q < p.numQuadTotal; q++)
221                                   rtmp1=0;                      {
222                                   for (q=0;q<p.numQuadTotal;q++) {                          rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];
223                                           rtmp=vol*S[INDEX2(s,q,p.row_numShapes)];                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
224                                           rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
225                                           rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                      }
226                                   }                      for (k = 0; k < p.numEqu; k++)
227                                   for (k=0;k<p.numEqu;k++) {                      {
228                                      for (m=0;m<p.numComp;m++) {                          for (m = 0; m < p.numComp; m++)
229                                             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)]                          {
230                                                                                               +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)];                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
231                                      }                              rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
232                                   }                              rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];
233                                }                          }
234                              }                      }
235                          }                      }
236                        }                  }
237                        /************************************************************* */                  }
238                        /* process D */              }
239                        /**************************************************************/                /**************************************************************/
240                        D_p=getSampleDataRO(D,e);              /*   process C: */
241                        if (NULL!=D_p) {                /**************************************************************/
242                          add_EM_S=TRUE;              C_p = getSampleDataRO(C, e);
243                          if (extendedD) {              if (NULL != C_p)
245                              for (s=0;s<p.row_numShapes;s++) {                  add_EM_S = TRUE;
246                                for (r=0;r<p.col_numShapes;r++) {                  if (extendedC)
247                                  for (k=0;k<p.numEqu;k++) {                  {
248                                    for (m=0;m<p.numComp;m++) {                  C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);
249                                      rtmp=0;                  for (s = 0; s < p.row_numShapes; s++)
251                                          rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*D_q[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_numShapes)];                      for (r = 0; r < p.col_numShapes; r++)
252                                      }                      {
253                                      EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                      for (k = 0; k < p.numEqu; k++)
254                                    }                      {
255                                  }                          for (m = 0; m < p.numComp; m++)
256                                }                          {
257                              }                          rtmp = 0;
258                          } else {                          for (q = 0; q < p.numQuadTotal; q++)
259                              for (s=0;s<p.row_numShapes;s++) {                          {
260                                for (r=0;r<p.col_numShapes;r++) {                              rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] *
261                                    rtmp=0;                              (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
262                                    for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(r,q,p.row_numShapes)];                               DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
263                                    for (k=0;k<p.numEqu;k++) {                               C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
264                                        for (m=0;m<p.numComp;m++) {                               DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);
265                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];                          }
266                                       }                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
267                                    }                              rtmp;
268                                }                          }
269                              }                      }
270                          }                      }
271                        }                  }
272                        /**************************************************************/                  } else
273                        /*   process X: */                  {
274                        /**************************************************************/                  for (s = 0; s < p.row_numShapes; s++)
275                        X_p=getSampleDataRO(X,e);                  {
276                        if (NULL!=X_p) {                      for (r = 0; r < p.col_numShapes; r++)
278                          if (extendedX) {                      rtmp0 = 0;
279                 X_q=&(X_p[INDEX4(0,0,0,0, p.numEqu,DIM,p.numQuadTotal)]);                      rtmp1 = 0;
280                             for (s=0;s<p.row_numShapes;s++) {                      for (q = 0; q < p.numQuadTotal; q++)
281                                for (k=0;k<p.numEqu;k++) {                      {
282                                  rtmp=0;                          rtmp = vol * S[INDEX2(s, q, p.row_numShapes)];
283                                  for (q=0;q<p.numQuadTotal;q++) {                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
284                                       rtmp+=vol*(DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,0,q,p.numEqu,DIM)]                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
285                                                    +DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,1,q,p.numEqu,DIM)]);                      }
286                                  }                      for (k = 0; k < p.numEqu; k++)
287                                  EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;                      {
288                                }                          for (m = 0; m < p.numComp; m++)
289                             }                          {
290                          } else {                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
291                             for (s=0;s<p.row_numShapes;s++) {                              rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
292                                   rtmp0=0;                              rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];
293                                   rtmp1=0;                          }
295                                       rtmp0+=vol*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                      }
296                                       rtmp1+=vol*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];                  }
297                                   }                  }
298                                   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)];              }
299                             }                /************************************************************* */
300                          }              /* process D */
301                       }                /**************************************************************/
302                       /**************************************************************/              D_p = getSampleDataRO(D, e);
303                       /*   process Y: */              if (NULL != D_p)
304                       /**************************************************************/              {
306                        if (NULL!=Y_p) {                  if (extendedD)
308                          if (extendedY) {                  D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);
309                 Y_q=&(Y_p[INDEX3(0,0,0, p.numEqu,p.numQuadTotal)]);                  for (s = 0; s < p.row_numShapes; s++)
310                             for (s=0;s<p.row_numShapes;s++) {                  {
311                                for (k=0;k<p.numEqu;k++) {                      for (r = 0; r < p.col_numShapes; r++)
312                                   rtmp=0;                      {
313                                   for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*Y_q[INDEX2(k,q,p.numEqu)];                      for (k = 0; k < p.numEqu; k++)
314                                   EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;                      {
315                                }                          for (m = 0; m < p.numComp; m++)
316                             }                          {
317                           } else {                          rtmp = 0;
318                             for (s=0;s<p.row_numShapes;s++) {                          for (q = 0; q < p.numQuadTotal; q++)
319                                 rtmp=0;                          {
320                                 for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)];                              rtmp +=
321                                 for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];                              vol * S[INDEX2(s, q, p.row_numShapes)] *
322                             }                              D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
323                           }                              S[INDEX2(r, q, p.row_numShapes)];
324                         }                          }
325                         /***********************************************************************************************/                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
326                         /* add the element matrices onto the matrix and right hand side                                */                              rtmp;
327                         /***********************************************************************************************/                          }
328                         for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];                      }
329                        }
332                    {
333                  } /* end color check */                  for (s = 0; s < p.row_numShapes; s++)
334               } /* end element loop */                  {
335           } /* end color loop */                      for (r = 0; r < p.col_numShapes; r++)
336                                  {
340                            vol * S[INDEX2(s, q, p.row_numShapes)] *
341                            S[INDEX2(r, q, p.row_numShapes)];
342                        for (k = 0; k < p.numEqu; k++)
343                        {
344                            for (m = 0; m < p.numComp; m++)
345                            {
346                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
347                                rtmp * D_p[INDEX2(k, m, p.numEqu)];
348                            }
349                        }
350                        }
351                    }
352                    }
353                }
354                  /**************************************************************/
355                /*   process X: */
356                  /**************************************************************/
357                X_p = getSampleDataRO(X, e);
358                if (NULL != X_p)
359                {
361                    if (extendedX)
362                    {
363                    X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuadTotal)]);
364                    for (s = 0; s < p.row_numShapes; s++)
365                    {
366                        for (k = 0; k < p.numEqu; k++)
367                        {
368                        rtmp = 0;
369                        for (q = 0; q < p.numQuadTotal; q++)
370                        {
371                            rtmp +=
372                            vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
373                                   X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
374                                   DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
375                                   X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);
376                        }
377                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
378                        }
379                    }
380                    } else
381                    {
382                    for (s = 0; s < p.row_numShapes; s++)
383                    {
384                        rtmp0 = 0;
385                        rtmp1 = 0;
386                        for (q = 0; q < p.numQuadTotal; q++)
387                        {
388                        rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
389                        rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
390                        }
391                        for (k = 0; k < p.numEqu; k++)
392                        EM_F[INDEX2(k, s, p.numEqu)] +=
393                            rtmp0 * X_p[INDEX2(k, 0, p.numEqu)] + rtmp1 * X_p[INDEX2(k, 1, p.numEqu)];
394                    }
395                    }
396                }
397                 /**************************************************************/
398                /*   process Y: */
399                 /**************************************************************/
400                Y_p = getSampleDataRO(Y, e);
401                if (NULL != Y_p)
402                {
404                    if (extendedY)
405                    {
406                    Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);
407                    for (s = 0; s < p.row_numShapes; s++)
408                    {
409                        for (k = 0; k < p.numEqu; k++)
410                        {
411                        rtmp = 0;
412                        for (q = 0; q < p.numQuadTotal; q++)
413                            rtmp +=
414                            vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
415                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
416                        }
417                    }
418                    } else
419                    {
420                    for (s = 0; s < p.row_numShapes; s++)
421                    {
422                        rtmp = 0;
423                        for (q = 0; q < p.numQuadTotal; q++)
424                        rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
425                        for (k = 0; k < p.numEqu; k++)
426                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
427                    }
428                    }
429                }
430                   /***********************************************************************************************/
431                /* add the element matrices onto the matrix and right hand side                                */
432                   /***********************************************************************************************/
433                for (q = 0; q < p.row_numShapesTotal; q++)
434                    row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
435
437                    Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
438                               p.row_DOF_UpperBound);
441                                      p.col_numShapesTotal, row_index, p.numComp, EM_S);
442
443                }       /* end color check */
444            }       /* end element loop */
445            }           /* end color loop */
446