/[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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26