# Diff of /branches/trilinos_from_5897/dudley/src/Assemble_PDE_System2_2D.cpp

branches/domexper/dudley/src/Assemble_PDE_System2_2D.c revision 3231 by jfenwick, Fri Oct 1 01:53:46 2010 UTC branches/trilinos_from_5897/dudley/src/Assemble_PDE_System2_2D.cpp revision 5963 by caltinay, Mon Feb 22 06:59:27 2016 UTC
# Line 1  Line 1
1
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2016 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
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,                      escript::ASM_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 52  void Dudley_Assemble_PDE_System2_2D(Asse Line 58  void Dudley_Assemble_PDE_System2_2D(Asse
58      __const double *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p, *A_q, *B_q, *C_q, *D_q, *X_q, *Y_q;      __const double *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p, *A_q, *B_q, *C_q, *D_q, *X_q, *Y_q;
59      double *EM_S, *EM_F, *DSDX;      double *EM_S, *EM_F, *DSDX;
60      index_t *row_index;      index_t *row_index;
61      register dim_t q, s, r, k, m;      dim_t q, s, r, k, m;
62      register double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11;      double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11;
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      const double *S = p.shapeFns;      const double *S = p.shapeFns;
73      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 76  void Dudley_Assemble_PDE_System2_2D(Asse
76  #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)
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.numShapes, 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 85  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 92  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);
double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
102              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)]);
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;
# Line 101  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              {              {
# Line 176  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              {              {
# Line 236  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              {              {
# Line 296  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              {              {
# Line 349  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              {              {
# Line 393  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              {              {
# Line 426  void Dudley_Assemble_PDE_System2_2D(Asse Line 432  void Dudley_Assemble_PDE_System2_2D(Asse
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.numShapes; 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
442                  Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound);                  Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound);
444                  Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,                  Dudley_Assemble_addToSystemMatrix(mat, p.numShapes, row_index, p.numEqu,
445                                    p.numShapes, 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 */
457  }  }
458
/*
* \$Log\$
*/

Legend:
 Removed from v.3231 changed lines Added in v.5963