/[escript]/trunk/finley/src/Assemble_PDE_System2_C.c
ViewVC logotype

Diff of /trunk/finley/src/Assemble_PDE_System2_C.c

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

trunk/finley/src/Assemble_PDE_System2_C.c revision 798 by gross, Fri Aug 4 01:05:36 2006 UTC temp/finley/src/Assemble_PDE_System2_C.c revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC
# Line 1  Line 1 
1  /*  
2   ************************************************************  /* $Id$ */
3   *          Copyright 2006 by ACcESS MNRF                   *  
4   *                                                          *  /*******************************************************
5   *              http://www.access.edu.au                    *   *
6   *       Primary Business: Queensland, Australia            *   *           Copyright 2003-2007 by ACceSS MNRF
7   *  Licensed under the Open Software License version 3.0    *   *       Copyright 2007 by University of Queensland
8   *     http://www.opensource.org/licenses/osl-3.0.php       *   *
9   *                                                          *   *                http://esscc.uq.edu.au
10   ************************************************************   *        Primary Business: Queensland, Australia
11  */   *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
# Line 35  Line 38 
38    
39  #include "Assemble.h"  #include "Assemble.h"
40  #include "Util.h"  #include "Util.h"
41    #ifdef _OPENMP
42    #include <omp.h>
43    #endif
44    
45    
46  /**************************************************************/  /**************************************************************/
47    
# Line 43  void  Finley_Assemble_PDE_System2_C(Asse Line 50  void  Finley_Assemble_PDE_System2_C(Asse
50    
51      index_t color;      index_t color;
52      dim_t e;      dim_t e;
53        double *EM_S, *EM_F, *Vol, *D_p, *Y_p;
54        index_t *row_index;
55        register dim_t q, s,r,k,m;
56        register double rtmp, rtmp_D;
57        bool_t add_EM_F, add_EM_S;
58    
59      bool_t extendedD=isExpanded(D);      bool_t extendedD=isExpanded(D);
60      bool_t extendedY=isExpanded(Y);      bool_t extendedY=isExpanded(Y);
61      double *F_p=getSampleData(F,0);      double *F_p=getSampleData(F,0);
62      double *S=p.row_jac->ReferenceElement->S;      double *S=p.row_jac->ReferenceElement->S;
63    
64    
65      #pragma omp parallel private(color)      #pragma omp parallel private(color,EM_S, EM_F, Vol, D_p, Y_p,row_index,q, s,r,k,m,rtmp, rtmp_D,add_EM_F, add_EM_S)
66      {      {
67         double EM_S[p.row_NN*p.col_NN*p.numEqu*p.numComp], EM_F[p.row_NN*p.numEqu];         EM_S=THREAD_MEMALLOC(p.row_NN*p.col_NN*p.numEqu*p.numComp,double);
68         index_t row_index[p.row_NN];         EM_F=THREAD_MEMALLOC(p.row_NN*p.numEqu,double);
69         register dim_t q, s,r,k,m;         row_index=THREAD_MEMALLOC(p.row_NN,index_t);
70         register double rtmp, rtmp_D;                                                                                                                                                                                                      
71         double *Vol, *D_p, *Y_p;         if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) {
72         bool_t add_EM_F, add_EM_S;  
73         #ifndef PASO_MPI            for (color=elements->minColor;color<=elements->maxColor;color++) {
74         for (color=elements->minColor;color<=elements->maxColor;color++) {               /*  open loop over all elements: */
75            /*  open loop over all elements: */               #pragma omp for private(e) schedule(static)
76            #pragma omp for private(e) schedule(static)               for(e=0;e<elements->numElements;e++){
77            for(e=0;e<elements->numElements;e++){                  if (elements->Color[e]==color) {
78               if (elements->Color[e]==color) {                     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
79         #else                     add_EM_F=FALSE;
80         {                     add_EM_S=FALSE;
81            for(e=0;e<elements->numElements;e++) {                     /************************************************************* */
82               {                     /* process D */
83         #endif                     /**************************************************************/
84                  Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);                     D_p=getSampleData(D,e);
85                  add_EM_F=FALSE;                     if (NULL!=D_p) {
86                  add_EM_S=FALSE;                       add_EM_S=TRUE;
87                  /************************************************************* */                       if (extendedD) {
88                  /* process D */                           for (s=0;s<p.row_NS;s++) {
89                  /**************************************************************/                             for (r=0;r<p.col_NS;r++) {
90                  D_p=getSampleData(D,e);                               for (k=0;k<p.numEqu;k++) {
91                  if (NULL!=D_p) {                                 for (m=0;m<p.numComp;m++) {
92                    add_EM_S=TRUE;                                   rtmp=0;
93                    if (extendedD) {                                   for (q=0;q<p.numQuad;q++) {
94                        for (s=0;s<p.row_NS;s++) {                                      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)];
95                          for (r=0;r<p.col_NS;r++) {                                   }
96                            for (k=0;k<p.numEqu;k++) {                                   EM_S[INDEX4(k,m,s         ,r         ,p.numEqu,p.numComp,p.row_NN)]= rtmp;
97                              for (m=0;m<p.numComp;m++) {                                   EM_S[INDEX4(k,m,s         ,r+p.col_NS,p.numEqu,p.numComp,p.row_NN)]=-rtmp;
98                                     EM_S[INDEX4(k,m,s+p.row_NS,r         ,p.numEqu,p.numComp,p.row_NN)]=-rtmp;
99                                     EM_S[INDEX4(k,m,s+p.row_NS,r+p.col_NS,p.numEqu,p.numComp,p.row_NN)]= rtmp;
100                                   }
101                                 }
102                               }
103                             }
104                         } else {
105                             for (s=0;s<p.row_NS;s++) {
106                               for (r=0;r<p.col_NS;r++) {
107                                   rtmp=0;
108                                   for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];
109                                   for (k=0;k<p.numEqu;k++) {
110                                       for (m=0;m<p.numComp;m++) {
111                                         rtmp_D=rtmp*D_p[INDEX2(k,m,p.numEqu)];
112                                         EM_S[INDEX4(k,m,s         ,r         ,p.numEqu,p.numComp,p.row_NN)]= rtmp_D;
113                                         EM_S[INDEX4(k,m,s         ,r+p.col_NS,p.numEqu,p.numComp,p.row_NN)]=-rtmp_D;
114                                         EM_S[INDEX4(k,m,s+p.row_NS,r         ,p.numEqu,p.numComp,p.row_NN)]=-rtmp_D;
115                                         EM_S[INDEX4(k,m,s+p.row_NS,r+p.col_NS,p.numEqu,p.numComp,p.row_NN)]= rtmp_D;
116                                      }
117                                   }
118                               }
119                             }
120                         }
121                       }
122                      /**************************************************************/
123                      /*   process Y: */
124                      /**************************************************************/
125                       Y_p=getSampleData(Y,e);
126                       if (NULL!=Y_p) {
127                         add_EM_F=TRUE;
128                         if (extendedY) {
129                            for (s=0;s<p.row_NS;s++) {
130                               for (k=0;k<p.numEqu;k++) {
131                                rtmp=0;                                rtmp=0;
132                                for (q=0;q<p.numQuad;q++) {                                for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[INDEX2(k,q,p.numEqu)];
133                                   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_F[INDEX2(k,s         ,p.numEqu)]=-rtmp;
134                                }                                EM_F[INDEX2(k,s+p.row_NS,p.numEqu)]= rtmp;
135                                EM_S[INDEX4(k,m,s         ,r         ,p.numEqu,p.numComp,p.row_NN)]= rtmp;                             }
                               EM_S[INDEX4(k,m,s         ,r+p.col_NS,p.numEqu,p.numComp,p.row_NN)]=-rtmp;  
                               EM_S[INDEX4(k,m,s+p.row_NS,r         ,p.numEqu,p.numComp,p.row_NN)]=-rtmp;  
                               EM_S[INDEX4(k,m,s+p.row_NS,r+p.col_NS,p.numEqu,p.numComp,p.row_NN)]= rtmp;  
                             }  
                           }  
136                          }                          }
137                        }                        } else {
138                    } else {                          for (s=0;s<p.row_NS;s++) {
                       for (s=0;s<p.row_NS;s++) {  
                         for (r=0;r<p.col_NS;r++) {  
139                              rtmp=0;                              rtmp=0;
140                              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 (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
141                              for (k=0;k<p.numEqu;k++) {                              for (k=0;k<p.numEqu;k++) {
142                                  for (m=0;m<p.numComp;m++) {                                 rtmp_D=rtmp*Y_p[k];
143                                    rtmp_D=rtmp*D_p[INDEX2(k,m,p.numEqu)];                                 EM_F[INDEX2(k,s         ,p.numEqu)]=-rtmp_D;
144                                    EM_S[INDEX4(k,m,s         ,r         ,p.numEqu,p.numComp,p.row_NN)]= rtmp_D;                                 EM_F[INDEX2(k,s+p.row_NS,p.numEqu)]= rtmp_D;
                                   EM_S[INDEX4(k,m,s         ,r+p.col_NS,p.numEqu,p.numComp,p.row_NN)]=-rtmp_D;  
                                   EM_S[INDEX4(k,m,s+p.row_NS,r         ,p.numEqu,p.numComp,p.row_NN)]=-rtmp_D;  
                                   EM_S[INDEX4(k,m,s+p.row_NS,r+p.col_NS,p.numEqu,p.numComp,p.row_NN)]= rtmp_D;  
                                }  
145                              }                              }
146                          }                          }
147                        }                        }
148                    }                      }
149                  }                      /***********************************************************************************************/
150                 /**************************************************************/                      /* add the element matrices onto the matrix and right hand side                                */
151                 /*   process Y: */                      /***********************************************************************************************/
152                 /**************************************************************/                      for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.row_NN)]];
153                  Y_p=getSampleData(Y,e);                      if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);
154                  if (NULL!=Y_p) {                      if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);
155                    add_EM_F=TRUE;    
156                    if (extendedY) {                  } /* end color check */
157                       for (s=0;s<p.row_NS;s++) {               } /* end element loop */
158                          for (k=0;k<p.numEqu;k++) {           } /* end color loop */
159                             rtmp=0;            
160                             for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[INDEX2(k,q,p.numEqu)];           THREAD_MEMFREE(EM_S);
161                             EM_F[INDEX2(k,s         ,p.numEqu)]=-rtmp;           THREAD_MEMFREE(EM_F);
162                             EM_F[INDEX2(k,s+p.row_NS,p.numEqu)]= rtmp;           THREAD_MEMFREE(row_index);
163                          }  
164                       }        } /* end of pointer check */
                    } 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++) {  
                             rtmp_D=rtmp*Y_p[k];  
                             EM_F[INDEX2(k,s         ,p.numEqu)]=-rtmp_D;  
                             EM_F[INDEX2(k,s+p.row_NS,p.numEqu)]= rtmp_D;  
                          }  
                      }  
                    }  
                  }  
                  /***********************************************************************************************/  
                  /* 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.row_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 */  
165     } /* end parallel region */     } /* end parallel region */
166  }  }
167  /*  /*

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

  ViewVC Help
Powered by ViewVC 1.1.26