/[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 3378 by jfenwick, Mon Oct 11 01:48:14 2010 UTC revision 3379 by gross, Wed Nov 24 04:48:49 2010 UTC
# Line 29  Line 29 
29  #endif  #endif
30    
31    
 /* Disabled until the tests pass */  
 /* #define NEW_LUMPING */  
   
32  /**************************************************************/  /**************************************************************/
33    
34  void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)  void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D, const bool_t useHRZ)
35  {  {
36    
37    bool_t reducedIntegrationOrder=FALSE, expandedD;    bool_t reducedIntegrationOrder=FALSE, expandedD;
# Line 47  void Finley_Assemble_LumpedSystem(Finley Line 44  void Finley_Assemble_LumpedSystem(Finley
44    double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *lumpedMat_p=NULL;    double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *lumpedMat_p=NULL;
45    register double rtmp;    register double rtmp;
46    size_t len_EM_lumpedMat_size;    size_t len_EM_lumpedMat_size;
47  #ifdef NEW_LUMPING    register double m_t=0., diagS=0., rtmp2=0.;
   register double m_t=0., diagS=0.;  
 #endif  
48    
49    Finley_resetError();    Finley_resetError();
50    
# Line 110  void Finley_Assemble_LumpedSystem(Finley Line 105  void Finley_Assemble_LumpedSystem(Finley
105      expandedD=isExpanded(D);      expandedD=isExpanded(D);
106      S=p.row_jac->BasisFunctions->S;      S=p.row_jac->BasisFunctions->S;
107    
108  #ifdef NEW_LUMPING      #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t, isub, rtmp2)
     #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t, isub)  
 #else  
     #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, isub)  
 #endif  
109      {      {
110         EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);         EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
111         row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);         row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);
# Line 129  void Finley_Assemble_LumpedSystem(Finley Line 120  void Finley_Assemble_LumpedSystem(Finley
120                              for (isub=0; isub<p.numSub; isub++) {                              for (isub=0; isub<p.numSub; isub++) {
121                                 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)]);
122                                 D_p=getSampleDataRO(D,e);                                                           D_p=getSampleDataRO(D,e);                          
123                     #ifdef NEW_LUMPING /* HRZ lumping */                     if (useHRZ) {
124    
125                                     m_t=0; /* mass of the element: m_t */                                     m_t=0; /* mass of the element: m_t */
126                                       #pragma ivdep
127                                     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) ];
128                                                        
129                                     diagS=0; /* diagonal sum: S */                                     diagS=0; /* diagonal sum: S */
130                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
131                                        rtmp=0;                                        rtmp=0;
132                                        for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                        #pragma ivdep
133                                          for (q=0;q<p.numQuadSub;q++) {
134                         rtmp2=S[INDEX2(s,q,p.row_numShapes)];
135                         rtmp+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub)]*rtmp2*rtmp2;
136                          }
137                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
138                                        diagS+=rtmp;                                        diagS+=rtmp;
139                                     }                                     }
140                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
141                     rtmp=m_t/diagS;                     rtmp=m_t/diagS;
142                       #pragma ivdep
143                                     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;
144                                                        
145                                 #else /* row-sum lumping */                     } else { /* row-sum lumping */
146                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
147                                         rtmp=0;                                         rtmp=0;
148                           #pragma ivdep
149                                         for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX2(q, isub,p.numQuadSub)];                                         for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX2(q, isub,p.numQuadSub)];
150                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
151                                     }                                     }
152    
153                                 #endif                     }
154                                 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)]];
155                                 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);
156                  } /* end of isub loop */                  } /* end of isub loop */
157                         } /* end color check */                                 } /* end color check */  
158                      } /* end element loop */                      } /* end element loop */
159                    } /* end color loop */                    } /* end color loop */
160    
161    
162               } else  {  /* with constant D */                 } else  {  /* with constant D */  
163    
164                   for (color=elements->minColor;color<=elements->maxColor;color++) {                   for (color=elements->minColor;color<=elements->maxColor;color++) {
# Line 169  void Finley_Assemble_LumpedSystem(Finley Line 169  void Finley_Assemble_LumpedSystem(Finley
169                              for (isub=0; isub<p.numSub; isub++) {                              for (isub=0; isub<p.numSub; isub++) {
170                                 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)]);                
171                                 D_p=getSampleDataRO(D,e);                                                           D_p=getSampleDataRO(D,e);                          
172                     #ifdef NEW_LUMPING /* HRZ lumping */                     if (useHRZ) { /* HRZ lumping */
173                                     m_t=0; /* mass of the element: m_t */                                     m_t=0; /* mass of the element: m_t */
174                                       #pragma ivdep
175                                     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];                                     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
176                                     diagS=0; /* diagonal sum: S */                                     diagS=0; /* diagonal sum: S */
177                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
178                                        rtmp=0;                                        rtmp=0;
179                                        for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                        #pragma ivdep
180                                          for (q=0;q<p.numQuadSub;q++){
181                         rtmp2=S[INDEX2(s,q,p.row_numShapes)];
182                          rtmp+=Vol[q]*rtmp2*rtmp2;
183                          }
184                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
185                                        diagS+=rtmp;                                        diagS+=rtmp;
186                                     }                                     }
187                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
188                     rtmp=m_t/diagS*D_p[0];                     rtmp=m_t/diagS*D_p[0];
189                       #pragma ivdep
190                                     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;
191                                 #else /* row-sum lumping */                     } else { /* row-sum lumping */
192                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
193                                         rtmp=0;                                         rtmp=0;
194                           #pragma ivdep
195                                         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)];
196                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];
197                                     }                                     }
198                                 #endif                     }
199                                 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)]];
200                                 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);
201                              } /* end of isub loop */                              } /* end of isub loop */
# Line 208  void Finley_Assemble_LumpedSystem(Finley Line 215  void Finley_Assemble_LumpedSystem(Finley
215                                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)]);
216                                D_p=getSampleDataRO(D,e);                                  D_p=getSampleDataRO(D,e);  
217                                    
218                                #ifdef NEW_LUMPING /* HRZ lumping */                                if (useHRZ) { /* HRZ lumping */
219                                    for (k=0;k<p.numEqu;k++) {                                    for (k=0;k<p.numEqu;k++) {
220                                        m_t=0; /* mass of the element: m_t */                                        m_t=0; /* mass of the element: m_t */
221                                          #pragma ivdep
222                                        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]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)];
223                                                                                                        
224                                        diagS=0; /* diagonal sum: S */                                        diagS=0; /* diagonal sum: S */
225                                        for (s=0;s<p.row_numShapes;s++) {                                        for (s=0;s<p.row_numShapes;s++) {
226                                            rtmp=0;                                            rtmp=0;
227                                            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)];                        #pragma ivdep
228                                              for (q=0;q<p.numQuadSub;q++) {
229                             rtmp2=S[INDEX2(s,q,p.row_numShapes)];
230                             rtmp+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)]*rtmp2*rtmp2;
231                          }
232                                            EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;                                            EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
233                                            diagS+=rtmp;                                            diagS+=rtmp;
234                                         }                                         }
235                         /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                         /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
236                         rtmp=m_t/diagS;                         rtmp=m_t/diagS;
237                           #pragma ivdep
238                         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;
239                      }                                      }                
240                  #else /* row-sum lumping */                    } else { /* row-sum lumping */
241                                    for (s=0;s<p.row_numShapes;s++) {                                    for (s=0;s<p.row_numShapes;s++) {
242                                        for (k=0;k<p.numEqu;k++) {                                        for (k=0;k<p.numEqu;k++) {
243                                           rtmp=0.;                                           rtmp=0.;
244                                            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)];                       #pragma ivdep
245                                             EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;                                           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)];
246                                             EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
247                                         }                                         }
248                    }                    }
249                    #endif                    }
250                    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)]];
251                    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);
252                 } /* end of isub loop */                 } /* end of isub loop */
253                         } /* end color check */                         } /* end color check */
254                      } /* end element loop */                      } /* end element loop */
# Line 249  void Finley_Assemble_LumpedSystem(Finley Line 263  void Finley_Assemble_LumpedSystem(Finley
263                                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)]);
264                                D_p=getSampleDataRO(D,e);                                D_p=getSampleDataRO(D,e);
265                                
266                                #ifdef NEW_LUMPING /* HRZ lumping */                                if (useHRZ) { /* HRZ lumping */
267                                        m_t=0; /* mass of the element: m_t */                                        m_t=0; /* mass of the element: m_t */
268                                          #pragma ivdep
269                                        for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];                                        for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
270                                        diagS=0; /* diagonal sum: S */                                        diagS=0; /* diagonal sum: S */
271                                        for (s=0;s<p.row_numShapes;s++) {                                        for (s=0;s<p.row_numShapes;s++) {
272                                            rtmp=0;                                            rtmp=0;
273                                            for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                        #pragma ivdep
274                                              for (q=0;q<p.numQuadSub;q++) {
275                             rtmp2=S[INDEX2(s,q,p.row_numShapes)];
276                             rtmp+=Vol[q]*rtmp2*rtmp2;
277                          }
278                          #pragma ivdep
279                        for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;                        for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
280                                            diagS+=rtmp;                                            diagS+=rtmp;
281                                        }                                        }
282                                          
283                        /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                        /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
284                        rtmp=m_t/diagS;                        rtmp=m_t/diagS;
285                        for (s=0;s<p.row_numShapes;s++) {                        for (s=0;s<p.row_numShapes;s++) {
286                            #pragma ivdep
287                          for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp*D_p[k];                          for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp*D_p[k];
288                        }                        }
289                     #else /* row-sum lumping */                    } else { /* row-sum lumping */
290                                   for (s=0;s<p.row_numShapes;s++) {                                   for (s=0;s<p.row_numShapes;s++) {
291                      for (k=0;k<p.numEqu;k++) {                      for (k=0;k<p.numEqu;k++) {
292                                          rtmp=0.;                                          rtmp=0.;
293                        #pragma ivdep
294                                          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)];
295                                          EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];                                          EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];
296                                       }                                       }
297                   }                   }
298                    #endif                    }
299                    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)]];
300                    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);
301                } /* end of isub loop */                } /* end of isub loop */

Legend:
Removed from v.3378  
changed lines
  Added in v.3379

  ViewVC Help
Powered by ViewVC 1.1.26