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

Diff of /branches/trilinos_from_5897/dudley/src/Assemble_PDE_System2_2D.cpp

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

revision 3187 by jfenwick, Thu Sep 16 02:57:17 2010 UTC revision 3231 by jfenwick, Fri Oct 1 01:53:46 2010 UTC
# Line 63  void Dudley_Assemble_PDE_System2_2D(Asse Line 63  void Dudley_Assemble_PDE_System2_2D(Asse
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;      const double *S = p.shapeFns;
67      dim_t len_EM_S = p.row_numShapesTotal * p.col_numShapesTotal * p.numEqu * p.numComp;      dim_t len_EM_S = p.numShapes * p.numShapes * p.numEqu * p.numComp;
68      dim_t len_EM_F = p.row_numShapesTotal * p.numEqu;      dim_t len_EM_F = p.numShapes * p.numEqu;
69    
70  #pragma omp parallel private(color,EM_S, EM_F, Vol, DSDX, A_p, B_p, C_p, D_p, X_p, Y_p, A_q, B_q, C_q, D_q, X_q, Y_q,row_index,q, s,r,k,m,rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11,add_EM_F, add_EM_S)  #pragma omp parallel private(color,EM_S, EM_F, DSDX, A_p, B_p, C_p, D_p, X_p, Y_p, A_q, B_q, C_q, D_q, X_q, Y_q,row_index,q, s,r,k,m,rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11,add_EM_F, add_EM_S)
71      {      {
72    
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 93  void Dudley_Assemble_PDE_System2_2D(Asse Line 93  void Dudley_Assemble_PDE_System2_2D(Asse
93              X_p = getSampleDataRO(X, e);              X_p = getSampleDataRO(X, e);
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              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.numQuad, 1)]);
97              for (q = 0; q < len_EM_S; ++q)              for (q = 0; q < len_EM_S; ++q)
98                  EM_S[q] = 0;                  EM_S[q] = 0;
99              for (q = 0; q < len_EM_F; ++q)              for (q = 0; q < len_EM_F; ++q)
# Line 110  void Dudley_Assemble_PDE_System2_2D(Asse Line 110  void Dudley_Assemble_PDE_System2_2D(Asse
110                  add_EM_S = TRUE;                  add_EM_S = TRUE;
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.numQuad)]);
114                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
115                  {                  {
116                      for (r = 0; r < p.col_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                      {                      {
120                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
121                          {                          {
122                          rtmp = 0;                          rtmp = 0;
123                          for (q = 0; q < p.numQuadTotal; q++)                          for (q = 0; q < p.numQuad; 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, 1, q, p.row_numShapesTotal, DIM)] *                                     DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
133                                     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)]
134                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +                                     * DSDX[INDEX3(r, 0, 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, 1, q, p.numEqu, DIM, p.numComp, DIM)]                                     A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
137                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);                                     * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
138                          }                          }
139                          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)] += rtmp;
                             rtmp;  
140                          }                          }
141                      }                      }
142                      }                      }
143                  }                  }
144                  } else                  }
145                    else
146                  {                  {
147                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
148                  {                  {
149                      for (r = 0; r < p.col_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
150                      {                      {
151                      rtmp00 = 0;                      rtmp00 = 0;
152                      rtmp01 = 0;                      rtmp01 = 0;
153                      rtmp10 = 0;                      rtmp10 = 0;
154                      rtmp11 = 0;                      rtmp11 = 0;
155                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
156                      {                      {
157                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
158                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
159                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
160                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
161                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
162                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
163                      }                      }
164                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
165                      {                      {
166                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
167                          {                          {
168                          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)] +=
169                              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)]
170                              + 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)]
171                              + rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)]                              + rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)]
# Line 185  void Dudley_Assemble_PDE_System2_2D(Asse Line 185  void Dudley_Assemble_PDE_System2_2D(Asse
185                  add_EM_S = TRUE;                  add_EM_S = TRUE;
186                  if (extendedB)                  if (extendedB)
187                  {                  {
188                  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.numQuad)]);
189                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
190                  {                  {
191                      for (r = 0; r < p.col_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
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                          rtmp = 0;                          rtmp = 0;
198                          for (q = 0; q < p.numQuadTotal; q++)                          for (q = 0; q < p.numQuad; q++)
199                          {                          {
200                              rtmp += vol * S[INDEX2(r, q, p.row_numShapes)] *                              rtmp += vol * S[INDEX2(r, q, p.numShapes)] *
201                              (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *                              (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
202                               B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +                               B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
203                               DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *                               DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
204                               B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);                               B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);
205                          }                          }
206                          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)] += rtmp;
                             rtmp;  
207                          }                          }
208                      }                      }
209                      }                      }
210                  }                  }
211                  } else                  }
212                    else
213                  {                  {
214                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
215                  {                  {
216                      for (r = 0; r < p.col_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
217                      {                      {
218                      rtmp0 = 0;                      rtmp0 = 0;
219                      rtmp1 = 0;                      rtmp1 = 0;
220                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
221                      {                      {
222                          rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];                          rtmp = vol * S[INDEX2(r, q, p.numShapes)];
223                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
224                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
225                      }                      }
226                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
227                      {                      {
228                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
229                          {                          {
230                          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)] +=
231                              rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +                              rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
232                              rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];                              rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];
233                          }                          }
# Line 245  void Dudley_Assemble_PDE_System2_2D(Asse Line 245  void Dudley_Assemble_PDE_System2_2D(Asse
245                  add_EM_S = TRUE;                  add_EM_S = TRUE;
246                  if (extendedC)                  if (extendedC)
247                  {                  {
248                  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.numQuad)]);
249                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
250                  {                  {
251                      for (r = 0; r < p.col_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
252                      {                      {
253                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
254                      {                      {
255                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
256                          {                          {
257                          rtmp = 0;                          rtmp = 0;
258                          for (q = 0; q < p.numQuadTotal; q++)                          for (q = 0; q < p.numQuad; q++)
259                          {                          {
260                              rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] *                              rtmp += vol * S[INDEX2(s, q, p.numShapes)] *
261                              (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *                              (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
262                               DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +                               DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
263                               C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *                               C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
264                               DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);                               DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
265                          }                          }
266                          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)] += rtmp;
                             rtmp;  
267                          }                          }
268                      }                      }
269                      }                      }
270                  }                  }
271                  } else                  }
272                    else
273                  {                  {
274                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
275                  {                  {
276                      for (r = 0; r < p.col_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
277                      {                      {
278                      rtmp0 = 0;                      rtmp0 = 0;
279                      rtmp1 = 0;                      rtmp1 = 0;
280                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
281                      {                      {
282                          rtmp = vol * S[INDEX2(s, q, p.row_numShapes)];                          rtmp = vol * S[INDEX2(s, q, p.numShapes)];
283                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
284                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
285                      }                      }
286                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
287                      {                      {
288                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
289                          {                          {
290                          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)] +=
291                              rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +                              rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
292                              rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];                              rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];
293                          }                          }
# Line 305  void Dudley_Assemble_PDE_System2_2D(Asse Line 305  void Dudley_Assemble_PDE_System2_2D(Asse
305                  add_EM_S = TRUE;                  add_EM_S = TRUE;
306                  if (extendedD)                  if (extendedD)
307                  {                  {
308                  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.numQuad)]);
309                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
310                  {                  {
311                      for (r = 0; r < p.col_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
312                      {                      {
313                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
314                      {                      {
315                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
316                          {                          {
317                          rtmp = 0;                          rtmp = 0;
318                          for (q = 0; q < p.numQuadTotal; q++)                          for (q = 0; q < p.numQuad; q++)
319                          {                          {
320                              rtmp +=                              rtmp +=
321                              vol * S[INDEX2(s, q, p.row_numShapes)] *                              vol * S[INDEX2(s, q, p.numShapes)] *
322                              D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *                              D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
323                              S[INDEX2(r, q, p.row_numShapes)];                              S[INDEX2(r, q, p.numShapes)];
324                          }                          }
325                          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)] += rtmp;
                             rtmp;  
326                          }                          }
327                      }                      }
328                      }                      }
329                  }                  }
330                  } else                  }
331                    else
332                  {                  {
333                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
334                  {                  {
335                      for (r = 0; r < p.col_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
336                      {                      {
337                      rtmp = 0;                      rtmp = 0;
338                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
339                          rtmp +=                          rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
                         vol * S[INDEX2(s, q, p.row_numShapes)] *  
                         S[INDEX2(r, q, p.row_numShapes)];  
340                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
341                      {                      {
342                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
343                          {                          {
344                          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)] +=
345                              rtmp * D_p[INDEX2(k, m, p.numEqu)];                              rtmp * D_p[INDEX2(k, m, p.numEqu)];
346                          }                          }
347                      }                      }
# Line 360  void Dudley_Assemble_PDE_System2_2D(Asse Line 358  void Dudley_Assemble_PDE_System2_2D(Asse
358                  add_EM_F = TRUE;                  add_EM_F = TRUE;
359                  if (extendedX)                  if (extendedX)
360                  {                  {
361                  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.numQuad)]);
362                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
363                  {                  {
364                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
365                      {                      {
366                      rtmp = 0;                      rtmp = 0;
367                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
368                      {                      {
369                          rtmp +=                          rtmp +=
370                          vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *                          vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
371                                 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +                                 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
372                                 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *                                 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
373                                 X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);                                 X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);
374                      }                      }
375                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp;                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
376                      }                      }
377                  }                  }
378                  } else                  }
379                    else
380                  {                  {
381                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
382                  {                  {
383                      rtmp0 = 0;                      rtmp0 = 0;
384                      rtmp1 = 0;                      rtmp1 = 0;
385                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
386                      {                      {
387                      rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];                      rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
388                      rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];                      rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
389                      }                      }
390                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
391                      EM_F[INDEX2(k, s, p.numEqu)] +=                      EM_F[INDEX2(k, s, p.numEqu)] +=
# Line 403  void Dudley_Assemble_PDE_System2_2D(Asse Line 402  void Dudley_Assemble_PDE_System2_2D(Asse
402                  add_EM_F = TRUE;                  add_EM_F = TRUE;
403                  if (extendedY)                  if (extendedY)
404                  {                  {
405                  Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);                  Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuad)]);
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                      {                      {
410                      rtmp = 0;                      rtmp = 0;
411                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
412                          rtmp +=                          rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
                         vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];  
413                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp;                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
414                      }                      }
415                  }                  }
416                  } else                  }
417                    else
418                  {                  {
419                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
420                  {                  {
421                      rtmp = 0;                      rtmp = 0;
422                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
423                      rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];                      rtmp += vol * S[INDEX2(s, q, p.numShapes)];
424                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
425                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
426                  }                  }
# Line 430  void Dudley_Assemble_PDE_System2_2D(Asse Line 429  void Dudley_Assemble_PDE_System2_2D(Asse
429                 /***********************************************************************************************/                 /***********************************************************************************************/
430              /* add the element matrices onto the matrix and right hand side                                */              /* add the element matrices onto the matrix and right hand side                                */
431                 /***********************************************************************************************/                 /***********************************************************************************************/
432              for (q = 0; q < p.row_numShapesTotal; q++)              for (q = 0; q < p.numShapes; q++)
433                  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)]];
434    
435              if (add_EM_F)              if (add_EM_F)
436                  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, p.row_DOF_UpperBound);
                            p.row_DOF_UpperBound);  
437              if (add_EM_S)              if (add_EM_S)
438                  Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,                  Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
439                                    p.col_numShapesTotal, row_index, p.numComp, EM_S);                                    p.numShapes, row_index, p.numComp, EM_S);
440    
441              }       /* end color check */              }       /* end color check */
442          }       /* end element loop */          }       /* end element loop */

Legend:
Removed from v.3187  
changed lines
  Added in v.3231

  ViewVC Help
Powered by ViewVC 1.1.26