# Diff of /trunk/dudley/src/Assemble_PDE_Single2_3D.c

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_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m  and -(X_,i)_i + Y */  /*      -(A_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m  and -(X_,i)_i + Y */
20
21  /*    in a 3D domain. The shape functions for test and solution must be identical  */  /*    in a 3D domain. The shape functions for test and solution must be identical  */
# Line 32  Line 30
30  /*      X = 3  */  /*      X = 3  */
31  /*      Y = scalar   */  /*      Y = scalar   */
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_Single2_3D(Assemble_Parameters p, Dudley_ElementFile* elements,  void Dudley_Assemble_PDE_Single2_3D(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 3  #define DIM 3
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;      register dim_t q, s, r;
56      register double rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2;      register double rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2;
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;      dim_t len_EM_S = p.row_numShapesTotal * p.col_numShapesTotal;
68      dim_t len_EM_F=p.row_numShapesTotal;      dim_t len_EM_F = p.row_numShapesTotal;
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,rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2,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,rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2,add_EM_F, add_EM_S)
71      {      {
75
76         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))
77        {
78            for (color=elements->minColor;color<=elements->maxColor;color++) {
79               /*  open loop over all elements: */          for (color = elements->minColor; color <= elements->maxColor; color++)
80               #pragma omp for private(e) schedule(static)          {
81               for(e=0;e<elements->numElements;e++){          /*  open loop over all elements: */
82                  if (elements->Color[e]==color) {  #pragma omp for private(e) schedule(static)
83            for (e = 0; e < elements->numElements; e++)
84                     A_p=getSampleDataRO(A,e);          {
85                     B_p=getSampleDataRO(B,e);              if (elements->Color[e] == color)
86                     C_p=getSampleDataRO(C,e);              {
87                     D_p=getSampleDataRO(D,e);
88                     X_p=getSampleDataRO(X,e);              A_p = getSampleDataRO(A, e);
89                     Y_p=getSampleDataRO(Y,e);              B_p = getSampleDataRO(B, e);
90                C_p = getSampleDataRO(C, e);
91                D_p = getSampleDataRO(D, e);
92                X_p = getSampleDataRO(X, e);
93                Y_p = getSampleDataRO(Y, e);
94
97                        DSDX=&(p.row_jac->DSDX[INDEX5(0,0,0,0,e,p.row_numShapesTotal,DIM,p.numQuadTotal,1)]);              DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapesTotal, DIM, p.numQuadTotal, 1)]);
98                        for (q=0;q<len_EM_S;++q) EM_S[q]=0;              for (q = 0; q < len_EM_S; ++q)
99                        for (q=0;q<len_EM_F;++q) EM_F[q]=0;                  EM_S[q] = 0;
100                        add_EM_F=FALSE;              for (q = 0; q < len_EM_F; ++q)
104                        /*   process A: */
105                        /**************************************************************/                /**************************************************************/
106                        if (NULL!=A_p) {              /*   process A: */
108                           if (extendedA) {              if (NULL != A_p)
110                              for (s=0;s<p.row_numShapes;s++) {                  add_EM_S = TRUE;
111                                for (r=0;r<p.col_numShapes;r++) {                  if (extendedA)
112                                  rtmp=0;                  {
113                                  for (q=0;q<p.numQuadTotal;q++) {                  A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuadTotal)]);
114                                        rtmp+=vol*( DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(0,0,q,DIM,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]                  for (s = 0; s < p.row_numShapes; s++)
115                                                     + DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(0,1,q,DIM,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]                  {
116                                                     + DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(0,2,q,DIM,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]                      for (r = 0; r < p.col_numShapes; r++)
117                                                     + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(1,0,q,DIM,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]                      {
118                                                     + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(1,1,q,DIM,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]                      rtmp = 0;
119                                                     + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(1,2,q,DIM,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]                      for (q = 0; q < p.numQuadTotal; q++)
120                                                     + DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(2,0,q,DIM,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]                      {
121                                                     + DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(2,1,q,DIM,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]                          rtmp +=
122                                                     + DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(2,2,q,DIM,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]);                          vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
123                                  }                                 A_q[INDEX3(0, 0, q, DIM, DIM)] *
124                                  EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                                 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
125                                }                                 DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
126                              }                                 A_q[INDEX3(0, 1, q, DIM, DIM)] *
127                           } else {                                 DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
128                              for (s=0;s<p.row_numShapes;s++) {                                 DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
129                                for (r=0;r<p.col_numShapes;r++) {                                 A_q[INDEX3(0, 2, q, DIM, DIM)] *
130                                         rtmp00=0;                                 DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)] +
131                                         rtmp01=0;                                 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
132                                         rtmp02=0;                                 A_q[INDEX3(1, 0, q, DIM, DIM)] *
133                                         rtmp10=0;                                 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
134                                         rtmp11=0;                                 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
135                                         rtmp12=0;                                 A_q[INDEX3(1, 1, q, DIM, DIM)] *
136                                         rtmp20=0;                                 DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
137                                         rtmp21=0;                                 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
138                                         rtmp22=0;                                 A_q[INDEX3(1, 2, q, DIM, DIM)] *
139                                         for (q=0;q<p.numQuadTotal;q++) {                                 DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)] +
140                                         DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
141                                             rtmp0=vol*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                                 A_q[INDEX3(2, 0, q, DIM, DIM)] *
142                                             rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                                 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
143                                             rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                                 DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
144                                             rtmp02+=rtmp0*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];                                 A_q[INDEX3(2, 1, q, DIM, DIM)] *
145                                         DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
146                                             rtmp1=vol*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];                                 DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
147                                             rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                                 A_q[INDEX3(2, 2, q, DIM, DIM)] *
148                                             rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                                 DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)]);
149                                             rtmp12+=rtmp1*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];                      }
150                              EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
151                                             rtmp2=vol*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];                      }
152                                             rtmp20+=rtmp2*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                  }
153                                             rtmp21+=rtmp2*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                  } else
154                                             rtmp22+=rtmp2*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];                  {
155                                         }                  for (s = 0; s < p.row_numShapes; s++)
156                                         EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp00*A_p[INDEX2(0,0,DIM)]                  {
157                                                                                           +rtmp01*A_p[INDEX2(0,1,DIM)]                      for (r = 0; r < p.col_numShapes; r++)
158                                                                                           +rtmp02*A_p[INDEX2(0,2,DIM)]                      {
159                                                                                           +rtmp10*A_p[INDEX2(1,0,DIM)]                      rtmp00 = 0;
160                                                                                           +rtmp11*A_p[INDEX2(1,1,DIM)]                      rtmp01 = 0;
161                                                                                           +rtmp12*A_p[INDEX2(1,2,DIM)]                      rtmp02 = 0;
162                                                                                           +rtmp20*A_p[INDEX2(2,0,DIM)]                      rtmp10 = 0;
163                                                                                           +rtmp21*A_p[INDEX2(2,1,DIM)]                      rtmp11 = 0;
164                                                                                           +rtmp22*A_p[INDEX2(2,2,DIM)];                      rtmp12 = 0;
165                                }                      rtmp20 = 0;
166                              }                      rtmp21 = 0;
167                            }                      rtmp22 = 0;
168                        }                      for (q = 0; q < p.numQuadTotal; q++)
169                        /**************************************************************/                      {
170                        /*   process B: */
171                        /**************************************************************/                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
172                        if (NULL!=B_p) {                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
173                           add_EM_S=TRUE;                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
174                           if (extendedB) {                          rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
176                              for (s=0;s<p.row_numShapes;s++) {                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
177                                for (r=0;r<p.col_numShapes;r++) {                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
178                                  rtmp=0;                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
179                                  for (q=0;q<p.numQuadTotal;q++) {                          rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
180                                      rtmp+=vol*S[INDEX2(r,q,p.row_numShapes)]*
181                                                   ( DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*B_q[INDEX2(0,q,DIM)]                          rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
182                                                   + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*B_q[INDEX2(1,q,DIM)]                          rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
183                                                   + DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*B_q[INDEX2(2,q,DIM)]);                          rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
184                                  }                          rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
185                                  EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                      }
186                                }                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
187                              }                          rtmp00 * A_p[INDEX2(0, 0, DIM)] + rtmp01 * A_p[INDEX2(0, 1, DIM)] +
188                           } else {                          rtmp02 * A_p[INDEX2(0, 2, DIM)] + rtmp10 * A_p[INDEX2(1, 0, DIM)] +
189                              for (s=0;s<p.row_numShapes;s++) {                          rtmp11 * A_p[INDEX2(1, 1, DIM)] + rtmp12 * A_p[INDEX2(1, 2, DIM)] +
190                                for (r=0;r<p.col_numShapes;r++) {                          rtmp20 * A_p[INDEX2(2, 0, DIM)] + rtmp21 * A_p[INDEX2(2, 1, DIM)] +
191                                    rtmp0=0;                          rtmp22 * A_p[INDEX2(2, 2, DIM)];
192                                    rtmp1=0;                      }
193                                    rtmp2=0;                  }
195                                        rtmp=vol*S[INDEX2(r,q,p.row_numShapes)];              }
196                                        rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                /**************************************************************/
197                                        rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];              /*   process B: */
198                                        rtmp2+=rtmp*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];                /**************************************************************/
199                                    }              if (NULL != B_p)
200                                    EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp0*B_p[0]+rtmp1*B_p[1]+rtmp2*B_p[2];              {
202                              }                  if (extendedB)
203                           }                  {
204                        }                  B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
205                        /**************************************************************/                  for (s = 0; s < p.row_numShapes; s++)
206                        /*   process C: */                  {
207                        /**************************************************************/                      for (r = 0; r < p.col_numShapes; r++)
208                        if (NULL!=C_p) {                      {
210                          if (extendedC) {                      for (q = 0; q < p.numQuadTotal; q++)
212                              for (s=0;s<p.row_numShapes;s++) {                          rtmp += vol * S[INDEX2(r, q, p.row_numShapes)] *
213                                for (r=0;r<p.col_numShapes;r++) {                          (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
214                                  rtmp=0;                           B_q[INDEX2(0, q, DIM)] +
215                                  for (q=0;q<p.numQuadTotal;q++) {                           DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
216                                     rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*                           B_q[INDEX2(1, q, DIM)] +
217                                           ( C_q[INDEX2(0,q,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]                           DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
218                                           + C_q[INDEX2(1,q,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]                           B_q[INDEX2(2, q, DIM)]);
219                                           + C_q[INDEX2(2,q,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]);                      }
220                                  }                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
221                                  EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                      }
222                                }                  }
223                              }                  } else
224                          } else {                  {
225                              for (s=0;s<p.row_numShapes;s++) {                  for (s = 0; s < p.row_numShapes; s++)
226                                for (r=0;r<p.col_numShapes;r++) {                  {
227                                    rtmp0=0;                      for (r = 0; r < p.col_numShapes; r++)
228                                    rtmp1=0;                      {
229                                    rtmp2=0;                      rtmp0 = 0;
230                                    for (q=0;q<p.numQuadTotal;q++) {                      rtmp1 = 0;
231                                        rtmp=vol*S[INDEX2(s,q,p.row_numShapes)];                      rtmp2 = 0;
232                                        rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                      for (q = 0; q < p.numQuadTotal; q++)
233                                        rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                      {
234                                        rtmp2+=rtmp*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];                          rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];
235                                    }                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
236                                    EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp0*C_p[0]+rtmp1*C_p[1]+rtmp2*C_p[2];                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
237                                }                          rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
238                              }                      }
239                          }                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
240                        }                          rtmp0 * B_p[0] + rtmp1 * B_p[1] + rtmp2 * B_p[2];
241                        /************************************************************* */                      }
242                        /* process D */                  }
243                        /**************************************************************/                  }
244                        if (NULL!=D_p) {              }
246                          if (extendedD) {              /*   process C: */
248                              for (s=0;s<p.row_numShapes;s++) {              if (NULL != C_p)
249                                for (r=0;r<p.col_numShapes;r++) {              {
251                                  for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*D_q[q]*S[INDEX2(r,q,p.row_numShapes)];                  if (extendedC)
252                                  EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                  {
253                                }                  C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
254                              }                  for (s = 0; s < p.row_numShapes; s++)
255                          } else {                  {
256                              for (s=0;s<p.row_numShapes;s++) {                      for (r = 0; r < p.col_numShapes; r++)
257                                for (r=0;r<p.col_numShapes;r++) {                      {
258                                    rtmp=0;                      rtmp = 0;
259                                    for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(r,q,p.row_numShapes)];                      for (q = 0; q < p.numQuadTotal; q++)
260                                    EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp*D_p[0];                      {
261                                }                          rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] *
262                              }                          (C_q[INDEX2(0, q, DIM)] *
263                          }                           DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
264                        }                           C_q[INDEX2(1, q, DIM)] *
265                        /**************************************************************/                           DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
266                        /*   process X: */                           C_q[INDEX2(2, q, DIM)] *
267                        /**************************************************************/                           DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)]);
268                        if (NULL!=X_p) {                      }
269                          add_EM_F=TRUE;                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
270                          if (extendedX) {                      }
272                             for (s=0;s<p.row_numShapes;s++) {                  } else
273                                rtmp=0;                  {
274                                for (q=0;q<p.numQuadTotal;q++) {                  for (s = 0; s < p.row_numShapes; s++)
275                                     rtmp+=vol*( DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*X_q[INDEX2(0,q,DIM)]                  {
276                                                  + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*X_q[INDEX2(1,q,DIM)]                      for (r = 0; r < p.col_numShapes; r++)
277                                                  + DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*X_q[INDEX2(2,q,DIM)]);                      {
278                                }                      rtmp0 = 0;
279                                EM_F[INDEX2(0,s,p.numEqu)]+=rtmp;                      rtmp1 = 0;
280                             }                      rtmp2 = 0;
281                           } else {                      for (q = 0; q < p.numQuadTotal; q++)
282                             for (s=0;s<p.row_numShapes;s++) {                      {
283                                rtmp0=0;                          rtmp = vol * S[INDEX2(s, q, p.row_numShapes)];
284                                rtmp1=0;                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
285                                rtmp2=0;                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
286                                for (q=0;q<p.numQuadTotal;q++) {                          rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
287                                     rtmp0+=vol*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                      }
288                                     rtmp1+=vol*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
289                                     rtmp2+=vol*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];                          rtmp0 * C_p[0] + rtmp1 * C_p[1] + rtmp2 * C_p[2];
290                                }                      }
291                                EM_F[INDEX2(0,s,p.numEqu)]+=rtmp0*X_p[0]+rtmp1*X_p[1]+rtmp2*X_p[2];                  }
292                             }                  }
293                          }              }
294                       }                /************************************************************* */
295                       /**************************************************************/              /* process D */
296                       /*   process Y: */                /**************************************************************/
297                       /**************************************************************/              if (NULL != D_p)
298                        if (NULL!=Y_p) {              {
300                          if (extendedY) {                  if (extendedD)
302                             for (s=0;s<p.row_numShapes;s++) {                  D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);
303                                rtmp=0;                  for (s = 0; s < p.row_numShapes; s++)
305                                EM_F[INDEX2(0,s,p.numEqu)]+=rtmp;                      for (r = 0; r < p.col_numShapes; r++)
306                             }                      {
307                           } else {                      rtmp = 0;
308                             for (s=0;s<p.row_numShapes;s++) {                      for (q = 0; q < p.numQuadTotal; q++)
309                                 rtmp=0;                          rtmp +=
310                                 for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)];                          vol * S[INDEX2(s, q, p.row_numShapes)] * D_q[q] *
311                                 EM_F[INDEX2(0,s,p.numEqu)]+=rtmp*Y_p[0];                          S[INDEX2(r, q, p.row_numShapes)];
312                             }                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
313                           }                      }
314                         }                  }
315                         /***********************************************************************************************/                  } else
316                         /* add the element matrices onto the matrix and right hand side                                */                  {
317                         /***********************************************************************************************/                  for (s = 0; s < p.row_numShapes; s++)
318                    {
319                         for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];                      for (r = 0; r < p.col_numShapes; r++)
320                        {
323                            rtmp +=
324                  } /* end color check */                          vol * S[INDEX2(s, q, p.row_numShapes)] *
325               } /* end element loop */                          S[INDEX2(r, q, p.row_numShapes)];
326           } /* end color loop */                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
327                                      rtmp * D_p[0];
331                }
332                  /**************************************************************/
333                /*   process X: */
334                  /**************************************************************/
335                if (NULL != X_p)
336                {
338                    if (extendedX)
339                    {
340                    X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
341                    for (s = 0; s < p.row_numShapes; s++)
342                    {
343                        rtmp = 0;
344                        for (q = 0; q < p.numQuadTotal; q++)
345                        {
346                        rtmp +=
347                            vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
348                               X_q[INDEX2(0, q, DIM)] +
349                               DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
350                               X_q[INDEX2(1, q, DIM)] +
351                               DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
352                               X_q[INDEX2(2, q, DIM)]);
353                        }
354                        EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
355                    }
356                    } else
357                    {
358                    for (s = 0; s < p.row_numShapes; s++)
359                    {
360                        rtmp0 = 0;
361                        rtmp1 = 0;
362                        rtmp2 = 0;
363                        for (q = 0; q < p.numQuadTotal; q++)
364                        {
365                        rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
366                        rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
367                        rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
368                        }
369                        EM_F[INDEX2(0, s, p.numEqu)] += rtmp0 * X_p[0] + rtmp1 * X_p[1] + rtmp2 * X_p[2];
370                    }
371                    }
372                }
373                 /**************************************************************/
374                /*   process Y: */
375                 /**************************************************************/
376                if (NULL != Y_p)
377                {
379                    if (extendedY)
380                    {
381                    Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);
382                    for (s = 0; s < p.row_numShapes; s++)
383                    {
384                        rtmp = 0;
385                        for (q = 0; q < p.numQuadTotal; q++)
386                        rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[q];
387                        EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
388                    }
389                    } else
390                    {
391                    for (s = 0; s < p.row_numShapes; s++)
392                    {
393                        rtmp = 0;
394                        for (q = 0; q < p.numQuadTotal; q++)
395                        rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
396                        EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
397                    }
398                    }
399                }
400                   /***********************************************************************************************/
401                /* add the element matrices onto the matrix and right hand side                                */
402                   /***********************************************************************************************/
403
404                for (q = 0; q < p.row_numShapesTotal; q++)
405                    row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
406
408                    Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
409                               p.row_DOF_UpperBound);
412                                      p.col_numShapesTotal, row_index, p.numComp, EM_S);
413
414                }       /* end color check */
415            }       /* end element loop */
416            }           /* end color loop */
417