/[escript]/branches/domexper/dudley/src/Assemble_PDE_System2_3D.c
ViewVC logotype

Diff of /branches/domexper/dudley/src/Assemble_PDE_System2_3D.c

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

trunk/finley/src/Assemble_PDE_System2_3D.c revision 853 by gross, Wed Sep 20 05:56:36 2006 UTC branches/domexper/dudley/src/Assemble_PDE_System2_3D.c revision 3202 by jfenwick, Thu Sep 23 23:24:09 2010 UTC
# Line 1  Line 1 
1  /*  
2   ************************************************************  /*******************************************************
3   *          Copyright 2006 by ACcESS MNRF                   *  *
4   *                                                          *  * Copyright (c) 2003-2010 by University of Queensland
5   *              http://www.access.edu.au                    *  * Earth Systems Science Computational Center (ESSCC)
6   *       Primary Business: Queensland, Australia            *  * http://www.uq.edu.au/esscc
7   *  Licensed under the Open Software License version 3.0    *  *
8   *     http://www.opensource.org/licenses/osl-3.0.php       *  * Primary Business: Queensland, Australia
9   *                                                          *  * Licensed under the Open Software License version 3.0
10   ************************************************************  * 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 3D domain. The shape functions for test and solution must be identical  */  /*    u has p.numComp components in a 3D domain. The shape functions for test and solution must be identical  */
# Line 30  Line 30 
30  /*      X = p.numEqu x 3  */  /*      X = p.numEqu x 3  */
31  /*      Y = p.numEqu   */  /*      Y = p.numEqu   */
32    
   
33  /**************************************************************/  /**************************************************************/
34    
 /*  Author: gross@access.edu.au */  
 /*  Version: $Id:$ */  
   
 /**************************************************************/  
   
   
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_3D(Assemble_Parameters p, Finley_ElementFile* elements,  void Dudley_Assemble_PDE_System2_3D(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 3  #define DIM 3
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, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22;      register double rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22;
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.col_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, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22,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, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22,add_EM_F, add_EM_S)
72      {      {
73         EM_S=THREAD_MEMALLOC(len_EM_S,double);      EM_S = THREAD_MEMALLOC(len_EM_S, double);
74         EM_F=THREAD_MEMALLOC(len_EM_F,double);      EM_F = THREAD_MEMALLOC(len_EM_F, double);
75         row_index=THREAD_MEMALLOC(p.row_NN,index_t);      row_index = THREAD_MEMALLOC(p.row_numShapesTotal, index_t);
76                                                                                                                                                                                                        
77         if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) {      if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
78        {
79            #ifndef PASO_MPI  
80            for (color=elements->minColor;color<=elements->maxColor;color++) {          for (color = elements->minColor; color <= elements->maxColor; color++)
81               /*  open loop over all elements: */          {
82               #pragma omp for private(e) schedule(static)          /*  open loop over all elements: */
83               for(e=0;e<elements->numElements;e++){  #pragma omp for private(e) schedule(static)
84                  if (elements->Color[e]==color) {          for (e = 0; e < elements->numElements; e++)
85            #else          {
86            {              if (elements->Color[e] == color)
87               for(e=0;e<elements->numElements;e++) {              {
88                  {  
89            #endif              A_p = getSampleDataRO(A, e);
90                     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);              B_p = getSampleDataRO(B, e);
91                     DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);              C_p = getSampleDataRO(C, e);
92                     for (q=0;q<len_EM_S;++q) EM_S[q]=0;              D_p = getSampleDataRO(D, e);
93                     for (q=0;q<len_EM_F;++q) EM_F[q]=0;              X_p = getSampleDataRO(X, e);
94                     add_EM_F=FALSE;              Y_p = getSampleDataRO(Y, e);
95                     add_EM_S=FALSE;              double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
96                     /**************************************************************/  
97                     /*   process A: */              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                     A_p=getSampleData(A,e);                  EM_S[q] = 0;
100                     if (NULL!=A_p) {              for (q = 0; q < len_EM_F; ++q)
101                        add_EM_S=TRUE;                  EM_F[q] = 0;
102                        if (extendedA) {              add_EM_F = FALSE;
103                           for (s=0;s<p.row_NS;s++) {              add_EM_S = FALSE;
104                             for (r=0;r<p.col_NS;r++) {  
105                               for (k=0;k<p.numEqu;k++) {                /**************************************************************/
106                                 for (m=0;m<p.numComp;m++) {              /*   process A: */
107                                    rtmp=0;                /**************************************************************/
108                                    for (q=0;q<p.numQuad;q++) {              if (NULL != A_p)
109                                       rtmp+=Vol[q]*(              {
110                                         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;
111                                        +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)
112                                        +DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX5(k,0,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]                  {
113                                        +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)]                  A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuadTotal)]);
114                                        +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)]                  for (s = 0; s < p.row_numShapes; s++)
115                                        +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX5(k,1,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]                  {
116                                        +DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*A_p[INDEX5(k,2,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]                      for (r = 0; r < p.col_numShapes; r++)
117                                        +DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*A_p[INDEX5(k,2,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]                      {
118                                        +DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*A_p[INDEX5(k,2,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]);                      for (k = 0; k < p.numEqu; k++)
119                                                  {
120                                    }                          for (m = 0; m < p.numComp; m++)
121                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                          {
122                                 }                          rtmp = 0;
123                               }                          for (q = 0; q < p.numQuadTotal; q++)
124                             }                          {
125                           }                              rtmp +=
126                        } else {                              vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
127                           for (s=0;s<p.row_NS;s++) {                                     A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
128                             for (r=0;r<p.col_NS;r++) {                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
129                                 rtmp00=0;                                     DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
130                                 rtmp01=0;                                     A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
131                                 rtmp02=0;                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
132                                 rtmp10=0;                                     DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
133                                 rtmp11=0;                                     A_q[INDEX5(k, 0, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
134                                 rtmp12=0;                                     * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)] +
135                                 rtmp20=0;                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
136                                 rtmp21=0;                                     A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
137                                 rtmp22=0;                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
138                                 for (q=0;q<p.numQuad;q++) {                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
139                                       rtmp0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                                     A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
140                                       rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
141                                       rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                                     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
142                                       rtmp02+=rtmp0*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];                                     A_q[INDEX5(k, 1, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
143                                         * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)] +
144                                       rtmp1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                                     DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
145                                       rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                                     A_q[INDEX5(k, 2, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
146                                       rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                                     * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
147                                       rtmp12+=rtmp1*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];                                     DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
148                                         A_q[INDEX5(k, 2, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
149                                       rtmp2=Vol[q]*DSDX[INDEX3(s,2,q,p.row_NN,DIM)];                                     * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
150                                       rtmp20+=rtmp2*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                                     DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
151                                       rtmp21+=rtmp2*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                                     A_q[INDEX5(k, 2, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
152                                       rtmp22+=rtmp2*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];                                     * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)]);
153                                 }  
154                                 for (k=0;k<p.numEqu;k++) {                          }
155                                    for (m=0;m<p.numComp;m++) {                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
156                                       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)]                              rtmp;
157                                                                                          +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]                          }
158                                                                                          +rtmp02*A_p[INDEX4(k,0,m,2,p.numEqu,DIM,p.numComp)]                      }
159                                                                                          +rtmp10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]                      }
160                                                                                          +rtmp11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)]                  }
161                                                                                          +rtmp12*A_p[INDEX4(k,1,m,2,p.numEqu,DIM,p.numComp)]                  } else
162                                                                                          +rtmp20*A_p[INDEX4(k,2,m,0,p.numEqu,DIM,p.numComp)]                  {
163                                                                                          +rtmp21*A_p[INDEX4(k,2,m,1,p.numEqu,DIM,p.numComp)]                  for (s = 0; s < p.row_numShapes; s++)
164                                                                                          +rtmp22*A_p[INDEX4(k,2,m,2,p.numEqu,DIM,p.numComp)];                  {
165                                    }                      for (r = 0; r < p.col_numShapes; r++)
166                                 }                      {
167                             }                      rtmp00 = 0;
168                           }                      rtmp01 = 0;
169                         }                      rtmp02 = 0;
170                     }                      rtmp10 = 0;
171                     /**************************************************************/                      rtmp11 = 0;
172                     /*   process B: */                      rtmp12 = 0;
173                     /**************************************************************/                      rtmp20 = 0;
174                     B_p=getSampleData(B,e);                      rtmp21 = 0;
175                     if (NULL!=B_p) {                      rtmp22 = 0;
176                        add_EM_S=TRUE;                      for (q = 0; q < p.numQuadTotal; q++)
177                        if (extendedB) {                      {
178                           for (s=0;s<p.row_NS;s++) {                          rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
179                             for (r=0;r<p.col_NS;r++) {                          rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
180                               for (k=0;k<p.numEqu;k++) {                          rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
181                                 for (m=0;m<p.numComp;m++) {                          rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
182                                    rtmp=0;  
183                                    for (q=0;q<p.numQuad;q++) {                          rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
184                                        rtmp+=Vol[q]*S[INDEX2(r,q,p.row_NS)]*(                          rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
185                                                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.row_numShapesTotal, DIM)];
186                                              + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*B_p[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]                          rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
187                                              + DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*B_p[INDEX4(k,2,m,q,p.numEqu,DIM,p.numComp)] );  
188                                    }                          rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
189                                    EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                          rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
190                                 }                          rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
191                               }                          rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
192                             }                      }
193                           }                      for (k = 0; k < p.numEqu; k++)
194                        } else {                      {
195                           for (s=0;s<p.row_NS;s++) {                          for (m = 0; m < p.numComp; m++)
196                             for (r=0;r<p.col_NS;r++) {                          {
197                                    rtmp0=0;                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
198                                    rtmp1=0;                              rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)] +
199                                    rtmp2=0;                              rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)] +
200                                    for (q=0;q<p.numQuad;q++) {                              rtmp02 * A_p[INDEX4(k, 0, m, 2, p.numEqu, DIM, p.numComp)] +
201                                         rtmp=Vol[q]*S[INDEX2(r,q,p.row_NS)];                              rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)] +
202                                         rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                              rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)] +
203                                         rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                              rtmp12 * A_p[INDEX4(k, 1, m, 2, p.numEqu, DIM, p.numComp)] +
204                                         rtmp2+=rtmp*DSDX[INDEX3(s,2,q,p.row_NN,DIM)];                              rtmp20 * A_p[INDEX4(k, 2, m, 0, p.numEqu, DIM, p.numComp)] +
205                                    }                              rtmp21 * A_p[INDEX4(k, 2, m, 1, p.numEqu, DIM, p.numComp)] +
206                                    for (k=0;k<p.numEqu;k++) {                              rtmp22 * A_p[INDEX4(k, 2, m, 2, p.numEqu, DIM, p.numComp)];
207                                       for (m=0;m<p.numComp;m++) {                          }
208                                          EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]                      }
209                                                                                            +rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)]                      }
210                                                                                            +rtmp2*B_p[INDEX3(k,2,m,p.numEqu,DIM)];                  }
211                                       }                  }
212                                    }              }
213                             }                /**************************************************************/
214                           }              /*   process B: */
215                        }                /**************************************************************/
216                     }              if (NULL != B_p)
217                     /**************************************************************/              {
218                     /*   process C: */                  add_EM_S = TRUE;
219                     /**************************************************************/                  if (extendedB)
220                     C_p=getSampleData(C,e);                  {
221                     if (NULL!=C_p) {                  B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);
222                       add_EM_S=TRUE;                  for (s = 0; s < p.row_numShapes; s++)
223                       if (extendedC) {                  {
224                           for (s=0;s<p.row_NS;s++) {                      for (r = 0; r < p.col_numShapes; r++)
225                             for (r=0;r<p.col_NS;r++) {                      {
226                               for (k=0;k<p.numEqu;k++) {                      for (k = 0; k < p.numEqu; k++)
227                                 for (m=0;m<p.numComp;m++) {                      {
228                                   rtmp=0;                          for (m = 0; m < p.numComp; m++)
229                                   for (q=0;q<p.numQuad;q++) {                          {
230                                        rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*(                          rtmp = 0;
231                                                   C_p[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]                          for (q = 0; q < p.numQuadTotal; q++)
232                                                  +C_p[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]                          {
233                                                  +C_p[INDEX4(k,m,2,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]);                              rtmp +=
234                                   }                              vol * S[INDEX2(r, q, p.row_numShapes)] *
235                                   EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                              (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
236                                 }                               B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
237                               }                               DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
238                             }                               B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)] +
239                           }                               DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
240                       } else {                               B_q[INDEX4(k, 2, m, q, p.numEqu, DIM, p.numComp)]);
241                           for (s=0;s<p.row_NS;s++) {                          }
242                             for (r=0;r<p.col_NS;r++) {                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
243                                    rtmp0=0;                              rtmp;
244                                    rtmp1=0;                          }
245                                    rtmp2=0;                      }
246                                    for (q=0;q<p.numQuad;q++) {                      }
247                                       rtmp=Vol[q]*S[INDEX2(s,q,p.row_NS)];                  }
248                                       rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                  } else
249                                       rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];                  {
250                                       rtmp2+=rtmp*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];                  for (s = 0; s < p.row_numShapes; s++)
251                                    }                  {
252                                    for (k=0;k<p.numEqu;k++) {                      for (r = 0; r < p.col_numShapes; r++)
253                                       for (m=0;m<p.numComp;m++) {                      {
254                                             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)]                      rtmp0 = 0;
255                                                                                               +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)]                      rtmp1 = 0;
256                                                                                               +rtmp2*C_p[INDEX3(k,m,2,p.numEqu,p.numComp)];                      rtmp2 = 0;
257                                        }                      for (q = 0; q < p.numQuadTotal; q++)
258                                    }                      {
259                             }                          rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];
260                           }                          rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
261                       }                          rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
262                     }                          rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
263                     /************************************************************* */                      }
264                     /* process D */                      for (k = 0; k < p.numEqu; k++)
265                     /**************************************************************/                      {
266                     D_p=getSampleData(D,e);                          for (m = 0; m < p.numComp; m++)
267                     if (NULL!=D_p) {                          {
268                       add_EM_S=TRUE;                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
269                       if (extendedD) {                              rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
270                           for (s=0;s<p.row_NS;s++) {                              rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)] +
271                             for (r=0;r<p.col_NS;r++) {                              rtmp2 * B_p[INDEX3(k, 2, m, p.numEqu, DIM)];
272                               for (k=0;k<p.numEqu;k++) {                          }
273                                 for (m=0;m<p.numComp;m++) {                      }
274                                   rtmp=0;                      }
275                                   for (q=0;q<p.numQuad;q++) {                  }
276                                       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)];                  }
277                                   }              }
278                                   EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                /**************************************************************/
279                                 }              /*   process C: */
280                               }                /**************************************************************/
281                             }              if (NULL != C_p)
282                           }              {
283                       } else {                  add_EM_S = TRUE;
284                           for (s=0;s<p.row_NS;s++) {                  if (extendedC)
285                             for (r=0;r<p.col_NS;r++) {                  {
286                                 rtmp=0;                  C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);
287                                 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 (s = 0; s < p.row_numShapes; s++)
288                                 for (k=0;k<p.numEqu;k++) {                  {
289                                     for (m=0;m<p.numComp;m++) {                      for (r = 0; r < p.col_numShapes; r++)
290                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];                      {
291                                    }                      for (k = 0; k < p.numEqu; k++)
292                                 }                      {
293                             }                          for (m = 0; m < p.numComp; m++)
294                           }                          {
295                       }                          rtmp = 0;
296                     }                          for (q = 0; q < p.numQuadTotal; q++)
297                     /**************************************************************/                          {
298                     /*   process X: */                              rtmp +=
299                     /**************************************************************/                              vol * S[INDEX2(s, q, p.row_numShapes)] *
300                     X_p=getSampleData(X,e);                              (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
301                     if (NULL!=X_p) {                               DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
302                       add_EM_F=TRUE;                               C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
303                       if (extendedX) {                               DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
304                          for (s=0;s<p.row_NS;s++) {                               C_q[INDEX4(k, m, 2, q, p.numEqu, p.numComp, DIM)] *
305                             for (k=0;k<p.numEqu;k++) {                               DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)]);
306                                rtmp=0;                          }
307                                for (q=0;q<p.numQuad;q++) {                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
308                                      rtmp+=Vol[q]* ( DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX3(k,0,q,p.numEqu,DIM)]                              rtmp;
309                                                    + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*X_p[INDEX3(k,1,q,p.numEqu,DIM)]                          }
310                                                    + DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*X_p[INDEX3(k,2,q,p.numEqu,DIM)]);                      }
311                                }                      }
312                                EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;                  }
313                             }                  } else
314                          }                  {
315                       } else {                  for (s = 0; s < p.row_numShapes; s++)
316                          for (s=0;s<p.row_NS;s++) {                  {
317                            rtmp0=0;                      for (r = 0; r < p.col_numShapes; r++)
318                            rtmp1=0;                      {
319                            rtmp2=0;                      rtmp0 = 0;
320                            for (q=0;q<p.numQuad;q++) {                      rtmp1 = 0;
321                               rtmp0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                      rtmp2 = 0;
322                               rtmp1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];                      for (q = 0; q < p.numQuadTotal; q++)
323                               rtmp2+=Vol[q]*DSDX[INDEX3(s,2,q,p.row_NN,DIM)];                      {
324                            }                          rtmp = vol * S[INDEX2(s, q, p.row_numShapes)];
325                            for (k=0;k<p.numEqu;k++) {                          rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
326                                      EM_F[INDEX2(k,s,p.numEqu)]+=rtmp0*X_p[INDEX2(k,0,p.numEqu)]                          rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
327                                                                 +rtmp1*X_p[INDEX2(k,1,p.numEqu)]                          rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
328                                                                 +rtmp2*X_p[INDEX2(k,2,p.numEqu)];                      }
329                            }                      for (k = 0; k < p.numEqu; k++)
330                          }                      {
331                       }                          for (m = 0; m < p.numComp; m++)
332                    }                          {
333                    /**************************************************************/                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
334                    /*   process Y: */                              rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
335                    /**************************************************************/                              rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)] +
336                     Y_p=getSampleData(Y,e);                              rtmp2 * C_p[INDEX3(k, m, 2, p.numEqu, p.numComp)];
337                     if (NULL!=Y_p) {                          }
338                       add_EM_F=TRUE;                      }
339                       if (extendedY) {                      }
340                          for (s=0;s<p.row_NS;s++) {                  }
341                             for (k=0;k<p.numEqu;k++) {                  }
342                                rtmp=0.;              }
343                                for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[INDEX2(k,q,p.numEqu)];                /************************************************************* */
344                                EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;              /* process D */
345                             }                /**************************************************************/
346                          }              if (NULL != D_p)
347                        } else {              {
348                          for (s=0;s<p.row_NS;s++) {                  add_EM_S = TRUE;
349                              rtmp=0;                  if (extendedD)
350                              for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];                  {
351                              for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];                  D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);
352                          }                  for (s = 0; s < p.row_numShapes; s++)
353                        }                  {
354                      }                      for (r = 0; r < p.col_numShapes; r++)
355                      /***********************************************************************************************/                      {
356                      /* add the element matrices onto the matrix and right hand side                                */                      for (k = 0; k < p.numEqu; k++)
357                      /***********************************************************************************************/                      {
358                      for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];                          for (m = 0; m < p.numComp; m++)
359                      if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);                          {
360                      if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);                          rtmp = 0;
361                              for (q = 0; q < p.numQuadTotal; q++)
362                  } /* end color check */                          {
363               } /* end element loop */                              rtmp +=
364           } /* end color loop */                              vol * S[INDEX2(s, q, p.row_numShapes)] *
365                                          D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
366           THREAD_MEMFREE(EM_S);                              S[INDEX2(r, q, p.row_numShapes)];
367           THREAD_MEMFREE(EM_F);                          }
368           THREAD_MEMFREE(row_index);                          EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
369                                rtmp;
370                            }
371                        }
372                        }
373                    }
374                    } else
375                    {
376                    for (s = 0; s < p.row_numShapes; s++)
377                    {
378                        for (r = 0; r < p.col_numShapes; r++)
379                        {
380                        rtmp = 0;
381                        for (q = 0; q < p.numQuadTotal; q++)
382                            rtmp +=
383                            vol * S[INDEX2(s, q, p.row_numShapes)] *
384                            S[INDEX2(r, q, p.row_numShapes)];
385                        for (k = 0; k < p.numEqu; k++)
386                        {
387                            for (m = 0; m < p.numComp; m++)
388                            {
389                            EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
390                                rtmp * D_p[INDEX2(k, m, p.numEqu)];
391                            }
392                        }
393                        }
394                    }
395                    }
396                }
397                  /**************************************************************/
398                /*   process X: */
399                  /**************************************************************/
400                if (NULL != X_p)
401                {
402                    add_EM_F = TRUE;
403                    if (extendedX)
404                    {
405                    X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuadTotal)]);
406                    for (s = 0; s < p.row_numShapes; s++)
407                    {
408                        for (k = 0; k < p.numEqu; k++)
409                        {
410                        rtmp = 0;
411                        for (q = 0; q < p.numQuadTotal; q++)
412                        {
413                            rtmp +=
414                            vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
415                                   X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
416                                   DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
417                                   X_q[INDEX3(k, 1, q, p.numEqu, DIM)] +
418                                   DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
419                                   X_q[INDEX3(k, 2, q, p.numEqu, DIM)]);
420                        }
421                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
422                        }
423                    }
424                    } else
425                    {
426                    for (s = 0; s < p.row_numShapes; s++)
427                    {
428                        rtmp0 = 0;
429                        rtmp1 = 0;
430                        rtmp2 = 0;
431                        for (q = 0; q < p.numQuadTotal; q++)
432                        {
433                        rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
434                        rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
435                        rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
436                        }
437                        for (k = 0; k < p.numEqu; k++)
438                        {
439                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp0 * X_p[INDEX2(k, 0, p.numEqu)]
440                            + rtmp1 * X_p[INDEX2(k, 1, p.numEqu)] + rtmp2 * X_p[INDEX2(k, 2, p.numEqu)];
441                        }
442                    }
443                    }
444                }
445                 /**************************************************************/
446                /*   process Y: */
447                 /**************************************************************/
448                if (NULL != Y_p)
449                {
450                    add_EM_F = TRUE;
451                    if (extendedY)
452                    {
453                    Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);
454                    for (s = 0; s < p.row_numShapes; s++)
455                    {
456                        for (k = 0; k < p.numEqu; k++)
457                        {
458                        rtmp = 0.;
459                        for (q = 0; q < p.numQuadTotal; q++)
460                            rtmp +=
461                            vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
462                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
463                        }
464                    }
465                    } else
466                    {
467                    for (s = 0; s < p.row_numShapes; s++)
468                    {
469                        rtmp = 0;
470                        for (q = 0; q < p.numQuadTotal; q++)
471                        rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
472                        for (k = 0; k < p.numEqu; k++)
473                        EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
474                    }
475                    }
476                }
477    
478                   /***********************************************************************************************/
479                /* add the element matrices onto the matrix and right hand side                                */
480                   /***********************************************************************************************/
481                for (q = 0; q < p.row_numShapesTotal; q++)
482                    row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
483    
484                if (add_EM_F)
485                    Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
486                               p.row_DOF_UpperBound);
487                if (add_EM_S)
488                    Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
489                                      p.col_numShapesTotal, row_index, p.numComp, EM_S);
490                }       /* end color check */
491            }       /* end element loop */
492            }           /* end color loop */
493    
494            THREAD_MEMFREE(EM_S);
495            THREAD_MEMFREE(EM_F);
496            THREAD_MEMFREE(row_index);
497    
498        } /* end of pointer check */      }           /* end of pointer check */
499     } /* end parallel region */      }               /* end parallel region */
500  }  }
 /*  
  * $Log$  
  */  

Legend:
Removed from v.853  
changed lines
  Added in v.3202

  ViewVC Help
Powered by ViewVC 1.1.26