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

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

  ViewVC Help
Powered by ViewVC 1.1.26