Diff of /branches/doubleplusgood/dudley/src/Assemble_PDE_System2_3D.cpp

revision 3186 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 3D domain. The shape functions for test and solution must be identical  */  /*    u has p.numComp components in a 3D domain. The shape functions for test and solution must be identical  */
# Line 32  Line 30
30  /*      X = p.numEqu x 3  */  /*      X = p.numEqu x 3  */
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_3D(Assemble_Parameters p, Dudley_ElementFile* elements,  void Dudley_Assemble_PDE_System2_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,k,m;      register dim_t q, s, r, k, m;
56      register double rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22;      register double rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22;
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, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22,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, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22,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            double vol=p.row_jac->absD[e]*p.row_jac->quadweight;              C_p = getSampleDataRO(C, e);
91                D_p = getSampleDataRO(D, e);
92                        DSDX=&(p.row_jac->DSDX[INDEX5(0,0,0,0,e,p.row_numShapesTotal,DIM,p.numQuadTotal,1)]);              X_p = getSampleDataRO(X, e);
93                        for (q=0;q<len_EM_S;++q) EM_S[q]=0;              Y_p = getSampleDataRO(Y, e);
94                        for (q=0;q<len_EM_F;++q) EM_F[q]=0;              double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
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;
101                        if (NULL!=A_p) {              add_EM_F = FALSE;
103                           if (extendedA) {
105                              for (s=0;s<p.row_numShapes;s++) {              /*   process A: */
106                                for (r=0;r<p.col_numShapes;r++) {                /**************************************************************/
107                                  for (k=0;k<p.numEqu;k++) {              if (NULL != A_p)
108                                    for (m=0;m<p.numComp;m++) {              {
110                                       for (q=0;q<p.numQuadTotal;q++) {                  if (extendedA)
111                                          rtmp+=vol*(                  {
112                                            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)]);
113                                           +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++)
114                                           +DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,2,q, p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]                  {
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)]                      for (r = 0; r < p.col_numShapes; r++)
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)]                      {
117                                           +DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,1,m,2,q, p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]                      for (k = 0; k < p.numEqu; k++)
118                                           +DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,2,m,0,q, p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]                      {
119                                           +DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,2,m,1,q, p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]                          for (m = 0; m < p.numComp; m++)
120                                           +DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,2,m,2,q, p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]);                          {
121                                                        rtmp = 0;
122                                       }                          for (q = 0; q < p.numQuadTotal; q++)
123                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                          {
124                                    }                              rtmp +=
125                                  }                              vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
126                                }                                     A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
127                              }                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
128                           } else {                                     DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
129                              for (s=0;s<p.row_numShapes;s++) {                                     A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
130                                for (r=0;r<p.col_numShapes;r++) {                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
131                                    rtmp00=0;                                     DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
132                                    rtmp01=0;                                     A_q[INDEX5(k, 0, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
133                                    rtmp02=0;                                     * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)] +
134                                    rtmp10=0;                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
135                                    rtmp11=0;                                     A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
136                                    rtmp12=0;                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
137                                    rtmp20=0;                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
138                                    rtmp21=0;                                     A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
139                                    rtmp22=0;                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
140                                    for (q=0;q<p.numQuadTotal;q++) {                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
141                                          rtmp0=vol*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                                     A_q[INDEX5(k, 1, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
142                                          rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                                     * DSDX[INDEX3(r, 2, 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[INDEX5(k, 2, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
145                                             * DSDX[INDEX3(r, 0, 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[INDEX5(k, 2, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
148                                          rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
149                                          rtmp12+=rtmp1*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];                                     DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
150                                             A_q[INDEX5(k, 2, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
151                                          rtmp2=vol*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];                                     * DSDX[INDEX3(r, 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)];                          }
154                                          rtmp22+=rtmp2*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
155                                    }                              rtmp;
156                                    for (k=0;k<p.numEqu;k++) {                          }
157                                       for (m=0;m<p.numComp;m++) {                      }
158                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+= rtmp00*A_p[INDEX4(k,0,m,0, p.numEqu,DIM,p.numComp)]                      }
159                                                                                             +rtmp01*A_p[INDEX4(k,0,m,1, p.numEqu,DIM,p.numComp)]                  }
160                                                                                             +rtmp02*A_p[INDEX4(k,0,m,2, p.numEqu,DIM,p.numComp)]                  } else
161                                                                                             +rtmp10*A_p[INDEX4(k,1,m,0, p.numEqu,DIM,p.numComp)]                  {
162                                                                                             +rtmp11*A_p[INDEX4(k,1,m,1, p.numEqu,DIM,p.numComp)]                  for (s = 0; s < p.row_numShapes; s++)
163                                                                                             +rtmp12*A_p[INDEX4(k,1,m,2, p.numEqu,DIM,p.numComp)]                  {
164                                                                                             +rtmp20*A_p[INDEX4(k,2,m,0, p.numEqu,DIM,p.numComp)]                      for (r = 0; r < p.col_numShapes; r++)
165                                                                                             +rtmp21*A_p[INDEX4(k,2,m,1, p.numEqu,DIM,p.numComp)]                      {
166                                                                                             +rtmp22*A_p[INDEX4(k,2,m,2, p.numEqu,DIM,p.numComp)];                      rtmp00 = 0;
167                                       }                      rtmp01 = 0;
168                                    }                      rtmp02 = 0;
169                                }                      rtmp10 = 0;
170                              }                      rtmp11 = 0;
171                            }                      rtmp12 = 0;
172                        }                      rtmp20 = 0;
173                        /**************************************************************/                      rtmp21 = 0;
174                        /*   process B: */                      rtmp22 = 0;
175                        /**************************************************************/                      for (q = 0; q < p.numQuadTotal; q++)
176                        if (NULL!=B_p) {                      {
177                           add_EM_S=TRUE;                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
178                           if (extendedB) {                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
179                      B_q=&(B_p[INDEX5(0,0,0,0,0, p.numEqu,DIM,p.numComp,p.numQuadTotal)]);                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
180                              for (s=0;s<p.row_numShapes;s++) {                          rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
181                                for (r=0;r<p.col_numShapes;r++) {
182                                  for (k=0;k<p.numEqu;k++) {                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
183                                    for (m=0;m<p.numComp;m++) {                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
184                                       rtmp=0;                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
185                                       for (q=0;q<p.numQuadTotal;q++) {                          rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
186                                           rtmp+=vol*S[INDEX2(r,q,p.row_numShapes)]*(
187                                                   DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]                          rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
188                                                 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]                          rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
189                                                 + DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,2,m,q,p.numEqu,DIM,p.numComp)] );                          rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
190                                       }                          rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
191                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                      }
192                                    }                      for (k = 0; k < p.numEqu; k++)
193                                  }                      {
194                                }                          for (m = 0; m < p.numComp; m++)
195                              }                          {
196                           } else {                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
197                              for (s=0;s<p.row_numShapes;s++) {                              rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)] +
198                                for (r=0;r<p.col_numShapes;r++) {                              rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)] +
199                                       rtmp0=0;                              rtmp02 * A_p[INDEX4(k, 0, m, 2, p.numEqu, DIM, p.numComp)] +
200                                       rtmp1=0;                              rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)] +
201                                       rtmp2=0;                              rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)] +
202                                       for (q=0;q<p.numQuadTotal;q++) {                              rtmp12 * A_p[INDEX4(k, 1, m, 2, p.numEqu, DIM, p.numComp)] +
203                                            rtmp=vol*S[INDEX2(r,q,p.row_numShapes)];                              rtmp20 * A_p[INDEX4(k, 2, m, 0, p.numEqu, DIM, p.numComp)] +
204                                            rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                              rtmp21 * A_p[INDEX4(k, 2, m, 1, p.numEqu, DIM, p.numComp)] +
205                                            rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];                              rtmp22 * A_p[INDEX4(k, 2, m, 2, p.numEqu, DIM, p.numComp)];
206                                            rtmp2+=rtmp*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];                          }
207                                       }                      }
208                                       for (k=0;k<p.numEqu;k++) {                      }
209                                          for (m=0;m<p.numComp;m++) {                  }
210                                             EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]                  }
211                                                                                               +rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)]              }
212                                                                                               +rtmp2*B_p[INDEX3(k,2,m,p.numEqu,DIM)];                /**************************************************************/
213                                          }              /*   process B: */
214                                       }                /**************************************************************/
215                                }              if (NULL != B_p)
216                              }              {
218                        }                  if (extendedB)
219                        /**************************************************************/                  {
220                        /*   process C: */                  B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);
221                        /**************************************************************/                  for (s = 0; s < p.row_numShapes; s++)
222                        if (NULL!=C_p) {                  {
223                          add_EM_S=TRUE;                      for (r = 0; r < p.col_numShapes; r++)
224                          if (extendedC) {                      {
225                  C_q=&(C_p[INDEX5(0,0,0,0,0 ,p.numEqu,p.numComp,DIM, p.numQuadTotal)]);                      for (k = 0; k < p.numEqu; k++)
226                              for (s=0;s<p.row_numShapes;s++) {                      {
227                                for (r=0;r<p.col_numShapes;r++) {                          for (m = 0; m < p.numComp; m++)
228                                  for (k=0;k<p.numEqu;k++) {                          {
229                                    for (m=0;m<p.numComp;m++) {                          rtmp = 0;
230                                      rtmp=0;                          for (q = 0; q < p.numQuadTotal; q++)
232                                           rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*(                              rtmp +=
233                                                      C_q[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]                              vol * S[INDEX2(r, q, p.row_numShapes)] *
234                                                     +C_q[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]                              (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
235                                                     +C_q[INDEX4(k,m,2,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]);                               B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
236                                      }                               DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
237                                      EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                               B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)] +
238                                    }                               DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
239                                  }                               B_q[INDEX4(k, 2, m, q, p.numEqu, DIM, p.numComp)]);
240                                }                          }
241                              }                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
242                          } else {                              rtmp;
243                              for (s=0;s<p.row_numShapes;s++) {                          }
244                                for (r=0;r<p.col_numShapes;r++) {                      }
245                                       rtmp0=0;                      }
246                                       rtmp1=0;                  }
247                                       rtmp2=0;                  } else
249                                          rtmp=vol*S[INDEX2(s,q,p.row_numShapes)];                  for (s = 0; s < p.row_numShapes; s++)
250                                          rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                  {
251                                          rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                      for (r = 0; r < p.col_numShapes; r++)
252                                          rtmp2+=rtmp*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];                      {
253                                       }                      rtmp0 = 0;
254                                       for (k=0;k<p.numEqu;k++) {                      rtmp1 = 0;
255                                          for (m=0;m<p.numComp;m++) {                      rtmp2 = 0;
256                                                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)]                      for (q = 0; q < p.numQuadTotal; q++)
257                                                                                                  +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)]                      {
258                                                                                                  +rtmp2*C_p[INDEX3(k,m,2,p.numEqu,p.numComp)];                          rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];
259                                           }                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
260                                       }                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
261                                }                          rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
262                              }                      }
263                          }                      for (k = 0; k < p.numEqu; k++)
264                        }                      {
265                        /************************************************************* */                          for (m = 0; m < p.numComp; m++)
266                        /* process D */                          {
267                        /**************************************************************/                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
268                        if (NULL!=D_p) {                              rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
269                          add_EM_S=TRUE;                              rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)] +
270                          if (extendedD) {                              rtmp2 * B_p[INDEX3(k, 2, m, p.numEqu, DIM)];
272                              for (s=0;s<p.row_numShapes;s++) {                      }
273                                for (r=0;r<p.col_numShapes;r++) {                      }
274                                  for (k=0;k<p.numEqu;k++) {                  }
275                                    for (m=0;m<p.numComp;m++) {                  }
276                                      rtmp=0;              }
278                                          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)];              /*   process C: */
279                                      }                /**************************************************************/
280                                      EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;              if (NULL != C_p)
281                                    }              {
283                                }                  if (extendedC)
284                              }                  {
285                          } else {                  C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);
286                              for (s=0;s<p.row_numShapes;s++) {                  for (s = 0; s < p.row_numShapes; s++)
287                                for (r=0;r<p.col_numShapes;r++) {                  {
288                                    rtmp=0;                      for (r = 0; r < p.col_numShapes; r++)
290                                    for (k=0;k<p.numEqu;k++) {                      for (k = 0; k < p.numEqu; k++)
291                                        for (m=0;m<p.numComp;m++) {                      {
292                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];                          for (m = 0; m < p.numComp; m++)
293                                       }                          {
294                                    }                          rtmp = 0;
295                                }                          for (q = 0; q < p.numQuadTotal; q++)
296                              }                          {
297                          }                              rtmp +=
298                        }                              vol * S[INDEX2(s, q, p.row_numShapes)] *
299                        /**************************************************************/                              (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
300                        /*   process X: */                               DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
301                        /**************************************************************/                               C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
302                        if (NULL!=X_p) {                               DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
303                          add_EM_F=TRUE;                               C_q[INDEX4(k, m, 2, q, p.numEqu, p.numComp, DIM)] *
304                          if (extendedX) {                               DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)]);
306                             for (s=0;s<p.row_numShapes;s++) {                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
307                                for (k=0;k<p.numEqu;k++) {                              rtmp;
308                                   rtmp=0;                          }
310                                         rtmp+=vol* ( DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,0,q,p.numEqu,DIM)]                      }
311                                                       + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,1,q,p.numEqu,DIM)]                  }
312                                                       + DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,2,q,p.numEqu,DIM)]);                  } else
313                                   }                  {
314                                   EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;                  for (s = 0; s < p.row_numShapes; s++)
315                                }                  {
316                             }                      for (r = 0; r < p.col_numShapes; r++)
317                          } else {                      {
318                             for (s=0;s<p.row_numShapes;s++) {                      rtmp0 = 0;
319                               rtmp0=0;                      rtmp1 = 0;
320                               rtmp1=0;                      rtmp2 = 0;
321                               rtmp2=0;                      for (q = 0; q < p.numQuadTotal; q++)
323                                  rtmp0+=vol*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                          rtmp = vol * S[INDEX2(s, q, p.row_numShapes)];
324                                  rtmp1+=vol*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
325                                  rtmp2+=vol*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
326                               }                          rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
327                               for (k=0;k<p.numEqu;k++) {                      }
328                                         EM_F[INDEX2(k,s,p.numEqu)]+=rtmp0*X_p[INDEX2(k,0,p.numEqu)]                      for (k = 0; k < p.numEqu; k++)
329                                                                    +rtmp1*X_p[INDEX2(k,1,p.numEqu)]                      {
330                                                                    +rtmp2*X_p[INDEX2(k,2,p.numEqu)];                          for (m = 0; m < p.numComp; m++)
331                               }                          {
332                             }                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
333                          }                              rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
334                       }                              rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)] +
335                       /**************************************************************/                              rtmp2 * C_p[INDEX3(k, m, 2, p.numEqu, p.numComp)];
336                       /*   process Y: */                          }
337                       /**************************************************************/                      }
338                        if (NULL!=Y_p) {                      }
340                          if (extendedY) {                  }
342                             for (s=0;s<p.row_numShapes;s++) {                /************************************************************* */
343                                for (k=0;k<p.numEqu;k++) {              /* process D */
344                                   rtmp=0.;                /**************************************************************/
345                                   for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*Y_q[INDEX2(k,q,p.numEqu)];              if (NULL != D_p)
346                                   EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;              {
348                             }                  if (extendedD)
349                           } else {                  {
350                             for (s=0;s<p.row_numShapes;s++) {                  D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);
351                                 rtmp=0;                  for (s = 0; s < p.row_numShapes; s++)
353                                 for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];                      for (r = 0; r < p.col_numShapes; r++)
354                             }                      {
355                           }                      for (k = 0; k < p.numEqu; k++)
356                         }                      {
357                            for (m = 0; m < p.numComp; m++)
358                         /***********************************************************************************************/                          {
359                         /* add the element matrices onto the matrix and right hand side                                */                                    rtmp = 0;
360                         /***********************************************************************************************/                          for (q = 0; q < p.numQuadTotal; q++)
361                         for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];                          {
362                                                      rtmp +=
365                  } /* end color check */                              S[INDEX2(r, q, p.row_numShapes)];
366               } /* end element loop */                          }
367           } /* end color loop */                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
368                                          rtmp;
372                    }
373                    } else
374                    {
375                    for (s = 0; s < p.row_numShapes; s++)
376                    {
377                        for (r = 0; r < p.col_numShapes; r++)
378                        {
379                        rtmp = 0;
380                        for (q = 0; q < p.numQuadTotal; q++)
381                            rtmp +=
382                            vol * S[INDEX2(s, q, p.row_numShapes)] *
383                            S[INDEX2(r, q, p.row_numShapes)];
384                        for (k = 0; k < p.numEqu; k++)
385                        {
386                            for (m = 0; m < p.numComp; m++)
387                            {
388                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
389                                rtmp * D_p[INDEX2(k, m, p.numEqu)];
390                            }
391                        }
392                        }
393                    }
394                    }
395                }
396                  /**************************************************************/
397                /*   process X: */
398                  /**************************************************************/
399                if (NULL != X_p)
400                {
402                    if (extendedX)
403                    {
404                    X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuadTotal)]);
405                    for (s = 0; s < p.row_numShapes; s++)
406                    {
407                        for (k = 0; k < p.numEqu; k++)
408                        {
409                        rtmp = 0;
410                        for (q = 0; q < p.numQuadTotal; q++)
411                        {
412                            rtmp +=
413                            vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
414                                   X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
415                                   DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
416                                   X_q[INDEX3(k, 1, q, p.numEqu, DIM)] +
417                                   DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
418                                   X_q[INDEX3(k, 2, q, p.numEqu, DIM)]);
419                        }
420                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
421                        }
422                    }
423                    } else
424                    {
425                    for (s = 0; s < p.row_numShapes; s++)
426                    {
427                        rtmp0 = 0;
428                        rtmp1 = 0;
429                        rtmp2 = 0;
430                        for (q = 0; q < p.numQuadTotal; q++)
431                        {
432                        rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
433                        rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
434                        rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
435                        }
436                        for (k = 0; k < p.numEqu; k++)
437                        {
438                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp0 * X_p[INDEX2(k, 0, p.numEqu)]
439                            + rtmp1 * X_p[INDEX2(k, 1, p.numEqu)] + rtmp2 * X_p[INDEX2(k, 2, p.numEqu)];
440                        }
441                    }
442                    }
443                }
444                 /**************************************************************/
445                /*   process Y: */
446                 /**************************************************************/
447                if (NULL != Y_p)
448                {
450                    if (extendedY)
451                    {
452                    Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);
453                    for (s = 0; s < p.row_numShapes; s++)
454                    {
455                        for (k = 0; k < p.numEqu; k++)
456                        {
457                        rtmp = 0.;
458                        for (q = 0; q < p.numQuadTotal; q++)
459                            rtmp +=
460                            vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
461                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
462                        }
463                    }
464                    } else
465                    {
466                    for (s = 0; s < p.row_numShapes; s++)
467                    {
468                        rtmp = 0;
469                        for (q = 0; q < p.numQuadTotal; q++)
470                        rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
471                        for (k = 0; k < p.numEqu; k++)
472                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
473                    }
474                    }
475                }
476
477                   /***********************************************************************************************/
478                /* add the element matrices onto the matrix and right hand side                                */
479                   /***********************************************************************************************/
480                for (q = 0; q < p.row_numShapesTotal; q++)
481                    row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
482
484                    Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
485                               p.row_DOF_UpperBound);
488                                      p.col_numShapesTotal, row_index, p.numComp, EM_S);
489                }       /* end color check */
490            }       /* end element loop */
491            }           /* end color loop */
492