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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26