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

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

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

revision 2881 by jfenwick, Thu Jan 28 02:03:15 2010 UTC revision 2988 by gross, Tue Mar 16 06:14:30 2010 UTC
# Line 30  Line 30 
30    
31    
32  /* Disabled until the tests pass */  /* Disabled until the tests pass */
33  /* #define NEW_LUMPING */ /* */  #define NEW_LUMPING
34    
35  /**************************************************************/  /**************************************************************/
36    
# Line 122  void Finley_Assemble_LumpedSystem(Finley Line 122  void Finley_Assemble_LumpedSystem(Finley
122         EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);         EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
123         row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);         row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);
124         if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {         if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
125            if (p.numEqu == 1) {            if (p.numEqu == 1) { /* single equation */
126               if (expandedD) {               if (expandedD) {   /* with expanded D */        
               
127                   for (color=elements->minColor;color<=elements->maxColor;color++) {                   for (color=elements->minColor;color<=elements->maxColor;color++) {
128                      /*  open loop over all elements: */                      /*  open loop over all elements: */
129                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
130                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){              
131                                    if (elements->Color[e]==color) {
                         if (elements->Color[e]==color) {  
132                              for (isub=0; isub<p.numSub; isub++) {                              for (isub=0; isub<p.numSub; isub++) {
133                                 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);                                 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);;
                                memset(EM_lumpedMat,0,len_EM_lumpedMat_size);  
134                                    
135                                 D_p=getSampleDataRO(D,e);                                                           D_p=getSampleDataRO(D,e);                          
136                                 #ifdef NEW_LUMPING /* HRZ lumping */                     #ifdef NEW_LUMPING /* HRZ lumping */
137    
138                                     m_t=0; /* mass of the element: m_t */                                     m_t=0; /* mass of the element: m_t */
139                                     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub) ];                                     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub) ];
140                                                      
141                                     diagS=0; /* diagonal sum: S */                                     diagS=0; /* diagonal sum: S */
142                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
143                                        rtmp=0;                                        rtmp=0;
# Line 148  void Finley_Assemble_LumpedSystem(Finley Line 146  void Finley_Assemble_LumpedSystem(Finley
146                                        diagS+=rtmp;                                        diagS+=rtmp;
147                                     }                                     }
148                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
149                                     rtmp=m_t/diagS;                     rtmp=m_t/diagS;
150                                     for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;                                     for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
151                                                        
152                                 #else /* row-sum lumping */                                 #else /* row-sum lumping */
# Line 161  void Finley_Assemble_LumpedSystem(Finley Line 159  void Finley_Assemble_LumpedSystem(Finley
159                                 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];                                 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
160                                 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                                 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
161                  } /* end of isub loop */                  } /* end of isub loop */
162                         } /* end color check */                         } /* end color check */        
                 
163                      } /* end element loop */                      } /* end element loop */
164                    } /* end color loop */                    } /* end color loop */
165               } else  {                         } else  {  /* with constant D */  
166                   for (color=elements->minColor;color<=elements->maxColor;color++) {                   for (color=elements->minColor;color<=elements->maxColor;color++) {
167                      /*  open loop over all elements: */                      /*  open loop over all elements: */
168                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
169                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){
170                                                if (elements->Color[e]==color) {
                         if (elements->Color[e]==color) {  
171                              for (isub=0; isub<p.numSub; isub++) {                              for (isub=0; isub<p.numSub; isub++) {
172                                                               Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);                
                                 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);  
                                 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);  
                   
173                                 D_p=getSampleDataRO(D,e);                                                           D_p=getSampleDataRO(D,e);                          
174                                 #ifdef NEW_LUMPING /* HRZ lumping */                     #ifdef NEW_LUMPING /* HRZ lumping */
175                                     m_t=0; /* mass of the element: m_t */                                     m_t=0; /* mass of the element: m_t */
176                                     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];                                     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
                                    m_t*=D_p[0];  
177                                        
178                                     diagS=0; /* diagonal sum: S */                                     diagS=0; /* diagonal sum: S */
179                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
180                                        rtmp=0;                                        rtmp=0;
181                                        for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                                        for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
                                       rtmp*=D_p[0];  
182                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
183                                        diagS+=rtmp;                                        diagS+=rtmp;
184                                     }                                     }
185                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
186                                     rtmp=m_t/diagS;                     rtmp=m_t/diagS*D_p[0];
187                                     for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;                                     for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
188                                                        
189                                 #else /* row-sum lumping */                                 #else /* row-sum lumping */
190                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
191                                         rtmp=0;                                         rtmp=0;
192                                         for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];                                         for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
193                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];
194                                     }                                     }
195                                 #endif                                 #endif
196                                 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];                                 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
# Line 210  void Finley_Assemble_LumpedSystem(Finley Line 201  void Finley_Assemble_LumpedSystem(Finley
201                    } /* end color loop */                    } /* end color loop */
202                    
203               }               }
204            } else {            } else { /* system of  equation */
205               if (expandedD) {               if (expandedD) { /* with expanded D */
206                   for (color=elements->minColor;color<=elements->maxColor;color++) {                   for (color=elements->minColor;color<=elements->maxColor;color++) {
207                      /*  open loop over all elements: */                      /*  open loop over all elements: */
208                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
209                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){
210                         if (elements->Color[e]==color) {                         if (elements->Color[e]==color) {
211                             for (isub=0; isub<p.numSub; isub++) {                    for (isub=0; isub<p.numSub; isub++) {    
212                                Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);                                Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
213                                memset(EM_lumpedMat,0,len_EM_lumpedMat_size);                                D_p=getSampleDataRO(D,e);  
214                                D_p=getSampleDataRO(D,e);                  
                 
215                                #ifdef NEW_LUMPING /* HRZ lumping */                                #ifdef NEW_LUMPING /* HRZ lumping */
216                                    for (k=0;k<p.numEqu;k++) {                                    for (k=0;k<p.numEqu;k++) {
217                                        m_t=0; /* mass of the element: m_t */                                        m_t=0; /* mass of the element: m_t */
# Line 233  void Finley_Assemble_LumpedSystem(Finley Line 223  void Finley_Assemble_LumpedSystem(Finley
223                                            for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                                            for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
224                                            EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;                                            EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
225                                            diagS+=rtmp;                                            diagS+=rtmp;
226                                        }                                         }
227                                        /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                         /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
228                                        rtmp=m_t/diagS;                         rtmp=m_t/diagS;
229                                        for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;                         for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;
230                                    }                      }                
231                                #else /* row-sum lumping */                  #else /* row-sum lumping */
232                                  for (s=0;s<p.row_numShapes;s++) {                                    for (s=0;s<p.row_numShapes;s++) {
233                                    for (k=0;k<p.numEqu;k++) {                                        for (k=0;k<p.numEqu;k++) {
234                                        rtmp=0.;                                           rtmp=0.;
235                                        for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)];                                            for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)];
236                                        EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;                                             EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
237                                    }                                         }
238                                  }                    }
239                                #endif                    #endif
240                                for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];                    for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
241                                Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                    Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
242                             } /* end of isub loop */                 } /* end of isub loop */
243                         } /* end color check */                         } /* end color check */
244                      } /* end element loop */                      } /* end element loop */
245                  } /* end color loop */                  } /* end color loop */
246               } else {               } else { /* with constant D */
247                   for (color=elements->minColor;color<=elements->maxColor;color++) {                   for (color=elements->minColor;color<=elements->maxColor;color++) {
248                      /*  open loop over all elements: */                      /*  open loop over all elements: */
249                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
250                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){
251                         if (elements->Color[e]==color) {                         if (elements->Color[e]==color) {
252                             for (isub=0; isub<p.numSub; isub++) {                     for (isub=0; isub<p.numSub; isub++) {    
253                                Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);                                Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
254                                memset(EM_lumpedMat,0,len_EM_lumpedMat_size);                                memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
255                                D_p=getSampleDataRO(D,e);                                D_p=getSampleDataRO(D,e);
256                                
257                                #ifdef NEW_LUMPING /* HRZ lumping */                                #ifdef NEW_LUMPING /* HRZ lumping */
                                   for (k=0;k<p.numEqu;k++) {  
258                                        m_t=0; /* mass of the element: m_t */                                        m_t=0; /* mass of the element: m_t */
259                                        for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)];                                        for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
                                       m_t*=D_p[k];            
260                                        diagS=0; /* diagonal sum: S */                                        diagS=0; /* diagonal sum: S */
261                                        for (s=0;s<p.row_numShapes;s++) {                                        for (s=0;s<p.row_numShapes;s++) {
262                                            rtmp=0;                                            rtmp=0;
263                                            for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                                            for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
264                                            rtmp*=D_p[k];                        for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
                                           EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;  
265                                            diagS+=rtmp;                                            diagS+=rtmp;
266                                        }                                        }
267                                        /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                        /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
268                                        rtmp=m_t/diagS;                        rtmp=m_t/diagS;
269                                        for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;                        for (s=0;s<p.row_numShapes;s++) {
270                                    }                          for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp*D_p[k];
271                               #else /* row-sum lumping */                        }
272                       #else /* row-sum lumping */
273                                   for (s=0;s<p.row_numShapes;s++) {                                   for (s=0;s<p.row_numShapes;s++) {
274                                       for (k=0;k<p.numEqu;k++) {                      for (k=0;k<p.numEqu;k++) {
275                                          rtmp=0.;                                          rtmp=0.;
276                                          for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];                                          for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
277                                          EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];                                          EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
278                                       }                                       }
279                                   }                   }
280                                #endif                    #endif
281                                for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];                    for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
282                                Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                    Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
283                            } /* end of isub loop */                } /* end of isub loop */
284                         } /* end color check */                         } /* end color check */
285                      } /* end element loop */                      } /* end element loop */
286                  } /* end color loop */                  } /* end color loop */

Legend:
Removed from v.2881  
changed lines
  Added in v.2988

  ViewVC Help
Powered by ViewVC 1.1.26