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

temp_trunk_copy/finley/src/Assemble_PDE_System2_2D.c revision 1384 by phornby, Fri Jan 11 02:29:38 2008 UTC branches/domexper/dudley/src/Assemble_PDE_System2_2D.c revision 3203 by jfenwick, Thu Sep 23 23:51:26 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;  //    double *S = p.row_jac->BasisFunctions->S;
67      dim_t len_EM_S=p.row_NN*p.col_NN*p.numEqu*p.numComp;      const double* S = p.shapeFns;
68      dim_t len_EM_F=p.row_NN*p.numEqu;      dim_t len_EM_S = p.row_numShapesTotal * p.row_numShapesTotal * p.numEqu * p.numComp;
69        dim_t len_EM_F = p.row_numShapesTotal * p.numEqu;
70    
71      #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, 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)
72      {      {
        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);  
73    
74        } /* end of pointer check */      EM_S = THREAD_MEMALLOC(len_EM_S, double);
75     } /* end parallel region */      EM_F = THREAD_MEMALLOC(len_EM_F, double);
76        row_index = THREAD_MEMALLOC(p.row_numShapesTotal, index_t);
77    
78        if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
79        {
80    
81            for (color = elements->minColor; color <= elements->maxColor; color++)
82            {
83            /*  open loop over all elements: */
84    #pragma omp for private(e) schedule(static)
85            for (e = 0; e < elements->numElements; e++)
86            {
87                if (elements->Color[e] == color)
88                {
89    
90                A_p = getSampleDataRO(A, e);
91                B_p = getSampleDataRO(B, e);
92                C_p = getSampleDataRO(C, e);
93                D_p = getSampleDataRO(D, e);
94                X_p = getSampleDataRO(X, e);
95                Y_p = getSampleDataRO(Y, e);
96                double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
97                DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapesTotal, DIM, p.numQuadTotal, 1)]);
98                for (q = 0; q < len_EM_S; ++q)
99                    EM_S[q] = 0;
100                for (q = 0; q < len_EM_F; ++q)
101                    EM_F[q] = 0;
102                add_EM_F = FALSE;
103                add_EM_S = FALSE;
104    
105                  /**************************************************************/
106                /*   process A: */
107                  /**************************************************************/
108                A_p = getSampleDataRO(A, e);
109                if (NULL != A_p)
110                {
111                    add_EM_S = TRUE;
112                    if (extendedA)
113                    {
114                    A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuadTotal)]);
115                    for (s = 0; s < p.row_numShapes; s++)
116                    {
117                        for (r = 0; r < p.row_numShapes; r++)
118                        {
119                        for (k = 0; k < p.numEqu; k++)
120                        {
121                            for (m = 0; m < p.numComp; m++)
122                            {
123                            rtmp = 0;
124                            for (q = 0; q < p.numQuadTotal; q++)
125                            {
126                                rtmp +=
127                                vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
128                                       A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
129                                       * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
130                                       DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
131                                       A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
132                                       * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
133                                       DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
134                                       A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
135                                       * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
136                                       DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
137                                       A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
138                                       * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);
139                            }
140                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
141                                rtmp;
142                            }
143                        }
144                        }
145                    }
146                    } else
147                    {
148                    for (s = 0; s < p.row_numShapes; s++)
149                    {
150                        for (r = 0; r < p.row_numShapes; r++)
151                        {
152                        rtmp00 = 0;
153                        rtmp01 = 0;
154                        rtmp10 = 0;
155                        rtmp11 = 0;
156                        for (q = 0; q < p.numQuadTotal; q++)
157                        {
158                            rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
159                            rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
160                            rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
161                            rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
162                            rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
163                            rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
164                        }
165                        for (k = 0; k < p.numEqu; k++)
166                        {
167                            for (m = 0; m < p.numComp; m++)
168                            {
169                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
170                                rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)]
171                                + rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)]
172                                + rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)]
173                                + rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)];
174                            }
175                        }
176                        }
177                    }
178                    }
179                }
180                  /**************************************************************/
181                /*   process B: */
182                  /**************************************************************/
183                B_p = getSampleDataRO(B, e);
184                if (NULL != B_p)
185                {
186                    add_EM_S = TRUE;
187                    if (extendedB)
188                    {
189                    B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);
190                    for (s = 0; s < p.row_numShapes; s++)
191                    {
192                        for (r = 0; r < p.row_numShapes; r++)
193                        {
194                        for (k = 0; k < p.numEqu; k++)
195                        {
196                            for (m = 0; m < p.numComp; m++)
197                            {
198                            rtmp = 0;
199                            for (q = 0; q < p.numQuadTotal; q++)
200                            {
201                                rtmp += vol * S[INDEX2(r, q, p.row_numShapes)] *
202                                (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
203                                 B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
204                                 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
205                                 B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);
206                            }
207                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
208                                rtmp;
209                            }
210                        }
211                        }
212                    }
213                    } else
214                    {
215                    for (s = 0; s < p.row_numShapes; s++)
216                    {
217                        for (r = 0; r < p.row_numShapes; r++)
218                        {
219                        rtmp0 = 0;
220                        rtmp1 = 0;
221                        for (q = 0; q < p.numQuadTotal; q++)
222                        {
223                            rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];
224                            rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
225                            rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
226                        }
227                        for (k = 0; k < p.numEqu; k++)
228                        {
229                            for (m = 0; m < p.numComp; m++)
230                            {
231                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
232                                rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
233                                rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];
234                            }
235                        }
236                        }
237                    }
238                    }
239                }
240                  /**************************************************************/
241                /*   process C: */
242                  /**************************************************************/
243                C_p = getSampleDataRO(C, e);
244                if (NULL != C_p)
245                {
246                    add_EM_S = TRUE;
247                    if (extendedC)
248                    {
249                    C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);
250                    for (s = 0; s < p.row_numShapes; s++)
251                    {
252                        for (r = 0; r < p.row_numShapes; r++)
253                        {
254                        for (k = 0; k < p.numEqu; k++)
255                        {
256                            for (m = 0; m < p.numComp; m++)
257                            {
258                            rtmp = 0;
259                            for (q = 0; q < p.numQuadTotal; q++)
260                            {
261                                rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] *
262                                (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
263                                 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
264                                 C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
265                                 DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);
266                            }
267                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
268                                rtmp;
269                            }
270                        }
271                        }
272                    }
273                    } else
274                    {
275                    for (s = 0; s < p.row_numShapes; s++)
276                    {
277                        for (r = 0; r < p.row_numShapes; r++)
278                        {
279                        rtmp0 = 0;
280                        rtmp1 = 0;
281                        for (q = 0; q < p.numQuadTotal; q++)
282                        {
283                            rtmp = vol * S[INDEX2(s, q, p.row_numShapes)];
284                            rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
285                            rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
286                        }
287                        for (k = 0; k < p.numEqu; k++)
288                        {
289                            for (m = 0; m < p.numComp; m++)
290                            {
291                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
292                                rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
293                                rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];
294                            }
295                        }
296                        }
297                    }
298                    }
299                }
300                  /************************************************************* */
301                /* process D */
302                  /**************************************************************/
303                D_p = getSampleDataRO(D, e);
304                if (NULL != D_p)
305                {
306                    add_EM_S = TRUE;
307                    if (extendedD)
308                    {
309                    D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);
310                    for (s = 0; s < p.row_numShapes; s++)
311                    {
312                        for (r = 0; r < p.row_numShapes; r++)
313                        {
314                        for (k = 0; k < p.numEqu; k++)
315                        {
316                            for (m = 0; m < p.numComp; m++)
317                            {
318                            rtmp = 0;
319                            for (q = 0; q < p.numQuadTotal; q++)
320                            {
321                                rtmp +=
322                                vol * S[INDEX2(s, q, p.row_numShapes)] *
323                                D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
324                                S[INDEX2(r, q, p.row_numShapes)];
325                            }
326                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
327                                rtmp;
328                            }
329                        }
330                        }
331                    }
332                    } else
333                    {
334                    for (s = 0; s < p.row_numShapes; s++)
335                    {
336                        for (r = 0; r < p.row_numShapes; r++)
337                        {
338                        rtmp = 0;
339                        for (q = 0; q < p.numQuadTotal; q++)
340                            rtmp +=
341                            vol * S[INDEX2(s, q, p.row_numShapes)] *
342                            S[INDEX2(r, q, p.row_numShapes)];
343                        for (k = 0; k < p.numEqu; k++)
344                        {
345                            for (m = 0; m < p.numComp; m++)
346                            {
347                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
348                                rtmp * D_p[INDEX2(k, m, p.numEqu)];
349                            }
350                        }
351                        }
352                    }
353                    }
354                }
355                  /**************************************************************/
356                /*   process X: */
357                  /**************************************************************/
358                X_p = getSampleDataRO(X, e);
359                if (NULL != X_p)
360                {
361                    add_EM_F = TRUE;
362                    if (extendedX)
363                    {
364                    X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuadTotal)]);
365                    for (s = 0; s < p.row_numShapes; s++)
366                    {
367                        for (k = 0; k < p.numEqu; k++)
368                        {
369                        rtmp = 0;
370                        for (q = 0; q < p.numQuadTotal; q++)
371                        {
372                            rtmp +=
373                            vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
374                                   X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
375                                   DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
376                                   X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);
377                        }
378                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
379                        }
380                    }
381                    } else
382                    {
383                    for (s = 0; s < p.row_numShapes; s++)
384                    {
385                        rtmp0 = 0;
386                        rtmp1 = 0;
387                        for (q = 0; q < p.numQuadTotal; q++)
388                        {
389                        rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
390                        rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, 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.numQuadTotal)]);
408                    for (s = 0; s < p.row_numShapes; s++)
409                    {
410                        for (k = 0; k < p.numEqu; k++)
411                        {
412                        rtmp = 0;
413                        for (q = 0; q < p.numQuadTotal; q++)
414                            rtmp +=
415                            vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
416                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
417                        }
418                    }
419                    } else
420                    {
421                    for (s = 0; s < p.row_numShapes; s++)
422                    {
423                        rtmp = 0;
424                        for (q = 0; q < p.numQuadTotal; q++)
425                        rtmp += vol * S[INDEX2(s, q, p.row_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.row_numShapesTotal; 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.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
439                               p.row_DOF_UpperBound);
440                if (add_EM_S)
441                    Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
442                                      p.row_numShapesTotal, row_index, p.numComp, EM_S);
443    
444                }       /* end color check */
445            }       /* end element loop */
446            }           /* end color loop */
447    
448            THREAD_MEMFREE(EM_S);
449            THREAD_MEMFREE(EM_F);
450            THREAD_MEMFREE(row_index);
451    
452        }           /* end of pointer check */
453        }               /* end parallel region */
454  }  }
455    
456  /*  /*
457   * $Log$   * $Log$
458   */   */

Legend:
Removed from v.1384  
changed lines
  Added in v.3203

  ViewVC Help
Powered by ViewVC 1.1.26