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

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

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_System2_1D(Ass Line 55  void  Finley_Assemble_PDE_System2_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,k,m;
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_System2_1D(Ass Line 73  void  Finley_Assemble_PDE_System2_1D(Ass
73      dim_t len_EM_F=p.row_NN*p.numEqu;      dim_t len_EM_F=p.row_NN*p.numEqu;
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,k,m,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,k,m;         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
85         #ifndef PASO_MPI            for (color=elements->minColor;color<=elements->maxColor;color++) {
86         for (color=elements->minColor;color<=elements->maxColor;color++) {               /*  open loop over all elements: */
87            /*  open loop over all elements: */               #pragma omp for private(e) schedule(static)
88            #pragma omp for private(e) schedule(static)               for(e=0;e<elements->numElements;e++){
89            for(e=0;e<elements->numElements;e++){                  if (elements->Color[e]==color) {
90               if (elements->Color[e]==color) {            #else
91         #else            {
92         {               for(e=0;e<elements->numElements;e++) {
93            for(e=0;e<elements->numElements;e++) {                  {
94               {            #endif
95         #endif                     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
96                  Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);                     DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);
97                  DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);                     for (q=0;q<len_EM_S;++q) EM_S[q]=0;
98                  for (q=0;q<len_EM_S;++q) EM_S[q]=0;                     for (q=0;q<len_EM_F;++q) EM_F[q]=0;
99                  for (q=0;q<len_EM_F;++q) EM_F[q]=0;                     add_EM_F=FALSE;
100                  add_EM_F=FALSE;                     add_EM_S=FALSE;
101                  add_EM_S=FALSE;                     /**************************************************************/
102                  /**************************************************************/                     /*   process A: */
103                  /*   process A: */                     /**************************************************************/
104                  /**************************************************************/                     A_p=getSampleData(A,e);
105                  A_p=getSampleData(A,e);                     if (NULL!=A_p) {
106                  if (NULL!=A_p) {                        add_EM_S=TRUE;
107                     add_EM_S=TRUE;                        if (extendedA) {
108                     if (extendedA) {                           for (s=0;s<p.row_NS;s++) {
109                        for (s=0;s<p.row_NS;s++) {                             for (r=0;r<p.col_NS;r++) {
110                          for (r=0;r<p.col_NS;r++) {                               for (k=0;k<p.numEqu;k++) {
111                            for (k=0;k<p.numEqu;k++) {                                 for (m=0;m<p.numComp;m++) {
112                              for (m=0;m<p.numComp;m++) {                                   rtmp=0.;
113                                rtmp=0.;                                   for (q=0;q<p.numQuad;q++) {
114                                for (q=0;q<p.numQuad;q++) {                                      rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*
115                                   rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*                                                   A_p[INDEX5(k,0,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
116                                                A_p[INDEX5(k,0,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                                   }
117                                }                                   EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
118                                EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                                 }
119                              }                               }
120                            }                             }
121                          }                           }
122                          } else {
123                             for (s=0;s<p.row_NS;s++) {
124                               for (r=0;r<p.col_NS;r++) {
125                                   rtmp=0;
126                                   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)];
127                                   for (k=0;k<p.numEqu;k++) {
128                                       for (m=0;m<p.numComp;m++) {
129                                            EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)];
130                                      }
131                                   }
132                               }
133                             }
134                           }
135                       }
136                       /**************************************************************/
137                       /*   process B: */
138                       /**************************************************************/
139                       B_p=getSampleData(B,e);
140                       if (NULL!=B_p) {
141                          add_EM_S=TRUE;
142                          if (extendedB) {
143                             for (s=0;s<p.row_NS;s++) {
144                               for (r=0;r<p.col_NS;r++) {
145                                 for (k=0;k<p.numEqu;k++) {
146                                   for (m=0;m<p.numComp;m++) {
147                                     rtmp=0.;
148                                     for (q=0;q<p.numQuad;q++) {
149                                         rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*
150                                                                         B_p[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]*S[INDEX2(r,q,p.row_NS)];
151                                     }
152                                     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
153                                   }
154                                 }
155                               }
156                             }
157                          } else {
158                             for (s=0;s<p.row_NS;s++) {
159                               for (r=0;r<p.col_NS;r++) {
160                                 rtmp=0;
161                                 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)];
162                                 for (k=0;k<p.numEqu;k++) {
163                                     for (m=0;m<p.numComp;m++) {
164                                         EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*B_p[INDEX3(k,0,m,p.numEqu,DIM)];
165                                     }
166                                 }
167                               }
168                             }
169                        }                        }
170                     } else {                     }
171                        for (s=0;s<p.row_NS;s++) {                     /**************************************************************/
172                          for (r=0;r<p.col_NS;r++) {                     /*   process C: */
173                              rtmp=0;                     /**************************************************************/
174                              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)];                     C_p=getSampleData(C,e);
175                              for (k=0;k<p.numEqu;k++) {                     if (NULL!=C_p) {
176                         add_EM_S=TRUE;
177                         if (extendedC) {
178                             for (s=0;s<p.row_NS;s++) {
179                               for (r=0;r<p.col_NS;r++) {
180                                 for (k=0;k<p.numEqu;k++) {
181                                   for (m=0;m<p.numComp;m++) {
182                                     rtmp=0;
183                                     for (q=0;q<p.numQuad;q++) {
184                                          rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*C_p[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
185                                     }
186                                     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
187                                   }
188                                 }
189                               }
190                             }
191                         } else {
192                             for (s=0;s<p.row_NS;s++) {
193                               for (r=0;r<p.col_NS;r++) {
194                                 rtmp=0;
195                                 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)];
196                                 for (k=0;k<p.numEqu;k++) {
197                                  for (m=0;m<p.numComp;m++) {                                  for (m=0;m<p.numComp;m++) {
198                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)];                                     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*C_p[INDEX3(k,m,0,p.numEqu,p.numComp)];
199                                    }
200                                 }
201                               }
202                             }
203                         }
204                       }
205                       /************************************************************* */
206                       /* process D */
207                       /**************************************************************/
208                       D_p=getSampleData(D,e);
209                       if (NULL!=D_p) {
210                         add_EM_S=TRUE;
211                         if (extendedD) {
212                             for (s=0;s<p.row_NS;s++) {
213                               for (r=0;r<p.col_NS;r++) {
214                                 for (k=0;k<p.numEqu;k++) {
215                                   for (m=0;m<p.numComp;m++) {
216                                     rtmp=0;
217                                     for (q=0;q<p.numQuad;q++) {
218                                        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)];
219                                     }
220                                     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
221      
222                                 }                                 }
223                              }                               }
224                          }                             }
225                        }                           }
226                      }                       } else {
227                  }                           for (s=0;s<p.row_NS;s++) {
228                  /**************************************************************/                             for (r=0;r<p.col_NS;r++) {
229                  /*   process B: */                                 rtmp=0;
230                  /**************************************************************/                                 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];
231                  B_p=getSampleData(B,e);                                 for (k=0;k<p.numEqu;k++) {
232                  if (NULL!=B_p) {                                     for (m=0;m<p.numComp;m++) {
233                     add_EM_S=TRUE;                                       EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];
234                     if (extendedB) {                                    }
235                        for (s=0;s<p.row_NS;s++) {                                 }
236                          for (r=0;r<p.col_NS;r++) {                             }
237                            for (k=0;k<p.numEqu;k++) {                           }
238                              for (m=0;m<p.numComp;m++) {                       }
                               rtmp=0.;  
                               for (q=0;q<p.numQuad;q++) {  
                                   rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*  
                                                                   B_p[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]*S[INDEX2(r,q,p.row_NS)];  
                               }  
                               EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;  
                             }  
                           }  
                         }  
                       }  
                    } else {  
                       for (s=0;s<p.row_NS;s++) {  
                         for (r=0;r<p.col_NS;r++) {  
                           rtmp=0;  
                           for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*S[INDEX2(r,q,p.row_NS)];  
                           for (k=0;k<p.numEqu;k++) {  
                               for (m=0;m<p.numComp;m++) {  
                                   EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*B_p[INDEX3(k,0,m,p.numEqu,DIM)];  
                               }  
                           }  
                         }  
                       }  
239                     }                     }
240                  }                     /**************************************************************/
241                  /**************************************************************/                     /*   process X: */
242                  /*   process C: */                     /**************************************************************/
243                  /**************************************************************/                     X_p=getSampleData(X,e);
244                  C_p=getSampleData(C,e);                     if (NULL!=X_p) {
245                  if (NULL!=C_p) {                       add_EM_F=TRUE;
246                    add_EM_S=TRUE;                       if (extendedX) {
247                    if (extendedC) {                          for (s=0;s<p.row_NS;s++) {
248                        for (s=0;s<p.row_NS;s++) {                             for (k=0;k<p.numEqu;k++) {
                         for (r=0;r<p.col_NS;r++) {  
                           for (k=0;k<p.numEqu;k++) {  
                             for (m=0;m<p.numComp;m++) {  
249                                rtmp=0;                                rtmp=0;
250                                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[INDEX3(k,0,q,p.numEqu,DIM)];
251                                     rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*C_p[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                                EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
252                                }                             }
                               EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;  
                             }  
                           }  
253                          }                          }
254                        }                       } else {
255                    } else {                          for (s=0;s<p.row_NS;s++) {
256                        for (s=0;s<p.row_NS;s++) {                             rtmp=0;
257                          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)];
258                            rtmp=0;                             for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*X_p[INDEX2(k,0,p.numEqu)];
                           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)];  
                           for (k=0;k<p.numEqu;k++) {  
                              for (m=0;m<p.numComp;m++) {  
                                 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*C_p[INDEX3(k,m,0,p.numEqu,p.numComp)];  
                              }  
                           }  
259                          }                          }
260                        }                       }
261                    }                    }
262                  }                    /**************************************************************/
263                  /************************************************************* */                    /*   process Y: */
264                  /* process D */                    /**************************************************************/
265                  /**************************************************************/                     Y_p=getSampleData(Y,e);
266                  D_p=getSampleData(D,e);                     if (NULL!=Y_p) {
267                  if (NULL!=D_p) {                       add_EM_F=TRUE;
268                    add_EM_S=TRUE;                       if (extendedY) {
269                    if (extendedD) {                          for (s=0;s<p.row_NS;s++) {
270                        for (s=0;s<p.row_NS;s++) {                             for (k=0;k<p.numEqu;k++) {
                         for (r=0;r<p.col_NS;r++) {  
                           for (k=0;k<p.numEqu;k++) {  
                             for (m=0;m<p.numComp;m++) {  
271                                rtmp=0;                                rtmp=0;
272                                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)];
273                                   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;
274                                }                             }
                               EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;  
   
                             }  
                           }  
275                          }                          }
276                        }                        } else {
277                    } 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++) {  
278                              rtmp=0;                              rtmp=0;
279                              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)];
280                              for (k=0;k<p.numEqu;k++) {                              for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];
                                 for (m=0;m<p.numComp;m++) {  
                                   EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];  
                                }  
                             }  
281                          }                          }
282                        }                        }
283                    }                      }
284                  }                      /***********************************************************************************************/
285                  /**************************************************************/                      /* add the element matrices onto the matrix and right hand side                                */
286                  /*   process X: */                      /***********************************************************************************************/
287                  /**************************************************************/                      for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
288                  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);
289                  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);
290                    add_EM_F=TRUE;    
291                    if (extendedX) {                  } /* end color check */
292                       for (s=0;s<p.row_NS;s++) {               } /* end element loop */
293                          for (k=0;k<p.numEqu;k++) {           } /* end color loop */
294                             rtmp=0;            
295                             for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX3(k,0,q,p.numEqu,DIM)];           THREAD_MEMFREE(EM_S);
296                             EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;           THREAD_MEMFREE(EM_F);
297                          }           THREAD_MEMFREE(row_index);
298                       }  
299                    } else {        } /* end of pointer check */
                      for (s=0;s<p.row_NS;s++) {  
                         rtmp=0;  
                         for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];  
                         for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*X_p[INDEX2(k,0,p.numEqu)];  
                      }  
                   }  
                }  
                /**************************************************************/  
                /*   process Y: */  
                /**************************************************************/  
                 Y_p=getSampleData(Y,e);  
                 if (NULL!=Y_p) {  
                   add_EM_F=TRUE;  
                   if (extendedY) {  
                      for (s=0;s<p.row_NS;s++) {  
                         for (k=0;k<p.numEqu;k++) {  
                            rtmp=0;  
                            for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[INDEX2(k,q,p.numEqu)];  
                            EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;  
                         }  
                      }  
                    } else {  
                      for (s=0;s<p.row_NS;s++) {  
                          rtmp=0;  
                          for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];  
                          for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];  
                      }  
                    }  
                  }  
                  /***********************************************************************************************/  
                  /* add the element matrices onto the matrix and right hand side                                */  
                  /***********************************************************************************************/  
                  for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];  
                  if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);  
                  if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);  
   
              } /* end color check */  
           } /* end element loop */  
       } /* end color loop */  
300     } /* end parallel region */     } /* end parallel region */
301  }  }
302  /*  /*

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

  ViewVC Help
Powered by ViewVC 1.1.26