/[escript]/branches/domexper/dudley/src/Assemble_LumpedSystem.c
ViewVC logotype

Diff of /branches/domexper/dudley/src/Assemble_LumpedSystem.c

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

revision 3135 by jfenwick, Thu Aug 5 05:07:58 2010 UTC revision 3136 by jfenwick, Thu Sep 2 03:48:29 2010 UTC
# Line 40  void Dudley_Assemble_LumpedSystem(Dudley Line 40  void Dudley_Assemble_LumpedSystem(Dudley
40    bool_t reducedIntegrationOrder=FALSE, expandedD;    bool_t reducedIntegrationOrder=FALSE, expandedD;
41    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
42    Assemble_Parameters p;    Assemble_Parameters p;
43    dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s, isub;    dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
44    type_t funcspace;    type_t funcspace;
45    index_t color,*row_index=NULL;    index_t color,*row_index=NULL;
46    __const double *D_p=NULL;    __const double *D_p=NULL;
# Line 111  void Dudley_Assemble_LumpedSystem(Dudley Line 111  void Dudley_Assemble_LumpedSystem(Dudley
111      S=p.row_jac->BasisFunctions->S;      S=p.row_jac->BasisFunctions->S;
112    
113  #ifdef NEW_LUMPING  #ifdef NEW_LUMPING
114      #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t, isub)      #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t)
115  #else  #else
116      #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, isub)      #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp)
117  #endif  #endif
118      {      {
119         EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);         EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
# Line 126  void Dudley_Assemble_LumpedSystem(Dudley Line 126  void Dudley_Assemble_LumpedSystem(Dudley
126                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
127                      for(e=0;e<elements->numElements;e++){                                    for(e=0;e<elements->numElements;e++){              
128              if (elements->Color[e]==color) {              if (elements->Color[e]==color) {
129                              for (isub=0; isub<p.numSub; isub++) {                                 Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadSub,1)]);
                                Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);  
130                                 D_p=getSampleDataRO(D,e);                                                           D_p=getSampleDataRO(D,e);                          
131                     #ifdef NEW_LUMPING /* HRZ lumping */                     #ifdef NEW_LUMPING /* HRZ lumping */
   
132                                     m_t=0; /* mass of the element: m_t */                                     m_t=0; /* mass of the element: m_t */
133                                     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, 0,p.numQuadSub) ];
134                                                        
135                                     diagS=0; /* diagonal sum: S */                                     diagS=0; /* diagonal sum: S */
136                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
137                                        rtmp=0;                                        rtmp=0;
138                                        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)];                                        for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*D_p[INDEX2(q, 0, p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
139                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
140                                        diagS+=rtmp;                                        diagS+=rtmp;
141                                     }                                     }
# Line 146  void Dudley_Assemble_LumpedSystem(Dudley Line 144  void Dudley_Assemble_LumpedSystem(Dudley
144                                     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;
145                                                        
146                                 #else /* row-sum lumping */                                 #else /* row-sum lumping */
147                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++)
148                       {
149                                         rtmp=0;                                         rtmp=0;
150                                         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, 0, p.numQuadSub)];
151                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
152                                     }                                     }
   
153                                 #endif                                 #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++)
155                       {
156                     row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,0, p.row_numShapesTotal)],e,p.NN)]];
157                       }
158                                 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                                 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
                 } /* end of isub loop */  
159                         } /* end color check */                                 } /* end color check */        
160                      } /* end element loop */                      } /* end element loop */
161                    } /* end color loop */                    } /* end color loop */
# Line 166  void Dudley_Assemble_LumpedSystem(Dudley Line 166  void Dudley_Assemble_LumpedSystem(Dudley
166                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
167                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){
168                          if (elements->Color[e]==color) {                          if (elements->Color[e]==color) {
169                              for (isub=0; isub<p.numSub; isub++) {                                 Vol=&(p.row_jac->volume[INDEX3(0, 0, e, p.numQuadSub,1)]);                
                                Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);                  
170                                 D_p=getSampleDataRO(D,e);                                                           D_p=getSampleDataRO(D,e);                          
171                     #ifdef NEW_LUMPING /* HRZ lumping */                     #ifdef NEW_LUMPING /* HRZ lumping */
172                                     m_t=0; /* mass of the element: m_t */                                     m_t=0; /* mass of the element: m_t */
# Line 175  void Dudley_Assemble_LumpedSystem(Dudley Line 174  void Dudley_Assemble_LumpedSystem(Dudley
174                                     diagS=0; /* diagonal sum: S */                                     diagS=0; /* diagonal sum: S */
175                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
176                                        rtmp=0;                                        rtmp=0;
177                                        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++)
178                          {
179                        rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
180                          }
181                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
182                                        diagS+=rtmp;                                        diagS+=rtmp;
183                                     }                                     }
# Line 189  void Dudley_Assemble_LumpedSystem(Dudley Line 191  void Dudley_Assemble_LumpedSystem(Dudley
191                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];
192                                     }                                     }
193                                 #endif                                 #endif
194                                 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, 0, p.row_numShapesTotal)],e,p.NN)]];
195                                 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                                 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
                             } /* end of isub loop */  
196                         } /* end color check */                         } /* end color check */
197                      } /* end element loop */                      } /* end element loop */
198                    } /* end color loop */                    } /* end color loop */
# Line 204  void Dudley_Assemble_LumpedSystem(Dudley Line 205  void Dudley_Assemble_LumpedSystem(Dudley
205                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
206                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){
207                         if (elements->Color[e]==color) {                         if (elements->Color[e]==color) {
208                for (isub=0; isub<p.numSub; isub++) {                                    Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadSub,1)]);
                               Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);  
209                                D_p=getSampleDataRO(D,e);                                  D_p=getSampleDataRO(D,e);  
210                                    
211                                #ifdef NEW_LUMPING /* HRZ lumping */                                #ifdef NEW_LUMPING /* HRZ lumping */
212                                    for (k=0;k<p.numEqu;k++) {                                    for (k=0;k<p.numEqu;k++) {
213                                        m_t=0; /* mass of the element: m_t */                                        m_t=0; /* mass of the element: m_t */
214                                        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,0,p.numEqu,p.numQuadSub)];
215                                                                                                        
216                                        diagS=0; /* diagonal sum: S */                                        diagS=0; /* diagonal sum: S */
217                                        for (s=0;s<p.row_numShapes;s++) {                                        for (s=0;s<p.row_numShapes;s++) {
218                                            rtmp=0;                                            rtmp=0;
219                                            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,0,p.numEqu,p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
220                                            EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;                                            EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
221                                            diagS+=rtmp;                                            diagS+=rtmp;
222                                         }                                         }
# Line 228  void Dudley_Assemble_LumpedSystem(Dudley Line 228  void Dudley_Assemble_LumpedSystem(Dudley
228                                    for (s=0;s<p.row_numShapes;s++) {                                    for (s=0;s<p.row_numShapes;s++) {
229                                        for (k=0;k<p.numEqu;k++) {                                        for (k=0;k<p.numEqu;k++) {
230                                           rtmp=0.;                                           rtmp=0.;
231                                            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,0,p.numEqu,p.numQuadSub)];
232                                             EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;                                             EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
233                                         }                                         }
234                    }                    }
235                    #endif                    #endif
236                    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,0,p.row_numShapesTotal)],e,p.NN)]];
237                    Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                    Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
                } /* end of isub loop */  
238                         } /* end color check */                         } /* end color check */
239                      } /* end element loop */                      } /* end element loop */
240                  } /* end color loop */                  } /* end color loop */
# Line 245  void Dudley_Assemble_LumpedSystem(Dudley Line 244  void Dudley_Assemble_LumpedSystem(Dudley
244                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
245                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){
246                         if (elements->Color[e]==color) {                         if (elements->Color[e]==color) {
247                 for (isub=0; isub<p.numSub; isub++) {                                    Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadSub,1)]);
                               Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);  
248                                D_p=getSampleDataRO(D,e);                                D_p=getSampleDataRO(D,e);
249                                
250                                #ifdef NEW_LUMPING /* HRZ lumping */                                #ifdef NEW_LUMPING /* HRZ lumping */
# Line 273  void Dudley_Assemble_LumpedSystem(Dudley Line 271  void Dudley_Assemble_LumpedSystem(Dudley
271                                       }                                       }
272                   }                   }
273                    #endif                    #endif
274                    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,0,p.row_numShapesTotal)],e,p.NN)]];
275                    Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                    Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
               } /* end of isub loop */  
276                         } /* end color check */                         } /* end color check */
277                      } /* end element loop */                      } /* end element loop */
278                  } /* end color loop */                  } /* end color loop */

Legend:
Removed from v.3135  
changed lines
  Added in v.3136

  ViewVC Help
Powered by ViewVC 1.1.26