/[escript]/branches/domexper/dudley/src/Assemble_PDE_Single2_3D.c
ViewVC logotype

Diff of /branches/domexper/dudley/src/Assemble_PDE_Single2_3D.c

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_Single2_3D(Asse Line 65  void Dudley_Assemble_PDE_Single2_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_numShapes * p.row_numShapesTotal;      dim_t len_EM_S = p.numShapes * p.numShapes;
69      dim_t len_EM_F = p.row_numShapes;      dim_t len_EM_F = p.numShapes;
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,rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2,add_EM_F, add_EM_S)  #pragma omp parallel private(color,EM_S, EM_F, Vol, DSDX, A_p, B_p, C_p, D_p, X_p, Y_p, A_q, B_q, C_q, D_q, X_q, Y_q,row_index,q, s,r,rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2,add_EM_F, add_EM_S)
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_numShapes, 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 95  void Dudley_Assemble_PDE_Single2_3D(Asse Line 95  void Dudley_Assemble_PDE_Single2_3D(Asse
95    
96  //                      Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);  //                      Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);
97              double vol = p.row_jac->absD[e] * p.row_jac->quadweight;              double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
98              DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapes, DIM, p.numQuadTotal, 1)]);              DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.numShapes, DIM, p.numQuadTotal, 1)]);
99              for (q = 0; q < len_EM_S; ++q)              for (q = 0; q < len_EM_S; ++q)
100                  EM_S[q] = 0;                  EM_S[q] = 0;
101              for (q = 0; q < len_EM_F; ++q)              for (q = 0; q < len_EM_F; ++q)
# Line 112  void Dudley_Assemble_PDE_Single2_3D(Asse Line 112  void Dudley_Assemble_PDE_Single2_3D(Asse
112                  if (extendedA)                  if (extendedA)
113                  {                  {
114                  A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuadTotal)]);                  A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuadTotal)]);
115                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
116                  {                  {
117                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
118                      {                      {
119                      rtmp = 0;                      rtmp = 0;
120                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
121                      {                      {
122                          rtmp +=                          rtmp +=
123                          vol * (DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *                          vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
124                                 A_q[INDEX3(0, 0, q, DIM, DIM)] *                                 A_q[INDEX3(0, 0, q, DIM, DIM)] *
125                                 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)] +                                 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
126                                 DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *                                 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
127                                 A_q[INDEX3(0, 1, q, DIM, DIM)] *                                 A_q[INDEX3(0, 1, q, DIM, DIM)] *
128                                 DSDX[INDEX3(r, 1, q, p.row_numShapes, DIM)] +                                 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
129                                 DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *                                 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
130                                 A_q[INDEX3(0, 2, q, DIM, DIM)] *                                 A_q[INDEX3(0, 2, q, DIM, DIM)] *
131                                 DSDX[INDEX3(r, 2, q, p.row_numShapes, DIM)] +                                 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
132                                 DSDX[INDEX3(s, 1, q, p.row_numShapes, DIM)] *                                 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
133                                 A_q[INDEX3(1, 0, q, DIM, DIM)] *                                 A_q[INDEX3(1, 0, q, DIM, DIM)] *
134                                 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)] +                                 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
135                                 DSDX[INDEX3(s, 1, q, p.row_numShapes, DIM)] *                                 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
136                                 A_q[INDEX3(1, 1, q, DIM, DIM)] *                                 A_q[INDEX3(1, 1, q, DIM, DIM)] *
137                                 DSDX[INDEX3(r, 1, q, p.row_numShapes, DIM)] +                                 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
138                                 DSDX[INDEX3(s, 1, q, p.row_numShapes, DIM)] *                                 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
139                                 A_q[INDEX3(1, 2, q, DIM, DIM)] *                                 A_q[INDEX3(1, 2, q, DIM, DIM)] *
140                                 DSDX[INDEX3(r, 2, q, p.row_numShapes, DIM)] +                                 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
141                                 DSDX[INDEX3(s, 2, q, p.row_numShapes, DIM)] *                                 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
142                                 A_q[INDEX3(2, 0, q, DIM, DIM)] *                                 A_q[INDEX3(2, 0, q, DIM, DIM)] *
143                                 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)] +                                 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
144                                 DSDX[INDEX3(s, 2, q, p.row_numShapes, DIM)] *                                 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
145                                 A_q[INDEX3(2, 1, q, DIM, DIM)] *                                 A_q[INDEX3(2, 1, q, DIM, DIM)] *
146                                 DSDX[INDEX3(r, 1, q, p.row_numShapes, DIM)] +                                 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
147                                 DSDX[INDEX3(s, 2, q, p.row_numShapes, DIM)] *                                 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
148                                 A_q[INDEX3(2, 2, q, DIM, DIM)] *                                 A_q[INDEX3(2, 2, q, DIM, DIM)] *
149                                 DSDX[INDEX3(r, 2, q, p.row_numShapes, DIM)]);                                 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
150                      }                      }
151                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
152                      }                      }
153                  }                  }
154                  } else                  } else
155                  {                  {
156                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
157                  {                  {
158                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
159                      {                      {
160                      rtmp00 = 0;                      rtmp00 = 0;
161                      rtmp01 = 0;                      rtmp01 = 0;
# Line 169  void Dudley_Assemble_PDE_Single2_3D(Asse Line 169  void Dudley_Assemble_PDE_Single2_3D(Asse
169                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
170                      {                      {
171    
172                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)];                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
173                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
174                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapes, DIM)];                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
175                          rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.row_numShapes, DIM)];                          rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
176    
177                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapes, DIM)];                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
178                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
179                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapes, DIM)];                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
180                          rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.row_numShapes, DIM)];                          rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
181    
182                          rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.row_numShapes, DIM)];                          rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
183                          rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];                          rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
184                          rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.row_numShapes, DIM)];                          rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
185                          rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.row_numShapes, DIM)];                          rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
186                      }                      }
187                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
188                          rtmp00 * A_p[INDEX2(0, 0, DIM)] + rtmp01 * A_p[INDEX2(0, 1, DIM)] +                          rtmp00 * A_p[INDEX2(0, 0, DIM)] + rtmp01 * A_p[INDEX2(0, 1, DIM)] +
189                          rtmp02 * A_p[INDEX2(0, 2, DIM)] + rtmp10 * A_p[INDEX2(1, 0, DIM)] +                          rtmp02 * A_p[INDEX2(0, 2, DIM)] + rtmp10 * A_p[INDEX2(1, 0, DIM)] +
190                          rtmp11 * A_p[INDEX2(1, 1, DIM)] + rtmp12 * A_p[INDEX2(1, 2, DIM)] +                          rtmp11 * A_p[INDEX2(1, 1, DIM)] + rtmp12 * A_p[INDEX2(1, 2, DIM)] +
# Line 203  void Dudley_Assemble_PDE_Single2_3D(Asse Line 203  void Dudley_Assemble_PDE_Single2_3D(Asse
203                  if (extendedB)                  if (extendedB)
204                  {                  {
205                  B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);                  B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
206                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
207                  {                  {
208                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
209                      {                      {
210                      rtmp = 0;                      rtmp = 0;
211                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
212                      {                      {
213                          rtmp += vol * S[INDEX2(r, q, p.row_numShapes)] *                          rtmp += vol * S[INDEX2(r, q, p.numShapes)] *
214                          (DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *                          (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
215                           B_q[INDEX2(0, q, DIM)] +                           B_q[INDEX2(0, q, DIM)] +
216                           DSDX[INDEX3(s, 1, q, p.row_numShapes, DIM)] *                           DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
217                           B_q[INDEX2(1, q, DIM)] +                           B_q[INDEX2(1, q, DIM)] +
218                           DSDX[INDEX3(s, 2, q, p.row_numShapes, DIM)] *                           DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
219                           B_q[INDEX2(2, q, DIM)]);                           B_q[INDEX2(2, q, DIM)]);
220                      }                      }
221                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
222                      }                      }
223                  }                  }
224                  } else                  } else
225                  {                  {
226                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
227                  {                  {
228                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
229                      {                      {
230                      rtmp0 = 0;                      rtmp0 = 0;
231                      rtmp1 = 0;                      rtmp1 = 0;
232                      rtmp2 = 0;                      rtmp2 = 0;
233                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
234                      {                      {
235                          rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];                          rtmp = vol * S[INDEX2(r, q, p.numShapes)];
236                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)];                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
237                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapes, DIM)];                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
238                          rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.row_numShapes, DIM)];                          rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
239                      }                      }
240                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
241                          rtmp0 * B_p[0] + rtmp1 * B_p[1] + rtmp2 * B_p[2];                          rtmp0 * B_p[0] + rtmp1 * B_p[1] + rtmp2 * B_p[2];
242                      }                      }
243                  }                  }
# Line 252  void Dudley_Assemble_PDE_Single2_3D(Asse Line 252  void Dudley_Assemble_PDE_Single2_3D(Asse
252                  if (extendedC)                  if (extendedC)
253                  {                  {
254                  C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);                  C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
255                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
256                  {                  {
257                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
258                      {                      {
259                      rtmp = 0;                      rtmp = 0;
260                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
261                      {                      {
262                          rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] *                          rtmp += vol * S[INDEX2(s, q, p.numShapes)] *
263                          (C_q[INDEX2(0, q, DIM)] *                          (C_q[INDEX2(0, q, DIM)] *
264                           DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)] +                           DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
265                           C_q[INDEX2(1, q, DIM)] *                           C_q[INDEX2(1, q, DIM)] *
266                           DSDX[INDEX3(r, 1, q, p.row_numShapes, DIM)] +                           DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
267                           C_q[INDEX2(2, q, DIM)] *                           C_q[INDEX2(2, q, DIM)] *
268                           DSDX[INDEX3(r, 2, q, p.row_numShapes, DIM)]);                           DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
269                      }                      }
270                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
271                      }                      }
272                  }                  }
273                  } else                  } else
274                  {                  {
275                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
276                  {                  {
277                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
278                      {                      {
279                      rtmp0 = 0;                      rtmp0 = 0;
280                      rtmp1 = 0;                      rtmp1 = 0;
281                      rtmp2 = 0;                      rtmp2 = 0;
282                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
283                      {                      {
284                          rtmp = vol * S[INDEX2(s, q, p.row_numShapes)];                          rtmp = vol * S[INDEX2(s, q, p.numShapes)];
285                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
286                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapes, DIM)];                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
287                          rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.row_numShapes, DIM)];                          rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
288                      }                      }
289                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
290                          rtmp0 * C_p[0] + rtmp1 * C_p[1] + rtmp2 * C_p[2];                          rtmp0 * C_p[0] + rtmp1 * C_p[1] + rtmp2 * C_p[2];
291                      }                      }
292                  }                  }
# Line 301  void Dudley_Assemble_PDE_Single2_3D(Asse Line 301  void Dudley_Assemble_PDE_Single2_3D(Asse
301                  if (extendedD)                  if (extendedD)
302                  {                  {
303                  D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);                  D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);
304                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
305                  {                  {
306                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
307                      {                      {
308                      rtmp = 0;                      rtmp = 0;
309                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
310                          rtmp +=                          rtmp +=
311                          vol * S[INDEX2(s, q, p.row_numShapes)] * D_q[q] *                          vol * S[INDEX2(s, q, p.numShapes)] * D_q[q] *
312                          S[INDEX2(r, q, p.row_numShapes)];                          S[INDEX2(r, q, p.numShapes)];
313                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
314                      }                      }
315                  }                  }
316                  } else                  } else
317                  {                  {
318                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
319                  {                  {
320                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
321                      {                      {
322                      rtmp = 0;                      rtmp = 0;
323                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
324                          rtmp +=                          rtmp +=
325                          vol * S[INDEX2(s, q, p.row_numShapes)] *                          vol * S[INDEX2(s, q, p.numShapes)] *
326                          S[INDEX2(r, q, p.row_numShapes)];                          S[INDEX2(r, q, p.numShapes)];
327                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
328                          rtmp * D_p[0];                          rtmp * D_p[0];
329                      }                      }
330                  }                  }
# Line 339  void Dudley_Assemble_PDE_Single2_3D(Asse Line 339  void Dudley_Assemble_PDE_Single2_3D(Asse
339                  if (extendedX)                  if (extendedX)
340                  {                  {
341                  X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);                  X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
342                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
343                  {                  {
344                      rtmp = 0;                      rtmp = 0;
345                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
346                      {                      {
347                      rtmp +=                      rtmp +=
348                          vol * (DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *                          vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
349                             X_q[INDEX2(0, q, DIM)] +                             X_q[INDEX2(0, q, DIM)] +
350                             DSDX[INDEX3(s, 1, q, p.row_numShapes, DIM)] *                             DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
351                             X_q[INDEX2(1, q, DIM)] +                             X_q[INDEX2(1, q, DIM)] +
352                             DSDX[INDEX3(s, 2, q, p.row_numShapes, DIM)] *                             DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
353                             X_q[INDEX2(2, q, DIM)]);                             X_q[INDEX2(2, q, DIM)]);
354                      }                      }
355                      EM_F[INDEX2(0, s, p.numEqu)] += rtmp;                      EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
356                  }                  }
357                  } else                  } else
358                  {                  {
359                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
360                  {                  {
361                      rtmp0 = 0;                      rtmp0 = 0;
362                      rtmp1 = 0;                      rtmp1 = 0;
363                      rtmp2 = 0;                      rtmp2 = 0;
364                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
365                      {                      {
366                      rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)];                      rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
367                      rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapes, DIM)];                      rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
368                      rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.row_numShapes, DIM)];                      rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
369                      }                      }
370                      EM_F[INDEX2(0, s, p.numEqu)] += rtmp0 * X_p[0] + rtmp1 * X_p[1] + rtmp2 * X_p[2];                      EM_F[INDEX2(0, s, p.numEqu)] += rtmp0 * X_p[0] + rtmp1 * X_p[1] + rtmp2 * X_p[2];
371                  }                  }
# Line 380  void Dudley_Assemble_PDE_Single2_3D(Asse Line 380  void Dudley_Assemble_PDE_Single2_3D(Asse
380                  if (extendedY)                  if (extendedY)
381                  {                  {
382                  Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);                  Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);
383                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
384                  {                  {
385                      rtmp = 0;                      rtmp = 0;
386                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
387                      rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[q];                      rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[q];
388                      EM_F[INDEX2(0, s, p.numEqu)] += rtmp;                      EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
389                  }                  }
390                  } else                  } else
391                  {                  {
392                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
393                  {                  {
394                      rtmp = 0;                      rtmp = 0;
395                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuadTotal; q++)
396                      rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];                      rtmp += vol * S[INDEX2(s, q, p.numShapes)];
397                      EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];                      EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
398                  }                  }
399                  }                  }
# Line 402  void Dudley_Assemble_PDE_Single2_3D(Asse Line 402  void Dudley_Assemble_PDE_Single2_3D(Asse
402              /* add the element matrices onto the matrix and right hand side                                */              /* add the element matrices onto the matrix and right hand side                                */
403                 /***********************************************************************************************/                 /***********************************************************************************************/
404    
405              for (q = 0; q < p.row_numShapes; q++)              for (q = 0; q < p.numShapes; q++)
406                  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)]];
407    
408              if (add_EM_F)              if (add_EM_F)
409                  Dudley_Util_AddScatter(p.row_numShapes, row_index, p.numEqu, EM_F, F_p,                  Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p,
410                             p.row_DOF_UpperBound);                             p.row_DOF_UpperBound);
411              if (add_EM_S)              if (add_EM_S)
412                  Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapes, row_index, p.numEqu,                  Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
413                                    p.row_numShapesTotal, row_index, p.numComp, EM_S);                                    p.numShapes, row_index, p.numComp, EM_S);
414    
415              }       /* end color check */              }       /* end color check */
416          }       /* end element loop */          }       /* end element loop */

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

  ViewVC Help
Powered by ViewVC 1.1.26