/[escript]/branches/doubleplusgood/dudley/src/Assemble_PDE_System2_3D.cpp
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3203 by jfenwick, Thu Sep 23 23:51:26 2010 UTC revision 3204 by jfenwick, Thu Sep 23 23:59:39 2010 UTC
# Line 65  void Dudley_Assemble_PDE_System2_3D(Asse Line 65  void Dudley_Assemble_PDE_System2_3D(Asse
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      const double* S=p.shapeFns;      const double* S=p.shapeFns;
68      dim_t len_EM_S = p.row_numShapesTotal * p.row_numShapesTotal * p.numEqu * p.numComp;      dim_t len_EM_S = p.numShapes * p.numShapes * p.numEqu * p.numComp;
69      dim_t len_EM_F = p.row_numShapesTotal * p.numEqu;      dim_t len_EM_F = p.numShapes * p.numEqu;
70    
71  #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)
72      {      {
73      EM_S = THREAD_MEMALLOC(len_EM_S, double);      EM_S = THREAD_MEMALLOC(len_EM_S, double);
74      EM_F = THREAD_MEMALLOC(len_EM_F, double);      EM_F = THREAD_MEMALLOC(len_EM_F, double);
75      row_index = THREAD_MEMALLOC(p.row_numShapesTotal, index_t);      row_index = THREAD_MEMALLOC(p.numShapes, index_t);
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      {      {
# Line 94  void Dudley_Assemble_PDE_System2_3D(Asse Line 94  void Dudley_Assemble_PDE_System2_3D(Asse
94              Y_p = getSampleDataRO(Y, e);              Y_p = getSampleDataRO(Y, e);
95              double vol = p.row_jac->absD[e] * p.row_jac->quadweight;              double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
96    
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.numShapes, DIM, p.numQuadTotal, 1)]);
98              for (q = 0; q < len_EM_S; ++q)              for (q = 0; q < len_EM_S; ++q)
99                  EM_S[q] = 0;                  EM_S[q] = 0;
100              for (q = 0; q < len_EM_F; ++q)              for (q = 0; q < len_EM_F; ++q)
# Line 111  void Dudley_Assemble_PDE_System2_3D(Asse Line 111  void Dudley_Assemble_PDE_System2_3D(Asse
111                  if (extendedA)                  if (extendedA)
112                  {                  {
113                  A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuadTotal)]);                  A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuadTotal)]);
114                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
115                  {                  {
116                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
117                      {                      {
118                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
119                      {                      {
# Line 123  void Dudley_Assemble_PDE_System2_3D(Asse Line 123  void Dudley_Assemble_PDE_System2_3D(Asse
123                          for (q = 0; q < p.numQuadTotal; q++)                          for (q = 0; q < p.numQuadTotal; q++)
124                          {                          {
125                              rtmp +=                              rtmp +=
126                              vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *                              vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
127                                     A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]                                     A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
128                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +                                     * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
129                                     DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *                                     DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
130                                     A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]                                     A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
131                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +                                     * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
132                                     DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *                                     DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
133                                     A_q[INDEX5(k, 0, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]                                     A_q[INDEX5(k, 0, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
134                                     * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)] +                                     * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
135                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *                                     DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
136                                     A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]                                     A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
137                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +                                     * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
138                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *                                     DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
139                                     A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]                                     A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
140                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +                                     * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
141                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *                                     DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
142                                     A_q[INDEX5(k, 1, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]                                     A_q[INDEX5(k, 1, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
143                                     * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)] +                                     * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
144                                     DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *                                     DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
145                                     A_q[INDEX5(k, 2, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]                                     A_q[INDEX5(k, 2, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
146                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +                                     * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
147                                     DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *                                     DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
148                                     A_q[INDEX5(k, 2, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]                                     A_q[INDEX5(k, 2, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
149                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +                                     * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
150                                     DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *                                     DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
151                                     A_q[INDEX5(k, 2, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]                                     A_q[INDEX5(k, 2, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
152                                     * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)]);                                     * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
153    
154                          }                          }
155                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
156                              rtmp;                              rtmp;
157                          }                          }
158                      }                      }
# Line 160  void Dudley_Assemble_PDE_System2_3D(Asse Line 160  void Dudley_Assemble_PDE_System2_3D(Asse
160                  }                  }
161                  } else                  } else
162                  {                  {
163                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
164                  {                  {
165                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
166                      {                      {
167                      rtmp00 = 0;                      rtmp00 = 0;
168                      rtmp01 = 0;                      rtmp01 = 0;
# Line 175  void Dudley_Assemble_PDE_System2_3D(Asse Line 175  void Dudley_Assemble_PDE_System2_3D(Asse
175                      rtmp22 = 0;                      rtmp22 = 0;
176                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
177                      {                      {
178                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
179                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
180                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
181                          rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];                          rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
182    
183                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
184                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
185                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
186                          rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];                          rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
187    
188                          rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];                          rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
189                          rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];                          rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
190                          rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];                          rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
191                          rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];                          rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
192                      }                      }
193                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
194                      {                      {
195                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
196                          {                          {
197                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
198                              rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)] +                              rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)] +
199                              rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)] +                              rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)] +
200                              rtmp02 * A_p[INDEX4(k, 0, m, 2, p.numEqu, DIM, p.numComp)] +                              rtmp02 * A_p[INDEX4(k, 0, m, 2, p.numEqu, DIM, p.numComp)] +
# Line 219  void Dudley_Assemble_PDE_System2_3D(Asse Line 219  void Dudley_Assemble_PDE_System2_3D(Asse
219                  if (extendedB)                  if (extendedB)
220                  {                  {
221                  B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);                  B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);
222                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
223                  {                  {
224                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
225                      {                      {
226                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
227                      {                      {
# Line 231  void Dudley_Assemble_PDE_System2_3D(Asse Line 231  void Dudley_Assemble_PDE_System2_3D(Asse
231                          for (q = 0; q < p.numQuadTotal; q++)                          for (q = 0; q < p.numQuadTotal; q++)
232                          {                          {
233                              rtmp +=                              rtmp +=
234                              vol * S[INDEX2(r, q, p.row_numShapes)] *                              vol * S[INDEX2(r, q, p.numShapes)] *
235                              (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *                              (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
236                               B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +                               B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
237                               DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *                               DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
238                               B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)] +                               B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)] +
239                               DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *                               DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
240                               B_q[INDEX4(k, 2, m, q, p.numEqu, DIM, p.numComp)]);                               B_q[INDEX4(k, 2, m, q, p.numEqu, DIM, p.numComp)]);
241                          }                          }
242                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
243                              rtmp;                              rtmp;
244                          }                          }
245                      }                      }
# Line 247  void Dudley_Assemble_PDE_System2_3D(Asse Line 247  void Dudley_Assemble_PDE_System2_3D(Asse
247                  }                  }
248                  } else                  } else
249                  {                  {
250                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
251                  {                  {
252                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
253                      {                      {
254                      rtmp0 = 0;                      rtmp0 = 0;
255                      rtmp1 = 0;                      rtmp1 = 0;
256                      rtmp2 = 0;                      rtmp2 = 0;
257                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
258                      {                      {
259                          rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];                          rtmp = vol * S[INDEX2(r, q, p.numShapes)];
260                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
261                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
262                          rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];                          rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
263                      }                      }
264                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
265                      {                      {
266                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
267                          {                          {
268                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
269                              rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +                              rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
270                              rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)] +                              rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)] +
271                              rtmp2 * B_p[INDEX3(k, 2, m, p.numEqu, DIM)];                              rtmp2 * B_p[INDEX3(k, 2, m, p.numEqu, DIM)];
# Line 284  void Dudley_Assemble_PDE_System2_3D(Asse Line 284  void Dudley_Assemble_PDE_System2_3D(Asse
284                  if (extendedC)                  if (extendedC)
285                  {                  {
286                  C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);                  C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);
287                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
288                  {                  {
289                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
290                      {                      {
291                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
292                      {                      {
# Line 296  void Dudley_Assemble_PDE_System2_3D(Asse Line 296  void Dudley_Assemble_PDE_System2_3D(Asse
296                          for (q = 0; q < p.numQuadTotal; q++)                          for (q = 0; q < p.numQuadTotal; q++)
297                          {                          {
298                              rtmp +=                              rtmp +=
299                              vol * S[INDEX2(s, q, p.row_numShapes)] *                              vol * S[INDEX2(s, q, p.numShapes)] *
300                              (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *                              (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
301                               DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +                               DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
302                               C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *                               C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
303                               DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +                               DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
304                               C_q[INDEX4(k, m, 2, q, p.numEqu, p.numComp, DIM)] *                               C_q[INDEX4(k, m, 2, q, p.numEqu, p.numComp, DIM)] *
305                               DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)]);                               DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
306                          }                          }
307                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
308                              rtmp;                              rtmp;
309                          }                          }
310                      }                      }
# Line 312  void Dudley_Assemble_PDE_System2_3D(Asse Line 312  void Dudley_Assemble_PDE_System2_3D(Asse
312                  }                  }
313                  } else                  } else
314                  {                  {
315                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
316                  {                  {
317                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
318                      {                      {
319                      rtmp0 = 0;                      rtmp0 = 0;
320                      rtmp1 = 0;                      rtmp1 = 0;
321                      rtmp2 = 0;                      rtmp2 = 0;
322                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
323                      {                      {
324                          rtmp = vol * S[INDEX2(s, q, p.row_numShapes)];                          rtmp = vol * S[INDEX2(s, q, p.numShapes)];
325                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
326                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
327                          rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];                          rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
328                      }                      }
329                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
330                      {                      {
331                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
332                          {                          {
333                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
334                              rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +                              rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
335                              rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)] +                              rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)] +
336                              rtmp2 * C_p[INDEX3(k, m, 2, p.numEqu, p.numComp)];                              rtmp2 * C_p[INDEX3(k, m, 2, p.numEqu, p.numComp)];
# Line 349  void Dudley_Assemble_PDE_System2_3D(Asse Line 349  void Dudley_Assemble_PDE_System2_3D(Asse
349                  if (extendedD)                  if (extendedD)
350                  {                  {
351                  D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);                  D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);
352                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
353                  {                  {
354                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
355                      {                      {
356                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
357                      {                      {
# Line 361  void Dudley_Assemble_PDE_System2_3D(Asse Line 361  void Dudley_Assemble_PDE_System2_3D(Asse
361                          for (q = 0; q < p.numQuadTotal; q++)                          for (q = 0; q < p.numQuadTotal; q++)
362                          {                          {
363                              rtmp +=                              rtmp +=
364                              vol * S[INDEX2(s, q, p.row_numShapes)] *                              vol * S[INDEX2(s, q, p.numShapes)] *
365                              D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *                              D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
366                              S[INDEX2(r, q, p.row_numShapes)];                              S[INDEX2(r, q, p.numShapes)];
367                          }                          }
368                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
369                              rtmp;                              rtmp;
370                          }                          }
371                      }                      }
# Line 373  void Dudley_Assemble_PDE_System2_3D(Asse Line 373  void Dudley_Assemble_PDE_System2_3D(Asse
373                  }                  }
374                  } else                  } else
375                  {                  {
376                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
377                  {                  {
378                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
379                      {                      {
380                      rtmp = 0;                      rtmp = 0;
381                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
382                          rtmp +=                          rtmp +=
383                          vol * S[INDEX2(s, q, p.row_numShapes)] *                          vol * S[INDEX2(s, q, p.numShapes)] *
384                          S[INDEX2(r, q, p.row_numShapes)];                          S[INDEX2(r, q, p.numShapes)];
385                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
386                      {                      {
387                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
388                          {                          {
389                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
390                              rtmp * D_p[INDEX2(k, m, p.numEqu)];                              rtmp * D_p[INDEX2(k, m, p.numEqu)];
391                          }                          }
392                      }                      }
# Line 403  void Dudley_Assemble_PDE_System2_3D(Asse Line 403  void Dudley_Assemble_PDE_System2_3D(Asse
403                  if (extendedX)                  if (extendedX)
404                  {                  {
405                  X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuadTotal)]);                  X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuadTotal)]);
406                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
407                  {                  {
408                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
409                      {                      {
# Line 411  void Dudley_Assemble_PDE_System2_3D(Asse Line 411  void Dudley_Assemble_PDE_System2_3D(Asse
411                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
412                      {                      {
413                          rtmp +=                          rtmp +=
414                          vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *                          vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
415                                 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +                                 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
416                                 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *                                 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
417                                 X_q[INDEX3(k, 1, q, p.numEqu, DIM)] +                                 X_q[INDEX3(k, 1, q, p.numEqu, DIM)] +
418                                 DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *                                 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
419                                 X_q[INDEX3(k, 2, q, p.numEqu, DIM)]);                                 X_q[INDEX3(k, 2, q, p.numEqu, DIM)]);
420                      }                      }
421                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp;                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
# Line 423  void Dudley_Assemble_PDE_System2_3D(Asse Line 423  void Dudley_Assemble_PDE_System2_3D(Asse
423                  }                  }
424                  } else                  } else
425                  {                  {
426                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
427                  {                  {
428                      rtmp0 = 0;                      rtmp0 = 0;
429                      rtmp1 = 0;                      rtmp1 = 0;
430                      rtmp2 = 0;                      rtmp2 = 0;
431                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
432                      {                      {
433                      rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];                      rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
434                      rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];                      rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
435                      rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];                      rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
436                      }                      }
437                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
438                      {                      {
# Line 451  void Dudley_Assemble_PDE_System2_3D(Asse Line 451  void Dudley_Assemble_PDE_System2_3D(Asse
451                  if (extendedY)                  if (extendedY)
452                  {                  {
453                  Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);                  Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);
454                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
455                  {                  {
456                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
457                      {                      {
458                      rtmp = 0.;                      rtmp = 0.;
459                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
460                          rtmp +=                          rtmp +=
461                          vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];                          vol * S[INDEX2(s, q, p.numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
462                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp;                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
463                      }                      }
464                  }                  }
465                  } else                  } else
466                  {                  {
467                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
468                  {                  {
469                      rtmp = 0;                      rtmp = 0;
470                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
471                      rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];                      rtmp += vol * S[INDEX2(s, q, p.numShapes)];
472                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
473                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
474                  }                  }
# Line 478  void Dudley_Assemble_PDE_System2_3D(Asse Line 478  void Dudley_Assemble_PDE_System2_3D(Asse
478                 /***********************************************************************************************/                 /***********************************************************************************************/
479              /* add the element matrices onto the matrix and right hand side                                */              /* add the element matrices onto the matrix and right hand side                                */
480                 /***********************************************************************************************/                 /***********************************************************************************************/
481              for (q = 0; q < p.row_numShapesTotal; q++)              for (q = 0; q < p.numShapes; q++)
482                  row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];                  row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
483    
484              if (add_EM_F)              if (add_EM_F)
485                  Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,                  Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p,
486                             p.row_DOF_UpperBound);                             p.row_DOF_UpperBound);
487              if (add_EM_S)              if (add_EM_S)
488                  Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,                  Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
489                                    p.row_numShapesTotal, row_index, p.numComp, EM_S);                                    p.numShapes, row_index, p.numComp, EM_S);
490              }       /* end color check */              }       /* end color check */
491          }       /* end element loop */          }       /* end element loop */
492          }           /* end color loop */          }           /* end color loop */

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

  ViewVC Help
Powered by ViewVC 1.1.26