/[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

branches/domexper/dudley/src/Assemble_PDE_System2_2D.c revision 3203 by jfenwick, Thu Sep 23 23:51:26 2010 UTC trunk/dudley/src/Assemble_PDE_System2_2D.c revision 4154 by jfenwick, Tue Jan 22 09:30:23 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2013 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16  /**************************************************************/  /************************************************************************************/
17    
18  /*    assembles the system of numEq PDEs into the stiffness matrix S right hand side F  */  /*    assembles the system of numEq PDEs into the stiffness matrix S right hand side F  */
19  /*    the shape functions for test and solution must be identical */  /*    the shape functions for test and solution must be identical */
# Line 30  Line 32 
32  /*      X = p.numEqu x 2  */  /*      X = p.numEqu x 2  */
33  /*      Y = p.numEqu   */  /*      Y = p.numEqu   */
34    
35  /**************************************************************/  /************************************************************************************/
36    
37  #include "Assemble.h"  #include "Assemble.h"
38  #include "Util.h"  #include "Util.h"
# Line 38  Line 40 
40  #include <omp.h>  #include <omp.h>
41  #endif  #endif
42    
43  /**************************************************************/  /************************************************************************************/
44    
45  void Dudley_Assemble_PDE_System2_2D(Assemble_Parameters p, Dudley_ElementFile * elements,  void Dudley_Assemble_PDE_System2_2D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements,
46                      Paso_SystemMatrix * Mat, escriptDataC * F,                      Paso_SystemMatrix * Mat, escriptDataC * F,
47                      escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,                      escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,
48                      escriptDataC * X, escriptDataC * Y)                      escriptDataC * X, escriptDataC * Y)
# Line 63  void Dudley_Assemble_PDE_System2_2D(Asse Line 65  void Dudley_Assemble_PDE_System2_2D(Asse
65      bool_t extendedX = isExpanded(X);      bool_t extendedX = isExpanded(X);
66      bool_t extendedY = isExpanded(Y);      bool_t extendedY = isExpanded(Y);
67      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 */
68  //    double *S = p.row_jac->BasisFunctions->S;      const double *S = p.shapeFns;
69      const double* S = p.shapeFns;      dim_t len_EM_S = p.numShapes * p.numShapes * p.numEqu * p.numComp;
70      dim_t len_EM_S = p.row_numShapesTotal * p.row_numShapesTotal * p.numEqu * p.numComp;      dim_t len_EM_F = p.numShapes * p.numEqu;
     dim_t len_EM_F = p.row_numShapesTotal * p.numEqu;  
71    
72  #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)
73      {      {
74    
75      EM_S = THREAD_MEMALLOC(len_EM_S, double);      EM_S = THREAD_MEMALLOC(len_EM_S, double);
76      EM_F = THREAD_MEMALLOC(len_EM_F, double);      EM_F = THREAD_MEMALLOC(len_EM_F, double);
77      row_index = THREAD_MEMALLOC(p.row_numShapesTotal, index_t);      row_index = THREAD_MEMALLOC(p.numShapes, index_t);
78    
79      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))
80      {      {
# Line 86  void Dudley_Assemble_PDE_System2_2D(Asse Line 87  void Dudley_Assemble_PDE_System2_2D(Asse
87          {          {
88              if (elements->Color[e] == color)              if (elements->Color[e] == color)
89              {              {
90                double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
91    
92              A_p = getSampleDataRO(A, e);              A_p = getSampleDataRO(A, e);
93              B_p = getSampleDataRO(B, e);              B_p = getSampleDataRO(B, e);
# Line 93  void Dudley_Assemble_PDE_System2_2D(Asse Line 95  void Dudley_Assemble_PDE_System2_2D(Asse
95              D_p = getSampleDataRO(D, e);              D_p = getSampleDataRO(D, e);
96              X_p = getSampleDataRO(X, e);              X_p = getSampleDataRO(X, e);
97              Y_p = getSampleDataRO(Y, e);              Y_p = getSampleDataRO(Y, e);
98              double vol = p.row_jac->absD[e] * p.row_jac->quadweight;              DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.numShapes, DIM, p.numQuad, 1)]);
             DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapesTotal, 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 102  void Dudley_Assemble_PDE_System2_2D(Asse Line 103  void Dudley_Assemble_PDE_System2_2D(Asse
103              add_EM_F = FALSE;              add_EM_F = FALSE;
104              add_EM_S = FALSE;              add_EM_S = FALSE;
105    
106                /**************************************************************/                /************************************************************************************/
107              /*   process A: */              /*   process A: */
108                /**************************************************************/                /************************************************************************************/
109              A_p = getSampleDataRO(A, e);              A_p = getSampleDataRO(A, e);
110              if (NULL != A_p)              if (NULL != A_p)
111              {              {
112                  add_EM_S = TRUE;                  add_EM_S = TRUE;
113                  if (extendedA)                  if (extendedA)
114                  {                  {
115                  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)]);
116                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
117                  {                  {
118                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
119                      {                      {
120                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
121                      {                      {
122                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
123                          {                          {
124                          rtmp = 0;                          rtmp = 0;
125                          for (q = 0; q < p.numQuadTotal; q++)                          for (q = 0; q < p.numQuad; q++)
126                          {                          {
127                              rtmp +=                              rtmp +=
128                              vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *                              vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
129                                     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)]
130                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +                                     * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
131                                     DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *                                     DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
132                                     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)]
133                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +                                     * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
134                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *                                     DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
135                                     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)]
136                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +                                     * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
137                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *                                     DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
138                                     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)]
139                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);                                     * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
140                          }                          }
141                          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;  
142                          }                          }
143                      }                      }
144                      }                      }
145                  }                  }
146                  } else                  }
147                    else
148                  {                  {
149                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
150                  {                  {
151                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
152                      {                      {
153                      rtmp00 = 0;                      rtmp00 = 0;
154                      rtmp01 = 0;                      rtmp01 = 0;
155                      rtmp10 = 0;                      rtmp10 = 0;
156                      rtmp11 = 0;                      rtmp11 = 0;
157                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
158                      {                      {
159                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
160                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
161                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
162                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
163                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
164                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
165                      }                      }
166                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
167                      {                      {
168                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
169                          {                          {
170                          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)] +=
171                              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)]
172                              + 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)]
173                              + 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 177  void Dudley_Assemble_PDE_System2_2D(Asse Line 178  void Dudley_Assemble_PDE_System2_2D(Asse
178                  }                  }
179                  }                  }
180              }              }
181                /**************************************************************/                /************************************************************************************/
182              /*   process B: */              /*   process B: */
183                /**************************************************************/                /************************************************************************************/
184              B_p = getSampleDataRO(B, e);              B_p = getSampleDataRO(B, e);
185              if (NULL != B_p)              if (NULL != B_p)
186              {              {
187                  add_EM_S = TRUE;                  add_EM_S = TRUE;
188                  if (extendedB)                  if (extendedB)
189                  {                  {
190                  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)]);
191                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
192                  {                  {
193                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
194                      {                      {
195                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
196                      {                      {
197                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
198                          {                          {
199                          rtmp = 0;                          rtmp = 0;
200                          for (q = 0; q < p.numQuadTotal; q++)                          for (q = 0; q < p.numQuad; q++)
201                          {                          {
202                              rtmp += vol * S[INDEX2(r, q, p.row_numShapes)] *                              rtmp += vol * S[INDEX2(r, q, p.numShapes)] *
203                              (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *                              (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
204                               B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +                               B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
205                               DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *                               DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
206                               B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);                               B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);
207                          }                          }
208                          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;  
209                          }                          }
210                      }                      }
211                      }                      }
212                  }                  }
213                  } else                  }
214                    else
215                  {                  {
216                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
217                  {                  {
218                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
219                      {                      {
220                      rtmp0 = 0;                      rtmp0 = 0;
221                      rtmp1 = 0;                      rtmp1 = 0;
222                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
223                      {                      {
224                          rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];                          rtmp = vol * S[INDEX2(r, q, p.numShapes)];
225                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
226                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
227                      }                      }
228                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
229                      {                      {
230                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
231                          {                          {
232                          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)] +=
233                              rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +                              rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
234                              rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];                              rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];
235                          }                          }
# Line 237  void Dudley_Assemble_PDE_System2_2D(Asse Line 238  void Dudley_Assemble_PDE_System2_2D(Asse
238                  }                  }
239                  }                  }
240              }              }
241                /**************************************************************/                /************************************************************************************/
242              /*   process C: */              /*   process C: */
243                /**************************************************************/                /************************************************************************************/
244              C_p = getSampleDataRO(C, e);              C_p = getSampleDataRO(C, e);
245              if (NULL != C_p)              if (NULL != C_p)
246              {              {
247                  add_EM_S = TRUE;                  add_EM_S = TRUE;
248                  if (extendedC)                  if (extendedC)
249                  {                  {
250                  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)]);
251                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
252                  {                  {
253                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
254                      {                      {
255                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
256                      {                      {
257                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
258                          {                          {
259                          rtmp = 0;                          rtmp = 0;
260                          for (q = 0; q < p.numQuadTotal; q++)                          for (q = 0; q < p.numQuad; q++)
261                          {                          {
262                              rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] *                              rtmp += vol * S[INDEX2(s, q, p.numShapes)] *
263                              (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *                              (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
264                               DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +                               DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
265                               C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *                               C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
266                               DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);                               DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
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)] += rtmp;
                             rtmp;  
269                          }                          }
270                      }                      }
271                      }                      }
272                  }                  }
273                  } else                  }
274                    else
275                  {                  {
276                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
277                  {                  {
278                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
279                      {                      {
280                      rtmp0 = 0;                      rtmp0 = 0;
281                      rtmp1 = 0;                      rtmp1 = 0;
282                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; 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_numShapesTotal, DIM)];                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
286                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
287                      }                      }
288                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
289                      {                      {
290                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
291                          {                          {
292                          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)] +=
293                              rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +                              rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
294                              rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];                              rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];
295                          }                          }
# Line 297  void Dudley_Assemble_PDE_System2_2D(Asse Line 298  void Dudley_Assemble_PDE_System2_2D(Asse
298                  }                  }
299                  }                  }
300              }              }
301                /************************************************************* */                /*********************************************************************************** */
302              /* process D */              /* process D */
303                /**************************************************************/                /************************************************************************************/
304              D_p = getSampleDataRO(D, e);              D_p = getSampleDataRO(D, e);
305              if (NULL != D_p)              if (NULL != D_p)
306              {              {
307                  add_EM_S = TRUE;                  add_EM_S = TRUE;
308                  if (extendedD)                  if (extendedD)
309                  {                  {
310                  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)]);
311                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
312                  {                  {
313                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
314                      {                      {
315                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
316                      {                      {
317                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
318                          {                          {
319                          rtmp = 0;                          rtmp = 0;
320                          for (q = 0; q < p.numQuadTotal; q++)                          for (q = 0; q < p.numQuad; q++)
321                          {                          {
322                              rtmp +=                              rtmp +=
323                              vol * S[INDEX2(s, q, p.row_numShapes)] *                              vol * S[INDEX2(s, q, p.numShapes)] *
324                              D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *                              D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
325                              S[INDEX2(r, q, p.row_numShapes)];                              S[INDEX2(r, q, p.numShapes)];
326                          }                          }
327                          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;  
328                          }                          }
329                      }                      }
330                      }                      }
331                  }                  }
332                  } else                  }
333                    else
334                  {                  {
335                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
336                  {                  {
337                      for (r = 0; r < p.row_numShapes; r++)                      for (r = 0; r < p.numShapes; r++)
338                      {                      {
339                      rtmp = 0;                      rtmp = 0;
340                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
341                          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)];  
342                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
343                      {                      {
344                          for (m = 0; m < p.numComp; m++)                          for (m = 0; m < p.numComp; m++)
345                          {                          {
346                          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)] +=
347                              rtmp * D_p[INDEX2(k, m, p.numEqu)];                              rtmp * D_p[INDEX2(k, m, p.numEqu)];
348                          }                          }
349                      }                      }
# Line 352  void Dudley_Assemble_PDE_System2_2D(Asse Line 351  void Dudley_Assemble_PDE_System2_2D(Asse
351                  }                  }
352                  }                  }
353              }              }
354                /**************************************************************/                /************************************************************************************/
355              /*   process X: */              /*   process X: */
356                /**************************************************************/                /************************************************************************************/
357              X_p = getSampleDataRO(X, e);              X_p = getSampleDataRO(X, e);
358              if (NULL != X_p)              if (NULL != X_p)
359              {              {
360                  add_EM_F = TRUE;                  add_EM_F = TRUE;
361                  if (extendedX)                  if (extendedX)
362                  {                  {
363                  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)]);
364                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
365                  {                  {
366                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
367                      {                      {
368                      rtmp = 0;                      rtmp = 0;
369                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
370                      {                      {
371                          rtmp +=                          rtmp +=
372                          vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *                          vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
373                                 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +                                 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
374                                 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *                                 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
375                                 X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);                                 X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);
376                      }                      }
377                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp;                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
378                      }                      }
379                  }                  }
380                  } else                  }
381                    else
382                  {                  {
383                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
384                  {                  {
385                      rtmp0 = 0;                      rtmp0 = 0;
386                      rtmp1 = 0;                      rtmp1 = 0;
387                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
388                      {                      {
389                      rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];                      rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
390                      rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];                      rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
391                      }                      }
392                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
393                      EM_F[INDEX2(k, s, p.numEqu)] +=                      EM_F[INDEX2(k, s, p.numEqu)] +=
# Line 395  void Dudley_Assemble_PDE_System2_2D(Asse Line 395  void Dudley_Assemble_PDE_System2_2D(Asse
395                  }                  }
396                  }                  }
397              }              }
398               /**************************************************************/               /************************************************************************************/
399              /*   process Y: */              /*   process Y: */
400               /**************************************************************/               /************************************************************************************/
401              Y_p = getSampleDataRO(Y, e);              Y_p = getSampleDataRO(Y, e);
402              if (NULL != Y_p)              if (NULL != Y_p)
403              {              {
404                  add_EM_F = TRUE;                  add_EM_F = TRUE;
405                  if (extendedY)                  if (extendedY)
406                  {                  {
407                  Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);                  Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuad)]);
408                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
409                  {                  {
410                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
411                      {                      {
412                      rtmp = 0;                      rtmp = 0;
413                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
414                          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)];  
415                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp;                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
416                      }                      }
417                  }                  }
418                  } else                  }
419                    else
420                  {                  {
421                  for (s = 0; s < p.row_numShapes; s++)                  for (s = 0; s < p.numShapes; s++)
422                  {                  {
423                      rtmp = 0;                      rtmp = 0;
424                      for (q = 0; q < p.numQuadTotal; q++)                      for (q = 0; q < p.numQuad; q++)
425                      rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];                      rtmp += vol * S[INDEX2(s, q, p.numShapes)];
426                      for (k = 0; k < p.numEqu; k++)                      for (k = 0; k < p.numEqu; k++)
427                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];                      EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
428                  }                  }
429                  }                  }
430              }              }
431                 /***********************************************************************************************/                 /*********************************************************************************************************************/
432              /* add the element matrices onto the matrix and right hand side                                */              /* add the element matrices onto the matrix and right hand side                                */
433                 /***********************************************************************************************/                 /*********************************************************************************************************************/
434              for (q = 0; q < p.row_numShapesTotal; q++)              for (q = 0; q < p.numShapes; q++)
435                  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)]];
436    
437              if (add_EM_F)              if (add_EM_F)
438                  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);  
439              if (add_EM_S)              if (add_EM_S)
440                  Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,                  Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
441                                    p.row_numShapesTotal, row_index, p.numComp, EM_S);                                    p.numShapes, row_index, p.numComp, EM_S);
442    
443              }       /* end color check */              }       /* end color check */
444          }       /* end element loop */          }       /* end element loop */

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

  ViewVC Help
Powered by ViewVC 1.1.26