/[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 3231 by jfenwick, Fri Oct 1 01:53:46 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      const double *S = p.shapeFns;      const double *S = p.shapeFns;
70      dim_t len_EM_S = p.numShapes * p.numShapes * p.numEqu * p.numComp;      dim_t len_EM_S = p.numShapes * p.numShapes * p.numEqu * p.numComp;
# Line 70  void Dudley_Assemble_PDE_System2_2D(Asse Line 73  void Dudley_Assemble_PDE_System2_2D(Asse
73  #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)  #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.numShapes, 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 85  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 92  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);
             double vol = p.row_jac->absD[e] * p.row_jac->quadweight;  
99              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.numShapes, DIM, p.numQuad, 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;
# Line 101  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              {              {
# Line 176  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              {              {
# Line 236  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              {              {
# Line 296  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              {              {
# Line 349  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              {              {
# Line 393  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              {              {
# Line 426  void Dudley_Assemble_PDE_System2_2D(Asse Line 429  void Dudley_Assemble_PDE_System2_2D(Asse
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.numShapes; 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    
# Line 442  void Dudley_Assemble_PDE_System2_2D(Asse Line 445  void Dudley_Assemble_PDE_System2_2D(Asse
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.3231  
changed lines
  Added in v.4657

  ViewVC Help
Powered by ViewVC 1.1.26