/[escript]/branches/doubleplusgood/dudley/src/Assemble_PDE_System2_2D.cpp
ViewVC logotype

Diff of /branches/doubleplusgood/dudley/src/Assemble_PDE_System2_2D.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/finley/src/Assemble_PDE_System2_2D.c revision 798 by gross, Fri Aug 4 01:05:36 2006 UTC branches/doubleplusgood/dudley/src/Assemble_PDE_System2_2D.cpp revision 4332 by jfenwick, Thu Mar 21 04:21:14 2013 UTC
# Line 1  Line 1 
 /*  
  ************************************************************  
  *          Copyright 2006 by ACcESS MNRF                   *  
  *                                                          *  
  *              http://www.access.edu.au                    *  
  *       Primary Business: Queensland, Australia            *  
  *  Licensed under the Open Software License version 3.0    *  
  *     http://www.opensource.org/licenses/osl-3.0.php       *  
  *                                                          *  
  ************************************************************  
 */  
1    
2  /**************************************************************/  /*****************************************************************************
3    *
4    * Copyright (c) 2003-2013 by University of Queensland
5    * http://www.uq.edu.au
6    *
7    * Primary Business: Queensland, Australia
8    * Licensed under the Open Software License version 3.0
9    * http://www.opensource.org/licenses/osl-3.0.php
10    *
11    * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16    /************************************************************************************/
17    
18  /*    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  */
19  /*    the shape functions for test and solution must be identical */  /*    the shape functions for test and solution must be identical */
20    
   
21  /*      -(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 */
22    
23  /*    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 30  Line 32 
32  /*      X = p.numEqu x 2  */  /*      X = p.numEqu x 2  */
33  /*      Y = p.numEqu   */  /*      Y = p.numEqu   */
34    
35    /************************************************************************************/
 /**************************************************************/  
   
 /*  Author: gross@access.edu.au */  
 /*  Version: $Id:$ */  
   
 /**************************************************************/  
   
36    
37  #include "Assemble.h"  #include "Assemble.h"
38  #include "Util.h"  #include "Util.h"
39    #ifdef _OPENMP
40    #include <omp.h>
41    #endif
42    
43    /************************************************************************************/
44    
45    void Dudley_Assemble_PDE_System2_2D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements,
46                        Paso_SystemMatrix * Mat, escriptDataC * F,
47                        escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,
48                        escriptDataC * X, escriptDataC * Y)
49    {
50    
51  /**************************************************************/  #define DIM 2
   
 void  Finley_Assemble_PDE_System2_2D(Assemble_Parameters p, Finley_ElementFile* elements,  
                                      Paso_SystemMatrix* Mat, escriptDataC* F,  
                                      escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y) {  
   
     #define DIM 2  
52      index_t color;      index_t color;
53      dim_t e;      dim_t e;
54      bool_t extendedA=isExpanded(A);      __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;
55      bool_t extendedB=isExpanded(B);      double *EM_S, *EM_F, *DSDX;
56      bool_t extendedC=isExpanded(C);      index_t *row_index;
57      bool_t extendedD=isExpanded(D);      register dim_t q, s, r, k, m;
58      bool_t extendedX=isExpanded(X);      register double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11;
59      bool_t extendedY=isExpanded(Y);      bool_t add_EM_F, add_EM_S;
60      double *F_p=getSampleData(F,0);  
61      double *S=p.row_jac->ReferenceElement->S;      bool_t extendedA = isExpanded(A);
62      dim_t len_EM_S=p.row_NN*p.col_NN*p.numEqu*p.numComp;      bool_t extendedB = isExpanded(B);
63      dim_t len_EM_F=p.row_NN*p.numEqu;      bool_t extendedC = isExpanded(C);
64        bool_t extendedD = isExpanded(D);
65        bool_t extendedX = isExpanded(X);
66        bool_t extendedY = isExpanded(Y);
67        double *F_p = (requireWrite(F), getSampleDataRW(F, 0)); /* use comma, to get around the mixed code and declarations thing */
68        const double *S = p.shapeFns;
69        dim_t len_EM_S = p.numShapes * p.numShapes * p.numEqu * p.numComp;
70        dim_t len_EM_F = p.numShapes * p.numEqu;
71    
72      #pragma omp parallel private(color)  #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)
73      {      {
74         double EM_S[len_EM_S], EM_F[len_EM_F];  
75         index_t row_index[p.row_NN];      EM_S = new  double[len_EM_S];
76         register dim_t q, s,r,k,m;      EM_F = new  double[len_EM_F];
77         register double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11;      row_index = new  index_t[p.numShapes];
78         double *Vol, *DSDX, *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p;  
79         bool_t add_EM_F, add_EM_S;      if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
80        {
81         #ifndef PASO_MPI  
82         for (color=elements->minColor;color<=elements->maxColor;color++) {          for (color = elements->minColor; color <= elements->maxColor; color++)
83            /*  open loop over all elements: */          {
84            #pragma omp for private(e) schedule(static)          /*  open loop over all elements: */
85            for(e=0;e<elements->numElements;e++){  #pragma omp for private(e) schedule(static)
86               if (elements->Color[e]==color) {          for (e = 0; e < elements->numElements; e++)
87         #else          {
88         {              if (elements->Color[e] == color)
89            for(e=0;e<elements->numElements;e++) {              {
90               {              double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
91         #endif  
92                  Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);              A_p = getSampleDataRO(A, e);
93                  DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);              B_p = getSampleDataRO(B, e);
94                  for (q=0;q<len_EM_S;++q) EM_S[q]=0;              C_p = getSampleDataRO(C, e);
95                  for (q=0;q<len_EM_F;++q) EM_F[q]=0;              D_p = getSampleDataRO(D, e);
96                  add_EM_F=FALSE;              X_p = getSampleDataRO(X, e);
97                  add_EM_S=FALSE;              Y_p = getSampleDataRO(Y, e);
98                  /**************************************************************/              DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.numShapes, DIM, p.numQuad, 1)]);
99                  /*   process A: */              for (q = 0; q < len_EM_S; ++q)
100                  /**************************************************************/                  EM_S[q] = 0;
101                  A_p=getSampleData(A,e);              for (q = 0; q < len_EM_F; ++q)
102                  if (NULL!=A_p) {                  EM_F[q] = 0;
103                     add_EM_S=TRUE;              add_EM_F = FALSE;
104                     if (extendedA) {              add_EM_S = FALSE;
105                        for (s=0;s<p.row_NS;s++) {  
106                          for (r=0;r<p.col_NS;r++) {                /************************************************************************************/
107                            for (k=0;k<p.numEqu;k++) {              /*   process A: */
108                              for (m=0;m<p.numComp;m++) {                /************************************************************************************/
109                                rtmp=0;              A_p = getSampleDataRO(A, e);
110                                for (q=0;q<p.numQuad;q++) {              if (NULL != A_p)
111                                   rtmp+=Vol[q]* (              {
112                                       DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX5(k,0,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]                  add_EM_S = TRUE;
113                                      +DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX5(k,0,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]                  if (extendedA)
114                                      +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX5(k,1,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]                  {
115                                      +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX5(k,1,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]);                  A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuad)]);
116                                }                  for (s = 0; s < p.numShapes; s++)
117                                EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                  {
118                              }                      for (r = 0; r < p.numShapes; r++)
119                            }                      {
120                          }                      for (k = 0; k < p.numEqu; k++)
121                        }                      {
122                     } else {                          for (m = 0; m < p.numComp; m++)
123                        for (s=0;s<p.row_NS;s++) {                          {
124                          for (r=0;r<p.col_NS;r++) {                          rtmp = 0;
125                             rtmp00=0;                          for (q = 0; q < p.numQuad; q++)
126                             rtmp01=0;                          {
127                             rtmp10=0;                              rtmp +=
128                             rtmp11=0;                              vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
129                             for (q=0;q<p.numQuad;q++) {                                     A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
130                                   rtmp0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                                     * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
131                                   rtmp1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                                     DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
132                                   rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                                     A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
133                                   rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                                     * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
134                                   rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                                     DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
135                                   rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                                     A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
136                             }                                     * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
137                             for (k=0;k<p.numEqu;k++) {                                     DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
138                                 for (m=0;m<p.numComp;m++) {                                     A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
139                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=                                     * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
140                                               rtmp00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]                          }
141                                              +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
142                                              +rtmp10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]                          }
143                                              +rtmp11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)];                      }
144                                 }                      }
145                             }                  }
146                          }                  }
147                        }                  else
148                    }                  {
149                  }                  for (s = 0; s < p.numShapes; s++)
150                  /**************************************************************/                  {
151                  /*   process B: */                      for (r = 0; r < p.numShapes; r++)
152                  /**************************************************************/                      {
153                  B_p=getSampleData(B,e);                      rtmp00 = 0;
154                  if (NULL!=B_p) {                      rtmp01 = 0;
155                     add_EM_S=TRUE;                      rtmp10 = 0;
156                     if (extendedB) {                      rtmp11 = 0;
157                        for (s=0;s<p.row_NS;s++) {                      for (q = 0; q < p.numQuad; q++)
158                          for (r=0;r<p.col_NS;r++) {                      {
159                            for (k=0;k<p.numEqu;k++) {                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
160                              for (m=0;m<p.numComp;m++) {                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
161                                 rtmp=0;                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
162                                 for (q=0;q<p.numQuad;q++) {                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
163                                    rtmp+=Vol[q]*S[INDEX2(r,q,p.row_NS)]*                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
164                                         ( DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*B_p[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]                          rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
165                                         + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*B_p[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]);                      }
166                                 }                      for (k = 0; k < p.numEqu; k++)
167                                 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                      {
168                              }                          for (m = 0; m < p.numComp; m++)
169                            }                          {
170                          }                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
171                        }                              rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)]
172                     } else {                              + rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)]
173                        for (s=0;s<p.row_NS;s++) {                              + rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)]
174                          for (r=0;r<p.col_NS;r++) {                              + rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)];
175                              rtmp0=0;                          }
176                              rtmp1=0;                      }
177                              for (q=0;q<p.numQuad;q++) {                      }
178                                  rtmp=Vol[q]*S[INDEX2(r,q,p.row_NS)];                  }
179                                  rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                  }
180                                  rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];              }
181                              }                /************************************************************************************/
182                              for (k=0;k<p.numEqu;k++) {              /*   process B: */
183                                 for (m=0;m<p.numComp;m++) {                /************************************************************************************/
184                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+= rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]              B_p = getSampleDataRO(B, e);
185                                                                                      + rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)];              if (NULL != B_p)
186                                 }              {
187                              }                  add_EM_S = TRUE;
188                          }                  if (extendedB)
189                        }                  {
190                     }                  B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuad)]);
191                  }                  for (s = 0; s < p.numShapes; s++)
192                  /**************************************************************/                  {
193                  /*   process C: */                      for (r = 0; r < p.numShapes; r++)
194                  /**************************************************************/                      {
195                  C_p=getSampleData(C,e);                      for (k = 0; k < p.numEqu; k++)
196                  if (NULL!=C_p) {                      {
197                    add_EM_S=TRUE;                          for (m = 0; m < p.numComp; m++)
198                    if (extendedC) {                          {
199                        for (s=0;s<p.row_NS;s++) {                          rtmp = 0;
200                          for (r=0;r<p.col_NS;r++) {                          for (q = 0; q < p.numQuad; q++)
201                            for (k=0;k<p.numEqu;k++) {                          {
202                              for (m=0;m<p.numComp;m++) {                              rtmp += vol * S[INDEX2(r, q, p.numShapes)] *
203                                 rtmp=0;                              (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
204                                 for (q=0;q<p.numQuad;q++) {                               B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
205                                     rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*                               DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
206                                          ( C_p[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]                               B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);
207                                          + C_p[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]);                          }
208                                 }                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
209                                 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                          }
210                              }                      }
211                            }                      }
212                          }                  }
213                        }                  }
214                    } else {                  else
215                        for (s=0;s<p.row_NS;s++) {                  {
216                          for (r=0;r<p.col_NS;r++) {                  for (s = 0; s < p.numShapes; s++)
217                             rtmp0=0;                  {
218                             rtmp1=0;                      for (r = 0; r < p.numShapes; r++)
219                             for (q=0;q<p.numQuad;q++) {                      {
220                                     rtmp=Vol[q]*S[INDEX2(s,q,p.row_NS)];                      rtmp0 = 0;
221                                     rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                      rtmp1 = 0;
222                                     rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                      for (q = 0; q < p.numQuad; q++)
223                             }                      {
224                             for (k=0;k<p.numEqu;k++) {                          rtmp = vol * S[INDEX2(r, q, p.numShapes)];
225                                for (m=0;m<p.numComp;m++) {                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
226                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp0*C_p[INDEX3(k,m,0,p.numEqu,p.numComp)]                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
227                                                                                         +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)];                      }
228                                }                      for (k = 0; k < p.numEqu; k++)
229                             }                      {
230                          }                          for (m = 0; m < p.numComp; m++)
231                        }                          {
232                    }                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
233                  }                              rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
234                  /************************************************************* */                              rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];
235                  /* process D */                          }
236                  /**************************************************************/                      }
237                  D_p=getSampleData(D,e);                      }
238                  if (NULL!=D_p) {                  }
239                    add_EM_S=TRUE;                  }
240                    if (extendedD) {              }
241                        for (s=0;s<p.row_NS;s++) {                /************************************************************************************/
242                          for (r=0;r<p.col_NS;r++) {              /*   process C: */
243                            for (k=0;k<p.numEqu;k++) {                /************************************************************************************/
244                              for (m=0;m<p.numComp;m++) {              C_p = getSampleDataRO(C, e);
245                                rtmp=0;              if (NULL != C_p)
246                                for (q=0;q<p.numQuad;q++) {              {
247                                    rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_NS)];                  add_EM_S = TRUE;
248                                }                  if (extendedC)
249                                EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                  {
250                              }                  C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuad)]);
251                            }                  for (s = 0; s < p.numShapes; s++)
252                          }                  {
253                        }                      for (r = 0; r < p.numShapes; r++)
254                    } else {                      {
255                        for (s=0;s<p.row_NS;s++) {                      for (k = 0; k < p.numEqu; k++)
256                          for (r=0;r<p.col_NS;r++) {                      {
257                              rtmp=0;                          for (m = 0; m < p.numComp; m++)
258                              for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];                          {
259                              for (k=0;k<p.numEqu;k++) {                          rtmp = 0;
260                                  for (m=0;m<p.numComp;m++) {                          for (q = 0; q < p.numQuad; q++)
261                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];                          {
262                                 }                              rtmp += vol * S[INDEX2(s, q, p.numShapes)] *
263                              }                              (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
264                          }                               DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
265                        }                               C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
266                    }                               DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
267                  }                          }
268                  /**************************************************************/                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
269                  /*   process X: */                          }
270                  /**************************************************************/                      }
271                  X_p=getSampleData(X,e);                      }
272                  if (NULL!=X_p) {                  }
273                    add_EM_F=TRUE;                  }
274                    if (extendedX) {                  else
275                       for (s=0;s<p.row_NS;s++) {                  {
276                          for (k=0;k<p.numEqu;k++) {                  for (s = 0; s < p.numShapes; s++)
277                            rtmp=0;                  {
278                            for (q=0;q<p.numQuad;q++) {                      for (r = 0; r < p.numShapes; r++)
279                                 rtmp+=Vol[q]*(DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX3(k,0,q,p.numEqu,DIM)]                      {
280                                              +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*X_p[INDEX3(k,1,q,p.numEqu,DIM)]);                      rtmp0 = 0;
281                            }                      rtmp1 = 0;
282                            EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;                      for (q = 0; q < p.numQuad; q++)
283                          }                      {
284                       }                          rtmp = vol * S[INDEX2(s, q, p.numShapes)];
285                    } else {                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
286                       for (s=0;s<p.row_NS;s++) {                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
287                             rtmp0=0;                      }
288                             rtmp1=0;                      for (k = 0; k < p.numEqu; k++)
289                             for (q=0;q<p.numQuad;q++) {                      {
290                                 rtmp0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                          for (m = 0; m < p.numComp; m++)
291                                 rtmp1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                          {
292                             }                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
293                             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)];                              rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
294                       }                              rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];
295                    }                          }
296                 }                      }
297                 /**************************************************************/                      }
298                 /*   process Y: */                  }
299                 /**************************************************************/                  }
300                  Y_p=getSampleData(Y,e);              }
301                  if (NULL!=Y_p) {                /*********************************************************************************** */
302                    add_EM_F=TRUE;              /* process D */
303                    if (extendedY) {                /************************************************************************************/
304                       for (s=0;s<p.row_NS;s++) {              D_p = getSampleDataRO(D, e);
305                          for (k=0;k<p.numEqu;k++) {              if (NULL != D_p)
306                             rtmp=0;              {
307                             for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[INDEX2(k,q,p.numEqu)];                  add_EM_S = TRUE;
308                             EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;                  if (extendedD)
309                          }                  {
310                       }                  D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuad)]);
311                     } else {                  for (s = 0; s < p.numShapes; s++)
312                       for (s=0;s<p.row_NS;s++) {                  {
313                           rtmp=0;                      for (r = 0; r < p.numShapes; r++)
314                           for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];                      {
315                           for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];                      for (k = 0; k < p.numEqu; k++)
316                       }                      {
317                     }                          for (m = 0; m < p.numComp; m++)
318                   }                          {
319                   /***********************************************************************************************/                          rtmp = 0;
320                   /* add the element matrices onto the matrix and right hand side                                */                          for (q = 0; q < p.numQuad; q++)
321                   /***********************************************************************************************/                          {
322                   for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];                              rtmp +=
323                   if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);                              vol * S[INDEX2(s, q, p.numShapes)] *
324                   if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);                              D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
325                                S[INDEX2(r, q, p.numShapes)];
326               } /* end color check */                          }
327            } /* end element loop */                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
328        } /* end color loop */                          }
329     } /* end parallel region */                      }
330                        }
331                    }
332                    }
333                    else
334                    {
335                    for (s = 0; s < p.numShapes; s++)
336                    {
337                        for (r = 0; r < p.numShapes; r++)
338                        {
339                        rtmp = 0;
340                        for (q = 0; q < p.numQuad; q++)
341                            rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
342                        for (k = 0; k < p.numEqu; k++)
343                        {
344                            for (m = 0; m < p.numComp; m++)
345                            {
346                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
347                                rtmp * D_p[INDEX2(k, m, p.numEqu)];
348                            }
349                        }
350                        }
351                    }
352                    }
353                }
354                  /************************************************************************************/
355                /*   process X: */
356                  /************************************************************************************/
357                X_p = getSampleDataRO(X, e);
358                if (NULL != X_p)
359                {
360                    add_EM_F = TRUE;
361                    if (extendedX)
362                    {
363                    X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuad)]);
364                    for (s = 0; s < p.numShapes; s++)
365                    {
366                        for (k = 0; k < p.numEqu; k++)
367                        {
368                        rtmp = 0;
369                        for (q = 0; q < p.numQuad; q++)
370                        {
371                            rtmp +=
372                            vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
373                                   X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
374                                   DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
375                                   X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);
376                        }
377                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
378                        }
379                    }
380                    }
381                    else
382                    {
383                    for (s = 0; s < p.numShapes; s++)
384                    {
385                        rtmp0 = 0;
386                        rtmp1 = 0;
387                        for (q = 0; q < p.numQuad; q++)
388                        {
389                        rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
390                        rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
391                        }
392                        for (k = 0; k < p.numEqu; k++)
393                        EM_F[INDEX2(k, s, p.numEqu)] +=
394                            rtmp0 * X_p[INDEX2(k, 0, p.numEqu)] + rtmp1 * X_p[INDEX2(k, 1, p.numEqu)];
395                    }
396                    }
397                }
398                 /************************************************************************************/
399                /*   process Y: */
400                 /************************************************************************************/
401                Y_p = getSampleDataRO(Y, e);
402                if (NULL != Y_p)
403                {
404                    add_EM_F = TRUE;
405                    if (extendedY)
406                    {
407                    Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuad)]);
408                    for (s = 0; s < p.numShapes; s++)
409                    {
410                        for (k = 0; k < p.numEqu; k++)
411                        {
412                        rtmp = 0;
413                        for (q = 0; q < p.numQuad; q++)
414                            rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
415                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
416                        }
417                    }
418                    }
419                    else
420                    {
421                    for (s = 0; s < p.numShapes; s++)
422                    {
423                        rtmp = 0;
424                        for (q = 0; q < p.numQuad; q++)
425                        rtmp += vol * S[INDEX2(s, q, p.numShapes)];
426                        for (k = 0; k < p.numEqu; k++)
427                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
428                    }
429                    }
430                }
431                   /*********************************************************************************************************************/
432                /* add the element matrices onto the matrix and right hand side                                */
433                   /*********************************************************************************************************************/
434                for (q = 0; q < p.numShapes; q++)
435                    row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
436    
437                if (add_EM_F)
438                    Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound);
439                if (add_EM_S)
440                    Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
441                                      p.numShapes, row_index, p.numComp, EM_S);
442    
443                }       /* end color check */
444            }       /* end element loop */
445            }           /* end color loop */
446    
447            delete[] EM_S;
448            delete[] EM_F;
449            delete[] row_index;
450    
451        }           /* end of pointer check */
452        }               /* end parallel region */
453  }  }
454    
455  /*  /*
456   * $Log$   * $Log$
457   */   */

Legend:
Removed from v.798  
changed lines
  Added in v.4332

  ViewVC Help
Powered by ViewVC 1.1.26