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

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

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

trunk/finley/src/Assemble_PDE_Single2_1D.c revision 798 by gross, Fri Aug 4 01:05:36 2006 UTC temp/finley/src/Assemble_PDE_Single2_1D.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 1.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 33  Line 36 
36    
37  /**************************************************************/  /**************************************************************/
38    
 /*  Author: gross@access.edu.au */  
 /*  Version: $Id:$ */  
   
 /**************************************************************/  
   
39    
40  #include "Assemble.h"  #include "Assemble.h"
41  #include "Util.h"  #include "Util.h"
42    #ifdef _OPENMP
43    #include <omp.h>
44    #endif
45    
46    
47  /**************************************************************/  /**************************************************************/
48    
# Line 51  void  Finley_Assemble_PDE_Single2_1D(Ass Line 53  void  Finley_Assemble_PDE_Single2_1D(Ass
53      #define DIM 1      #define DIM 1
54      index_t color;      index_t color;
55      dim_t e;      dim_t e;
56        double *EM_S, *EM_F, *Vol, *DSDX, *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p;
57        index_t *row_index;
58        register dim_t q, s,r;
59        register double rtmp;
60        bool_t add_EM_F, add_EM_S;
61    
62      bool_t extendedA=isExpanded(A);      bool_t extendedA=isExpanded(A);
63      bool_t extendedB=isExpanded(B);      bool_t extendedB=isExpanded(B);
64      bool_t extendedC=isExpanded(C);      bool_t extendedC=isExpanded(C);
# Line 63  void  Finley_Assemble_PDE_Single2_1D(Ass Line 71  void  Finley_Assemble_PDE_Single2_1D(Ass
71      dim_t len_EM_F=p.row_NN;      dim_t len_EM_F=p.row_NN;
72    
73    
74      #pragma omp parallel private(color)      #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,rtmp,add_EM_F, add_EM_S)
75      {      {
76         double EM_S[len_EM_S], EM_F[len_EM_F];         EM_S=THREAD_MEMALLOC(len_EM_S,double);
77         index_t row_index[p.row_NN];         EM_F=THREAD_MEMALLOC(len_EM_F,double);
78         register dim_t q, s,r;         row_index=THREAD_MEMALLOC(p.row_NN,index_t);
79         register double rtmp;  
80         double *Vol, *DSDX, *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p;         if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) {
81         bool_t add_EM_F, add_EM_S;  
82         #ifndef PASO_MPI            for (color=elements->minColor;color<=elements->maxColor;color++) {
83         for (color=elements->minColor;color<=elements->maxColor;color++) {               /*  open loop over all elements: */
84            /*  open loop over all elements: */               #pragma omp for private(e) schedule(static)
85            #pragma omp for private(e) schedule(static)               for(e=0;e<elements->numElements;e++){
86            for(e=0;e<elements->numElements;e++){                  if (elements->Color[e]==color) {
87               if (elements->Color[e]==color) {                     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
88         #else                     DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);
89         {                     for (q=0;q<len_EM_S;++q) EM_S[q]=0;
90            for(e=0;e<elements->numElements;e++) {                     for (q=0;q<len_EM_F;++q) EM_F[q]=0;
91               {                     add_EM_F=FALSE;
92         #endif                     add_EM_S=FALSE;
93                  Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);                     /**************************************************************/
94                  DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);                     /*   process A: */
95                  for (q=0;q<len_EM_S;++q) EM_S[q]=0;                     /**************************************************************/
96                  for (q=0;q<len_EM_F;++q) EM_F[q]=0;                     A_p=getSampleData(A,e);
97                  add_EM_F=FALSE;                     if (NULL!=A_p) {
98                  add_EM_S=FALSE;                        add_EM_S=TRUE;
99                  /**************************************************************/                        if (extendedA) {
100                  /*   process A: */                           for (s=0;s<p.row_NS;s++) {
101                  /**************************************************************/                             for (r=0;r<p.col_NS;r++) {
102                  A_p=getSampleData(A,e);                                rtmp=0;
103                  if (NULL!=A_p) {                                for (q=0;q<p.numQuad;q++) {
104                     add_EM_S=TRUE;                                   rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX3(0,0,q,DIM,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
105                     if (extendedA) {                               }
106                        for (s=0;s<p.row_NS;s++) {                               EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
107                          for (r=0;r<p.col_NS;r++) {                             }
108                             rtmp=0;                           }
109                             for (q=0;q<p.numQuad;q++) {                        } else {
110                                rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX3(0,0,q,DIM,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                           for (s=0;s<p.row_NS;s++) {
111                            }                             for (r=0;r<p.col_NS;r++) {
112                            EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                                 rtmp=0;
113                          }                                 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
114                        }                                 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*A_p[INDEX2(0,0,DIM)];
115                     } else {                             }
116                        for (s=0;s<p.row_NS;s++) {                           }
117                          for (r=0;r<p.col_NS;r++) {                         }
118                              rtmp=0;                     }
119                              for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                     /**************************************************************/
120                              EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*A_p[INDEX2(0,0,DIM)];                     /*   process B: */
121                          }                     /**************************************************************/
122                        }                     B_p=getSampleData(B,e);
123                      }                     if (NULL!=B_p) {
124                  }                        add_EM_S=TRUE;
125                  /**************************************************************/                        if (extendedB) {
126                  /*   process B: */                           for (s=0;s<p.row_NS;s++) {
127                  /**************************************************************/                             for (r=0;r<p.col_NS;r++) {
128                  B_p=getSampleData(B,e);                               rtmp=0;
129                  if (NULL!=B_p) {                               for (q=0;q<p.numQuad;q++) {
130                     add_EM_S=TRUE;                                  rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*B_p[INDEX2(0,q,DIM)]*S[INDEX2(r,q,p.row_NS)];
131                     if (extendedB) {                               }
132                        for (s=0;s<p.row_NS;s++) {                               EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
133                          for (r=0;r<p.col_NS;r++) {                             }
134                            rtmp=0;                           }
135                            for (q=0;q<p.numQuad;q++) {                        } else {
136                               rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*B_p[INDEX2(0,q,DIM)]*S[INDEX2(r,q,p.row_NS)];                           for (s=0;s<p.row_NS;s++) {
137                            }                             for (r=0;r<p.col_NS;r++) {
138                            EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                                 rtmp=0;
139                          }                                 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*S[INDEX2(r,q,p.row_NS)];
140                        }                                 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*B_p[0];
141                     } else {                             }
142                        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]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*S[INDEX2(r,q,p.row_NS)];  
                             EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*B_p[0];  
                         }  
143                        }                        }
144                     }                     }
145                  }                     /**************************************************************/
146                  /**************************************************************/                     /*   process C: */
147                  /*   process C: */                     /**************************************************************/
148                  /**************************************************************/                     C_p=getSampleData(C,e);
149                  C_p=getSampleData(C,e);                     if (NULL!=C_p) {
150                  if (NULL!=C_p) {                        add_EM_S=TRUE;
151                     add_EM_S=TRUE;                       if (extendedC) {
152                    if (extendedC) {                           for (s=0;s<p.row_NS;s++) {
153                        for (s=0;s<p.row_NS;s++) {                             for (r=0;r<p.col_NS;r++) {
154                          for (r=0;r<p.col_NS;r++) {                               rtmp=0;
155                                 for (q=0;q<p.numQuad;q++) {
156                                    rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*C_p[INDEX2(0,q,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
157                                 }
158                                 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
159                               }
160                             }
161                         } else {
162                             for (s=0;s<p.row_NS;s++) {
163                               for (r=0;r<p.col_NS;r++) {
164                                  rtmp=0;
165                                  for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
166                                  EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*C_p[0];
167                               }
168                             }
169                         }
170                       }
171                       /************************************************************* */
172                       /* process D */
173                       /**************************************************************/
174                       D_p=getSampleData(D,e);
175                       if (NULL!=D_p) {
176                         add_EM_S=TRUE;
177                         if (extendedD) {
178                             for (s=0;s<p.row_NS;s++) {
179                               for (r=0;r<p.col_NS;r++) {
180                                  rtmp=0;
181                                  for (q=0;q<p.numQuad;q++) {
182                                     rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q]*S[INDEX2(r,q,p.row_NS)];
183                                 }
184                                 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
185                               }
186                             }
187                         } else {
188                             for (s=0;s<p.row_NS;s++) {
189                               for (r=0;r<p.col_NS;r++) {
190                                   rtmp=0;
191                                   for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];
192                                   EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[0];
193                               }
194                             }
195                         }
196                       }
197                       /**************************************************************/
198                       /*   process X: */
199                       /**************************************************************/
200                       X_p=getSampleData(X,e);
201                       if (NULL!=X_p) {
202                         add_EM_F=TRUE;
203                         if (extendedX) {
204                            for (s=0;s<p.row_NS;s++) {
205                            rtmp=0;                            rtmp=0;
206                            for (q=0;q<p.numQuad;q++) {                            for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX2(0,q,DIM)];
207                               rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*C_p[INDEX2(0,q,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                            EM_F[INDEX2(0,s,p.numEqu)]+=rtmp;
                           }  
                           EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;  
208                          }                          }
209                        }                       } else {
210                    } else {                          for (s=0;s<p.row_NS;s++) {
211                        for (s=0;s<p.row_NS;s++) {                            rtmp=0;
212                          for (r=0;r<p.col_NS;r++) {                            for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
213                             rtmp=0;                            EM_F[INDEX2(0,s,p.numEqu)]+=rtmp*X_p[0];
                            for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];  
                            EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*C_p[0];  
214                          }                          }
215                        }                       }
216                    }                    }
217                  }                    /**************************************************************/
218                  /************************************************************* */                    /*   process Y: */
219                  /* process D */                    /**************************************************************/
220                  /**************************************************************/                     Y_p=getSampleData(Y,e);
221                  D_p=getSampleData(D,e);                     if (NULL!=Y_p) {
222                  if (NULL!=D_p) {                       add_EM_F=TRUE;
223                    add_EM_S=TRUE;                       if (extendedY) {
224                    if (extendedD) {                          for (s=0;s<p.row_NS;s++) {
                       for (s=0;s<p.row_NS;s++) {  
                         for (r=0;r<p.col_NS;r++) {  
225                             rtmp=0;                             rtmp=0;
226                             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[q];
227                                rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q]*S[INDEX2(r,q,p.row_NS)];                             EM_F[INDEX2(0,s,p.numEqu)]+=rtmp;
                           }  
                           EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;  
228                          }                          }
229                        }                        } else {
230                    } 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++) {  
231                              rtmp=0;                              rtmp=0;
232                              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)];
233                              EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[0];                              EM_F[INDEX2(0,s,p.numEqu)]+=rtmp*Y_p[0];
234                          }                          }
235                        }                        }
236                    }                      }
237                  }                      /***********************************************************************************************/
238                  /**************************************************************/                      /* add the element matrices onto the matrix and right hand side                                */
239                  /*   process X: */                      /***********************************************************************************************/
240                  /**************************************************************/                      for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
241                  X_p=getSampleData(X,e);                      if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);
242                  if (NULL!=X_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);
243                    add_EM_F=TRUE;    
244                    if (extendedX) {                  } /* end color check */
245                       for (s=0;s<p.row_NS;s++) {               } /* end element loop */
246                         rtmp=0;           } /* end color loop */
247                         for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX2(0,q,DIM)];            
248                         EM_F[INDEX2(0,s,p.numEqu)]+=rtmp;           THREAD_MEMFREE(EM_S);
249                       }           THREAD_MEMFREE(EM_F);
250                    } else {           THREAD_MEMFREE(row_index);
251                       for (s=0;s<p.row_NS;s++) {  
252                         rtmp=0;        } /* end of pointer check */
                        for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];  
                        EM_F[INDEX2(0,s,p.numEqu)]+=rtmp*X_p[0];  
                      }  
                   }  
                }  
                /**************************************************************/  
                /*   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++) {  
                         rtmp=0;  
                         for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[q];  
                         EM_F[INDEX2(0,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)];  
                          EM_F[INDEX2(0,s,p.numEqu)]+=rtmp*Y_p[0];  
                      }  
                    }  
                  }  
                  /***********************************************************************************************/  
                  /* 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 */  
253     } /* end parallel region */     } /* end parallel region */
254  }  }
255  /*  /*

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

  ViewVC Help
Powered by ViewVC 1.1.26