/[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 3136 by jfenwick, Thu Sep 2 03:48:29 2010 UTC revision 3137 by jfenwick, Thu Sep 2 05:03:17 2010 UTC
# Line 80  void Dudley_Assemble_LumpedSystem(Dudley Line 80  void Dudley_Assemble_LumpedSystem(Dudley
80    
81    /* check if all function spaces are the same */    /* check if all function spaces are the same */
82    if (! numSamplesEqual(D,p.numQuadTotal,elements->numElements) ) {    if (! numSamplesEqual(D,p.numQuadTotal,elements->numElements) ) {
83          sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuadSub,elements->numElements);          sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuadTotal,elements->numElements);
84          Dudley_setError(TYPE_ERROR,error_msg);          Dudley_setError(TYPE_ERROR,error_msg);
85    }    }
86    
# 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                                 Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadSub,1)]);                                 Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);
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, 0,p.numQuadSub) ];                                     for (q=0;q<p.numQuadTotal;q++) m_t+=Vol[q]*D_p[INDEX2(q, 0,p.numQuadTotal) ];
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, 0, p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                                        for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*D_p[INDEX2(q, 0, p.numQuadTotal)]*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 147  void Dudley_Assemble_LumpedSystem(Dudley Line 147  void Dudley_Assemble_LumpedSystem(Dudley
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, 0, p.numQuadSub)];                                         for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX2(q, 0, p.numQuadTotal)];
151                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
152                                     }                                     }
153                                 #endif                                 #endif
# 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                                 Vol=&(p.row_jac->volume[INDEX3(0, 0, e, p.numQuadSub,1)]);                                                 Vol=&(p.row_jac->volume[INDEX3(0, 0, e, p.numQuadTotal,1)]);              
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 */
173                                     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];                                     for (q=0;q<p.numQuadTotal;q++) m_t+=Vol[q];
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++)                                        for (q=0;q<p.numQuadTotal;q++)
178                        {                        {
179                      rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                      rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
180                        }                        }
# Line 187  void Dudley_Assemble_LumpedSystem(Dudley Line 187  void Dudley_Assemble_LumpedSystem(Dudley
187                                 #else /* row-sum lumping */                                 #else /* row-sum lumping */
188                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
189                                         rtmp=0;                                         rtmp=0;
190                                         for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];                                         for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
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
# Line 205  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                                Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadSub,1)]);                                Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);
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,0,p.numEqu,p.numQuadSub)];                                        for (q=0;q<p.numQuadTotal;q++) m_t+=Vol[q]*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadTotal)];
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,0,p.numEqu,p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                                            for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadTotal)]*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,0,p.numEqu,p.numQuadSub)];                                            for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadTotal)];
232                                             EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;                                             EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
233                                         }                                         }
234                    }                    }
# Line 244  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                                Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadSub,1)]);                                Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);
248                                D_p=getSampleDataRO(D,e);                                D_p=getSampleDataRO(D,e);
249                                
250                                #ifdef NEW_LUMPING /* HRZ lumping */                                #ifdef NEW_LUMPING /* HRZ lumping */
251                                        m_t=0; /* mass of the element: m_t */                                        m_t=0; /* mass of the element: m_t */
252                                        for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];                                        for (q=0;q<p.numQuadTotal;q++) m_t+=Vol[q];
253                                        diagS=0; /* diagonal sum: S */                                        diagS=0; /* diagonal sum: S */
254                                        for (s=0;s<p.row_numShapes;s++) {                                        for (s=0;s<p.row_numShapes;s++) {
255                                            rtmp=0;                                            rtmp=0;
256                                            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.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
257                        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;
258                                            diagS+=rtmp;                                            diagS+=rtmp;
259                                        }                                        }
# Line 266  void Dudley_Assemble_LumpedSystem(Dudley Line 266  void Dudley_Assemble_LumpedSystem(Dudley
266                                   for (s=0;s<p.row_numShapes;s++) {                                   for (s=0;s<p.row_numShapes;s++) {
267                      for (k=0;k<p.numEqu;k++) {                      for (k=0;k<p.numEqu;k++) {
268                                          rtmp=0.;                                          rtmp=0.;
269                                          for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];                                          for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
270                                          EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];                                          EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];
271                                       }                                       }
272                   }                   }

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

  ViewVC Help
Powered by ViewVC 1.1.26