/[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 3203 by jfenwick, Thu Sep 23 23:24:09 2010 UTC revision 3204 by jfenwick, Thu Sep 23 23:59:39 2010 UTC
# Line 106  void Dudley_Assemble_LumpedSystem(Dudley Line 106  void Dudley_Assemble_LumpedSystem(Dudley
106    if (Dudley_noError()) {    if (Dudley_noError()) {
107      requireWrite(lumpedMat);      requireWrite(lumpedMat);
108      lumpedMat_p=getSampleDataRW(lumpedMat,0);      lumpedMat_p=getSampleDataRW(lumpedMat,0);
109      len_EM_lumpedMat=p.row_numShapesTotal*p.numEqu;      len_EM_lumpedMat=p.numShapes*p.numEqu;
110      len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);      len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
111            
112      expandedD=isExpanded(D);      expandedD=isExpanded(D);
# Line 123  void Dudley_Assemble_LumpedSystem(Dudley Line 123  void Dudley_Assemble_LumpedSystem(Dudley
123  #endif  #endif
124      {      {
125         EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);         EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
126         row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);         row_index=THREAD_MEMALLOC(p.numShapes,index_t);
127         if ( !Dudley_checkPtr(EM_lumpedMat) && !Dudley_checkPtr(row_index) ) {         if ( !Dudley_checkPtr(EM_lumpedMat) && !Dudley_checkPtr(row_index) ) {
128            if (p.numEqu == 1) { /* single equation */            if (p.numEqu == 1) { /* single equation */
129               if (expandedD) {   /* with expanded D */                       if (expandedD) {   /* with expanded D */        
# Line 140  void Dudley_Assemble_LumpedSystem(Dudley Line 140  void Dudley_Assemble_LumpedSystem(Dudley
140                                     for (q=0;q<p.numQuadTotal;q++) m_t+=vol*D_p[INDEX2(q, 0,p.numQuadTotal) ];                                     for (q=0;q<p.numQuadTotal;q++) m_t+=vol*D_p[INDEX2(q, 0,p.numQuadTotal) ];
141                                                        
142                                     diagS=0; /* diagonal sum: S */                                     diagS=0; /* diagonal sum: S */
143                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.numShapes;s++) {
144                                        rtmp=0;                                        rtmp=0;
145                                        for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*D_p[INDEX2(q, 0, p.numQuadTotal)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                                        for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*D_p[INDEX2(q, 0, p.numQuadTotal)]*S[INDEX2(s,q,p.numShapes)]*S[INDEX2(s,q,p.numShapes)];
146                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
147                                        diagS+=rtmp;                                        diagS+=rtmp;
148                                     }                                     }
149                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
150                     rtmp=m_t/diagS;                     rtmp=m_t/diagS;
151                                     for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;                                     for (s=0;s<p.numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
152                                                        
153                                 #else /* row-sum lumping */                                 #else /* row-sum lumping */
154                                     for (s=0;s<p.row_numShapes;s++)                                     for (s=0;s<p.numShapes;s++)
155                     {                     {
156                                         rtmp=0;                                         rtmp=0;
157                                         for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX2(q, 0, p.numQuadTotal)];                                         for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.numShapes)]*D_p[INDEX2(q, 0, p.numQuadTotal)];
158                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
159                                     }                                     }
160                                 #endif                                 #endif
161                                 for (q=0;q<p.row_numShapesTotal;q++)                                 for (q=0;q<p.numShapes;q++)
162                     {                     {
163                   row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];                   row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
164                     }                     }
165                                 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                                 Dudley_Util_AddScatter(p.numShapes,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
166                         } /* end color check */                                 } /* end color check */        
167                      } /* end element loop */                      } /* end element loop */
168                    } /* end color loop */                    } /* end color loop */
# Line 179  void Dudley_Assemble_LumpedSystem(Dudley Line 179  void Dudley_Assemble_LumpedSystem(Dudley
179                                     m_t=0; /* mass of the element: m_t */                                     m_t=0; /* mass of the element: m_t */
180                                     for (q=0;q<p.numQuadTotal;q++) m_t+=vol;                                     for (q=0;q<p.numQuadTotal;q++) m_t+=vol;
181                                     diagS=0; /* diagonal sum: S */                                     diagS=0; /* diagonal sum: S */
182                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.numShapes;s++) {
183                                        rtmp=0;                                        rtmp=0;
184                                        for (q=0;q<p.numQuadTotal;q++)                                        for (q=0;q<p.numQuadTotal;q++)
185                        {                        {
186                      rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                      rtmp+=vol*S[INDEX2(s,q,p.numShapes)]*S[INDEX2(s,q,p.numShapes)];
187                        }                        }
188                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
189                                        diagS+=rtmp;                                        diagS+=rtmp;
190                                     }                                     }
191                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
192                     rtmp=m_t/diagS*D_p[0];                     rtmp=m_t/diagS*D_p[0];
193                                     for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;                                     for (s=0;s<p.numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
194                                 #else /* row-sum lumping */                                 #else /* row-sum lumping */
195                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.numShapes;s++) {
196                                         rtmp=0;                                         rtmp=0;
197                                         for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)];                                         for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.numShapes)];
198                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];
199                                     }                                     }
200                                 #endif                                 #endif
201                                 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];                                 for (q=0;q<p.numShapes;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
202                                 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                                 Dudley_Util_AddScatter(p.numShapes,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
203                         } /* end color check */                         } /* end color check */
204                      } /* end element loop */                      } /* end element loop */
205                    } /* end color loop */                    } /* end color loop */
# Line 221  void Dudley_Assemble_LumpedSystem(Dudley Line 221  void Dudley_Assemble_LumpedSystem(Dudley
221                                        for (q=0;q<p.numQuadTotal;q++) m_t+=vol*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadTotal)];                                        for (q=0;q<p.numQuadTotal;q++) m_t+=vol*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadTotal)];
222                                                                                                        
223                                        diagS=0; /* diagonal sum: S */                                        diagS=0; /* diagonal sum: S */
224                                        for (s=0;s<p.row_numShapes;s++) {                                        for (s=0;s<p.numShapes;s++) {
225                                            rtmp=0;                                            rtmp=0;
226                                            for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadTotal)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                                            for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadTotal)]*S[INDEX2(s,q,p.numShapes)]*S[INDEX2(s,q,p.numShapes)];
227                                            EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;                                            EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
228                                            diagS+=rtmp;                                            diagS+=rtmp;
229                                         }                                         }
230                         /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                         /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
231                         rtmp=m_t/diagS;                         rtmp=m_t/diagS;
232                         for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;                         for (s=0;s<p.numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;
233                      }                                      }                
234                  #else /* row-sum lumping */                  #else /* row-sum lumping */
235                                    for (s=0;s<p.row_numShapes;s++) {                                    for (s=0;s<p.numShapes;s++) {
236                                        for (k=0;k<p.numEqu;k++) {                                        for (k=0;k<p.numEqu;k++) {
237                                           rtmp=0.;                                           rtmp=0.;
238                                            for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadTotal)];                                            for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.numShapes)]*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadTotal)];
239                                             EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;                                             EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
240                                         }                                         }
241                    }                    }
242                    #endif                    #endif
243                    for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];                    for (q=0;q<p.numShapes;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
244                    Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                    Dudley_Util_AddScatter(p.numShapes,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
245                         } /* end color check */                         } /* end color check */
246                      } /* end element loop */                      } /* end element loop */
247                  } /* end color loop */                  } /* end color loop */
# Line 257  void Dudley_Assemble_LumpedSystem(Dudley Line 257  void Dudley_Assemble_LumpedSystem(Dudley
257                                        m_t=0; /* mass of the element: m_t */                                        m_t=0; /* mass of the element: m_t */
258                                        for (q=0;q<p.numQuadTotal;q++) m_t+=vol;                                        for (q=0;q<p.numQuadTotal;q++) m_t+=vol;
259                                        diagS=0; /* diagonal sum: S */                                        diagS=0; /* diagonal sum: S */
260                                        for (s=0;s<p.row_numShapes;s++) {                                        for (s=0;s<p.numShapes;s++) {
261                                            rtmp=0;                                            rtmp=0;
262                                            for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                                            for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.numShapes)]*S[INDEX2(s,q,p.numShapes)];
263                        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;
264                                            diagS+=rtmp;                                            diagS+=rtmp;
265                                        }                                        }
266                        /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                        /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
267                        rtmp=m_t/diagS;                        rtmp=m_t/diagS;
268                        for (s=0;s<p.row_numShapes;s++) {                        for (s=0;s<p.numShapes;s++) {
269                          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];
270                        }                        }
271                     #else /* row-sum lumping */                     #else /* row-sum lumping */
272                                   for (s=0;s<p.row_numShapes;s++) {                                   for (s=0;s<p.numShapes;s++) {
273                      for (k=0;k<p.numEqu;k++) {                      for (k=0;k<p.numEqu;k++) {
274                                          rtmp=0.;                                          rtmp=0.;
275                                          for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)];                                          for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.numShapes)];
276                                          EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];                                          EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];
277                                       }                                       }
278                   }                   }
279                    #endif                    #endif
280                    for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];                    for (q=0;q<p.numShapes;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
281                    Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                    Dudley_Util_AddScatter(p.numShapes,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
282                         } /* end color check */                         } /* end color check */
283                      } /* end element loop */                      } /* end element loop */
284                  } /* end color loop */                  } /* end color loop */

Legend:
Removed from v.3203  
changed lines
  Added in v.3204

  ViewVC Help
Powered by ViewVC 1.1.26