/[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 3184 by jfenwick, Wed Sep 15 00:23:42 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 */
21    
   
22  /*      -(A_{k,i,m,j} u_m,j)_i-(B_{k,i,m} u_m)_i+C_{k,m,j} u_m,j-D_{k,m} u_m  and -(X_{k,i})_i + Y_k */  /*      -(A_{k,i,m,j} u_m,j)_i-(B_{k,i,m} u_m)_i+C_{k,m,j} u_m,j-D_{k,m} u_m  and -(X_{k,i})_i + Y_k */
23    
24  /*    u has p.numComp components in a 2D domain. The shape functions for test and solution must be identical  */  /*    u has p.numComp components in a 2D domain. The shape functions for test and solution must be identical  */
# Line 32  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 42  Line 44 
44  #include <omp.h>  #include <omp.h>
45  #endif  #endif
46    
47    /************************************************************************************/
48    
49  /**************************************************************/  void Dudley_Assemble_PDE_System2_2D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements,
50                        paso::SystemMatrix_ptr Mat, escript::Data* F,
51  void  Dudley_Assemble_PDE_System2_2D(Assemble_Parameters p, Dudley_ElementFile* elements,                      const escript::Data* A, const escript::Data* B, const escript::Data* C, const escript::Data* D,
52                                       Paso_SystemMatrix* Mat, escriptDataC* F,                      const escript::Data* X, const escript::Data* Y)
53                                       escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y) {  {
54    
55      #define DIM 2  #define DIM 2
56      index_t color;      index_t color;
57      dim_t e;      dim_t e;
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;      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      dim_t len_EM_S=p.row_numShapesTotal*p.col_numShapesTotal*p.numEqu*p.numComp;      dim_t len_EM_S = p.numShapes * p.numShapes * p.numEqu * p.numComp;
74      dim_t len_EM_F=p.row_numShapesTotal*p.numEqu;      dim_t len_EM_F = p.numShapes * 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        {
85            for (color=elements->minColor;color<=elements->maxColor;color++) {  
86               /*  open loop over all elements: */          for (color = elements->minColor; color <= elements->maxColor; color++)
87               #pragma omp for private(e) schedule(static)          {
88               for(e=0;e<elements->numElements;e++){          /*  open loop over all elements: */
89                  if (elements->Color[e]==color) {  #pragma omp for private(e) schedule(static)
90                                  for (e = 0; e < elements->numElements; e++)
91             A_p=getSampleDataRO(A,e);          {
92             B_p=getSampleDataRO(B,e);              if (elements->Color[e] == color)
93                     C_p=getSampleDataRO(C,e);                  {
94                     D_p=getSampleDataRO(D,e);              double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
95                     X_p=getSampleDataRO(X,e);  
96                     Y_p=getSampleDataRO(Y,e);              A_p = getSampleDataRO(A, e);
97            double vol=p.row_jac->absD[e]*p.row_jac->quadweight;                      B_p = getSampleDataRO(B, e);
98                        DSDX=&(p.row_jac->DSDX[INDEX5(0,0,0,0,e,p.row_numShapesTotal,DIM,p.numQuadTotal,1)]);              C_p = getSampleDataRO(C, e);
99                        for (q=0;q<len_EM_S;++q) EM_S[q]=0;              D_p = getSampleDataRO(D, e);
100                        for (q=0;q<len_EM_F;++q) EM_F[q]=0;              X_p = getSampleDataRO(X, e);
101                        add_EM_F=FALSE;              Y_p = getSampleDataRO(Y, e);
102                        add_EM_S=FALSE;              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)
104                        /**************************************************************/                  EM_S[q] = 0;
105                        /*   process A: */              for (q = 0; q < len_EM_F; ++q)
106                        /**************************************************************/                  EM_F[q] = 0;
107                        A_p=getSampleDataRO(A,e);              add_EM_F = FALSE;
108                        if (NULL!=A_p) {              add_EM_S = FALSE;
109                           add_EM_S=TRUE;  
110                           if (extendedA) {                /************************************************************************************/
111                  A_q=&(A_p[INDEX6(0,0,0,0,0,0, p.numEqu,DIM,p.numComp,DIM,p.numQuadTotal)]);              /*   process A: */
112                              for (s=0;s<p.row_numShapes;s++) {                /************************************************************************************/
113                                for (r=0;r<p.col_numShapes;r++) {              A_p = getSampleDataRO(A, e);
114                                  for (k=0;k<p.numEqu;k++) {              if (NULL != A_p)
115                                    for (m=0;m<p.numComp;m++) {              {
116                                      rtmp=0;                  add_EM_S = TRUE;
117                                      for (q=0;q<p.numQuadTotal;q++) {                  if (extendedA)
118                                         rtmp+=vol* (                  {
119                                             DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]                  A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuad)]);
120                                            +DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]                  for (s = 0; s < p.numShapes; s++)
121                                            +DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,1,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]                  {
122                                            +DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,1,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]);                      for (r = 0; r < p.numShapes; r++)
123                                      }                      {
124                                      EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                      for (k = 0; k < p.numEqu; k++)
125                                    }                      {
126                                  }                          for (m = 0; m < p.numComp; m++)
127                                }                          {
128                              }                          rtmp = 0;
129                           } else {                          for (q = 0; q < p.numQuad; q++)
130                              for (s=0;s<p.row_numShapes;s++) {                          {
131                                for (r=0;r<p.col_numShapes;r++) {                              rtmp +=
132                                   rtmp00=0;                              vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
133                                   rtmp01=0;                                     A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
134                                   rtmp10=0;                                     * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
135                                   rtmp11=0;                                     DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
136                                   for (q=0;q<p.numQuadTotal;q++) {                                     A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
137                                         rtmp0=vol*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                                     * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
138                                         rtmp1=vol*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];                                     DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
139                                         rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                                     A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
140                                         rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                                     * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
141                                         rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                                     DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
142                                         rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                                     A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
143                                   }                                     * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
144                                   for (k=0;k<p.numEqu;k++) {                          }
145                                       for (m=0;m<p.numComp;m++) {                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
146                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=                          }
147                                                     rtmp00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]                      }
148                                                    +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]                      }
149                                                    +rtmp10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]                  }
150                                                    +rtmp11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)];                  }
151                                       }                  else
152                                   }                  {
153                                }                  for (s = 0; s < p.numShapes; s++)
154                              }                  {
155                          }                      for (r = 0; r < p.numShapes; r++)
156                        }                      {
157                        /**************************************************************/                      rtmp00 = 0;
158                        /*   process B: */                      rtmp01 = 0;
159                        /**************************************************************/                      rtmp10 = 0;
160                        B_p=getSampleDataRO(B,e);                      rtmp11 = 0;
161                        if (NULL!=B_p) {                      for (q = 0; q < p.numQuad; q++)
162                           add_EM_S=TRUE;                      {
163                           if (extendedB) {                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
164                      B_q=&(B_p[INDEX5(0,0,0,0,0 ,p.numEqu,DIM,p.numComp,p.numQuadTotal)]);                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
165                              for (s=0;s<p.row_numShapes;s++) {                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
166                                for (r=0;r<p.col_numShapes;r++) {                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
167                                  for (k=0;k<p.numEqu;k++) {                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
168                                    for (m=0;m<p.numComp;m++) {                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
169                                       rtmp=0;                      }
170                                       for (q=0;q<p.numQuadTotal;q++) {                      for (k = 0; k < p.numEqu; k++)
171                                          rtmp+=vol*S[INDEX2(r,q,p.row_numShapes)]*                      {
172                                               ( DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]                          for (m = 0; m < p.numComp; m++)
173                                               + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]);                          {
174                                       }                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
175                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                              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)]
177                                  }                              + rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)]
178                                }                              + rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)];
179                              }                          }
180                           } else {                      }
181                              for (s=0;s<p.row_numShapes;s++) {                      }
182                                for (r=0;r<p.col_numShapes;r++) {                  }
183                                    rtmp0=0;                  }
184                                    rtmp1=0;              }
185                                    for (q=0;q<p.numQuadTotal;q++) {                /************************************************************************************/
186                                        rtmp=vol*S[INDEX2(r,q,p.row_numShapes)];              /*   process B: */
187                                        rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                /************************************************************************************/
188                                        rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];              B_p = getSampleDataRO(B, e);
189                                    }              if (NULL != B_p)
190                                    for (k=0;k<p.numEqu;k++) {              {
191                                       for (m=0;m<p.numComp;m++) {                  add_EM_S = TRUE;
192                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+= rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]                  if (extendedB)
193                                                                                            + rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)];                  {
194                                       }                  B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuad)]);
195                                    }                  for (s = 0; s < p.numShapes; s++)
196                                }                  {
197                              }                      for (r = 0; r < p.numShapes; r++)
198                           }                      {
199                        }                      for (k = 0; k < p.numEqu; k++)
200                        /**************************************************************/                      {
201                        /*   process C: */                          for (m = 0; m < p.numComp; m++)
202                        /**************************************************************/                          {
203                        C_p=getSampleDataRO(C,e);                          rtmp = 0;
204                        if (NULL!=C_p) {                          for (q = 0; q < p.numQuad; q++)
205                          add_EM_S=TRUE;                          {
206                          if (extendedC) {                              rtmp += vol * S[INDEX2(r, q, p.numShapes)] *
207                  C_q=&(C_p[INDEX5(0,0,0,0,0, p.numEqu,p.numComp,DIM,p.numQuadTotal)]);                              (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
208                              for (s=0;s<p.row_numShapes;s++) {                               B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
209                                for (r=0;r<p.col_numShapes;r++) {                               DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
210                                  for (k=0;k<p.numEqu;k++) {                               B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);
211                                    for (m=0;m<p.numComp;m++) {                          }
212                                       rtmp=0;                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
213                                       for (q=0;q<p.numQuadTotal;q++) {                          }
214                                           rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*                      }
215                                                ( C_q[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]                      }
216                                                + C_q[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]);                  }
217                                       }                  }
218                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                  else
219                                    }                  {
220                                  }                  for (s = 0; s < p.numShapes; s++)
221                                }                  {
222                              }                      for (r = 0; r < p.numShapes; r++)
223                          } else {                      {
224                              for (s=0;s<p.row_numShapes;s++) {                      rtmp0 = 0;
225                                for (r=0;r<p.col_numShapes;r++) {                      rtmp1 = 0;
226                                   rtmp0=0;                      for (q = 0; q < p.numQuad; q++)
227                                   rtmp1=0;                      {
228                                   for (q=0;q<p.numQuadTotal;q++) {                          rtmp = vol * S[INDEX2(r, q, p.numShapes)];
229                                           rtmp=vol*S[INDEX2(s,q,p.row_numShapes)];                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
230                                           rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
231                                           rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];                      }
232                                   }                      for (k = 0; k < p.numEqu; k++)
233                                   for (k=0;k<p.numEqu;k++) {                      {
234                                      for (m=0;m<p.numComp;m++) {                          for (m = 0; m < p.numComp; m++)
235                                             EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp0*C_p[INDEX3(k,m,0,p.numEqu,p.numComp)]                          {
236                                                                                               +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)];                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
237                                      }                              rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
238                                   }                              rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];
239                                }                          }
240                              }                      }
241                          }                      }
242                        }                  }
243                        /************************************************************* */                  }
244                        /* process D */              }
245                        /**************************************************************/                /************************************************************************************/
246                        D_p=getSampleDataRO(D,e);              /*   process C: */
247                        if (NULL!=D_p) {                /************************************************************************************/
248                          add_EM_S=TRUE;              C_p = getSampleDataRO(C, e);
249                          if (extendedD) {              if (NULL != C_p)
250                  D_q=&(D_p[INDEX4(0,0,0,0, p.numEqu,p.numComp,p.numQuadTotal)]);              {
251                              for (s=0;s<p.row_numShapes;s++) {                  add_EM_S = TRUE;
252                                for (r=0;r<p.col_numShapes;r++) {                  if (extendedC)
253                                  for (k=0;k<p.numEqu;k++) {                  {
254                                    for (m=0;m<p.numComp;m++) {                  C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuad)]);
255                                      rtmp=0;                  for (s = 0; s < p.numShapes; s++)
256                                      for (q=0;q<p.numQuadTotal;q++) {                  {
257                                          rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*D_q[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_numShapes)];                      for (r = 0; r < p.numShapes; r++)
258                                      }                      {
259                                      EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp;                      for (k = 0; k < p.numEqu; k++)
260                                    }                      {
261                                  }                          for (m = 0; m < p.numComp; m++)
262                                }                          {
263                              }                          rtmp = 0;
264                          } else {                          for (q = 0; q < p.numQuad; q++)
265                              for (s=0;s<p.row_numShapes;s++) {                          {
266                                for (r=0;r<p.col_numShapes;r++) {                              rtmp += vol * S[INDEX2(s, q, p.numShapes)] *
267                                    rtmp=0;                              (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
268                                    for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(r,q,p.row_numShapes)];                               DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
269                                    for (k=0;k<p.numEqu;k++) {                               C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
270                                        for (m=0;m<p.numComp;m++) {                               DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
271                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];                          }
272                                       }                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
273                                    }                          }
274                                }                      }
275                              }                      }
276                          }                  }
277                        }                  }
278                        /**************************************************************/                  else
279                        /*   process X: */                  {
280                        /**************************************************************/                  for (s = 0; s < p.numShapes; s++)
281                        X_p=getSampleDataRO(X,e);                  {
282                        if (NULL!=X_p) {                      for (r = 0; r < p.numShapes; r++)
283                          add_EM_F=TRUE;                      {
284                          if (extendedX) {                      rtmp0 = 0;
285                 X_q=&(X_p[INDEX4(0,0,0,0, p.numEqu,DIM,p.numQuadTotal)]);                      rtmp1 = 0;
286                             for (s=0;s<p.row_numShapes;s++) {                      for (q = 0; q < p.numQuad; q++)
287                                for (k=0;k<p.numEqu;k++) {                      {
288                                  rtmp=0;                          rtmp = vol * S[INDEX2(s, q, p.numShapes)];
289                                  for (q=0;q<p.numQuadTotal;q++) {                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
290                                       rtmp+=vol*(DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,0,q,p.numEqu,DIM)]                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
291                                                    +DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,1,q,p.numEqu,DIM)]);                      }
292                                  }                      for (k = 0; k < p.numEqu; k++)
293                                  EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;                      {
294                                }                          for (m = 0; m < p.numComp; m++)
295                             }                          {
296                          } else {                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
297                             for (s=0;s<p.row_numShapes;s++) {                              rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
298                                   rtmp0=0;                              rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];
299                                   rtmp1=0;                          }
300                                   for (q=0;q<p.numQuadTotal;q++) {                      }
301                                       rtmp0+=vol*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];                      }
302                                       rtmp1+=vol*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];                  }
303                                   }                  }
304                                   for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp0*X_p[INDEX2(k,0,p.numEqu)]+rtmp1*X_p[INDEX2(k,1,p.numEqu)];              }
305                             }                /*********************************************************************************** */
306                          }              /* process D */
307                       }                /************************************************************************************/
308                       /**************************************************************/              D_p = getSampleDataRO(D, e);
309                       /*   process Y: */              if (NULL != D_p)
310                       /**************************************************************/              {
311                        Y_p=getSampleDataRO(Y,e);                  add_EM_S = TRUE;
312                        if (NULL!=Y_p) {                  if (extendedD)
313                          add_EM_F=TRUE;                  {
314                          if (extendedY) {                  D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuad)]);
315                 Y_q=&(Y_p[INDEX3(0,0,0, p.numEqu,p.numQuadTotal)]);                  for (s = 0; s < p.numShapes; s++)
316                             for (s=0;s<p.row_numShapes;s++) {                  {
317                                for (k=0;k<p.numEqu;k++) {                      for (r = 0; r < p.numShapes; r++)
318                                   rtmp=0;                      {
319                                   for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*Y_q[INDEX2(k,q,p.numEqu)];                      for (k = 0; k < p.numEqu; k++)
320                                   EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;                      {
321                                }                          for (m = 0; m < p.numComp; m++)
322                             }                          {
323                           } else {                          rtmp = 0;
324                             for (s=0;s<p.row_numShapes;s++) {                          for (q = 0; q < p.numQuad; q++)
325                                 rtmp=0;                          {
326                                 for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)];                              rtmp +=
327                                 for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];                              vol * S[INDEX2(s, q, p.numShapes)] *
328                             }                              D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
329                           }                              S[INDEX2(r, q, p.numShapes)];
330                         }                          }
331                         /***********************************************************************************************/                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
332                         /* add the element matrices onto the matrix and right hand side                                */                          }
333                         /***********************************************************************************************/                      }
334                         for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];                      }
335                    }
336                         if (add_EM_F) Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);                  }
337                         if (add_EM_S) Dudley_Assemble_addToSystemMatrix(Mat,p.row_numShapesTotal,row_index,p.numEqu,p.col_numShapesTotal,row_index,p.numComp,EM_S);                  else
338                    {
339                  } /* end color check */                  for (s = 0; s < p.numShapes; s++)
340               } /* end element loop */                  {
341           } /* end color loop */                      for (r = 0; r < p.numShapes; r++)
342                                  {
343           THREAD_MEMFREE(EM_S);                      rtmp = 0;
344           THREAD_MEMFREE(EM_F);                      for (q = 0; q < p.numQuad; q++)
345           THREAD_MEMFREE(row_index);                          rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
346                        for (k = 0; k < p.numEqu; k++)
347                        {
348                            for (m = 0; m < p.numComp; m++)
349                            {
350                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
351                                rtmp * D_p[INDEX2(k, m, p.numEqu)];
352                            }
353                        }
354                        }
355                    }
356                    }
357                }
358                  /************************************************************************************/
359                /*   process X: */
360                  /************************************************************************************/
361                X_p = getSampleDataRO(X, e);
362                if (NULL != X_p)
363                {
364                    add_EM_F = TRUE;
365                    if (extendedX)
366                    {
367                    X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuad)]);
368                    for (s = 0; s < p.numShapes; s++)
369                    {
370                        for (k = 0; k < p.numEqu; k++)
371                        {
372                        rtmp = 0;
373                        for (q = 0; q < p.numQuad; q++)
374                        {
375                            rtmp +=
376                            vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
377                                   X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
378                                   DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
379                                   X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);
380                        }
381                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
382                        }
383                    }
384                    }
385                    else
386                    {
387                    for (s = 0; s < p.numShapes; s++)
388                    {
389                        rtmp0 = 0;
390                        rtmp1 = 0;
391                        for (q = 0; q < p.numQuad; q++)
392                        {
393                        rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
394                        rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
395                        }
396                        for (k = 0; k < p.numEqu; k++)
397                        EM_F[INDEX2(k, s, p.numEqu)] +=
398                            rtmp0 * X_p[INDEX2(k, 0, p.numEqu)] + rtmp1 * X_p[INDEX2(k, 1, p.numEqu)];
399                    }
400                    }
401                }
402                 /************************************************************************************/
403                /*   process Y: */
404                 /************************************************************************************/
405                Y_p = getSampleDataRO(Y, e);
406                if (NULL != Y_p)
407                {
408                    add_EM_F = TRUE;
409                    if (extendedY)
410                    {
411                    Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuad)]);
412                    for (s = 0; s < p.numShapes; s++)
413                    {
414                        for (k = 0; k < p.numEqu; k++)
415                        {
416                        rtmp = 0;
417                        for (q = 0; q < p.numQuad; q++)
418                            rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
419                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
420                        }
421                    }
422                    }
423                    else
424                    {
425                    for (s = 0; s < p.numShapes; s++)
426                    {
427                        rtmp = 0;
428                        for (q = 0; q < p.numQuad; q++)
429                        rtmp += vol * S[INDEX2(s, q, p.numShapes)];
430                        for (k = 0; k < p.numEqu; k++)
431                        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                                */
437                   /*********************************************************************************************************************/
438                for (q = 0; q < p.numShapes; q++)
439                    row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
440    
441                if (add_EM_F)
442                    Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound);
443                if (add_EM_S)
444                    Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
445                                      p.numShapes, row_index, p.numComp, EM_S);
446    
447                }       /* end color check */
448            }       /* end element loop */
449            }           /* end color loop */
450    
451            delete[] EM_S;
452            delete[] EM_F;
453            delete[] row_index;
454    
455        } /* end of pointer check */      }           /* end of pointer check */
456     } /* end parallel region */      }               /* end parallel region */
457  }  }
458    
459  /*  /*
460   * $Log$   * $Log$
461   */   */

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

  ViewVC Help
Powered by ViewVC 1.1.26