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

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

  ViewVC Help
Powered by ViewVC 1.1.26