/[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 3183 by jfenwick, Fri Sep 3 00:49:02 2010 UTC revision 3184 by jfenwick, Wed Sep 15 00:23:42 2010 UTC
# Line 44  void Dudley_Assemble_LumpedSystem(Dudley Line 44  void Dudley_Assemble_LumpedSystem(Dudley
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;
47    double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *lumpedMat_p=NULL;    double *S=NULL, *EM_lumpedMat=NULL, *lumpedMat_p=NULL;
48    register double rtmp;    register double rtmp;
49    size_t len_EM_lumpedMat_size;    size_t len_EM_lumpedMat_size;
50  #ifdef NEW_LUMPING  #ifdef NEW_LUMPING
# 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.numQuadTotal,1)]);  //                               Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);
130                    double vol=p.row_jac->absD[e]*p.row_jac->quadweight;
131                                 D_p=getSampleDataRO(D,e);                                                           D_p=getSampleDataRO(D,e);                          
132                     #ifdef NEW_LUMPING /* HRZ lumping */                     #ifdef NEW_LUMPING /* HRZ lumping */
133                                     m_t=0; /* mass of the element: m_t */                                     m_t=0; /* mass of the element: m_t */
134                                     for (q=0;q<p.numQuadTotal;q++) m_t+=Vol[q]*D_p[INDEX2(q, 0,p.numQuadTotal) ];                                     for (q=0;q<p.numQuadTotal;q++) m_t+=vol*D_p[INDEX2(q, 0,p.numQuadTotal) ];
135                                                        
136                                     diagS=0; /* diagonal sum: S */                                     diagS=0; /* diagonal sum: S */
137                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
138                                        rtmp=0;                                        rtmp=0;
139                                        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)];                                        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)];
140                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
141                                        diagS+=rtmp;                                        diagS+=rtmp;
142                                     }                                     }
# Line 147  void Dudley_Assemble_LumpedSystem(Dudley Line 148  void Dudley_Assemble_LumpedSystem(Dudley
148                                     for (s=0;s<p.row_numShapes;s++)                                     for (s=0;s<p.row_numShapes;s++)
149                     {                     {
150                                         rtmp=0;                                         rtmp=0;
151                                         for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*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.row_numShapes)]*D_p[INDEX2(q, 0, p.numQuadTotal)];
152                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
153                                     }                                     }
154                                 #endif                                 #endif
# Line 166  void Dudley_Assemble_LumpedSystem(Dudley Line 167  void Dudley_Assemble_LumpedSystem(Dudley
167                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
168                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){
169                          if (elements->Color[e]==color) {                          if (elements->Color[e]==color) {
170                                 Vol=&(p.row_jac->volume[INDEX3(0, 0, e, p.numQuadTotal,1)]);                                double vol=p.row_jac->absD[e]*p.row_jac->quadweight;
171                                 D_p=getSampleDataRO(D,e);                                                           D_p=getSampleDataRO(D,e);                          
172                     #ifdef NEW_LUMPING /* HRZ lumping */                     #ifdef NEW_LUMPING /* HRZ lumping */
173                                     m_t=0; /* mass of the element: m_t */                                     m_t=0; /* mass of the element: m_t */
174                                     for (q=0;q<p.numQuadTotal;q++) m_t+=Vol[q];                                     for (q=0;q<p.numQuadTotal;q++) m_t+=vol;
175                                     diagS=0; /* diagonal sum: S */                                     diagS=0; /* diagonal sum: S */
176                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
177                                        rtmp=0;                                        rtmp=0;
178                                        for (q=0;q<p.numQuadTotal;q++)                                        for (q=0;q<p.numQuadTotal;q++)
179                        {                        {
180                      rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];                      rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
181                        }                        }
182                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;                                        EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
183                                        diagS+=rtmp;                                        diagS+=rtmp;
# Line 187  void Dudley_Assemble_LumpedSystem(Dudley Line 188  void Dudley_Assemble_LumpedSystem(Dudley
188                                 #else /* row-sum lumping */                                 #else /* row-sum lumping */
189                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
190                                         rtmp=0;                                         rtmp=0;
191                                         for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];                                         for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*S[INDEX2(s,q,p.row_numShapes)];
192                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];
193                                     }                                     }
194                                 #endif                                 #endif
# Line 205  void Dudley_Assemble_LumpedSystem(Dudley Line 206  void Dudley_Assemble_LumpedSystem(Dudley
206                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
207                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){
208                         if (elements->Color[e]==color) {                         if (elements->Color[e]==color) {
209                                Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);                  double vol=p.row_jac->absD[e]*p.row_jac->quadweight;
210                                D_p=getSampleDataRO(D,e);                                  D_p=getSampleDataRO(D,e);  
211                                    
212                                #ifdef NEW_LUMPING /* HRZ lumping */                                #ifdef NEW_LUMPING /* HRZ lumping */
213                                    for (k=0;k<p.numEqu;k++) {                                    for (k=0;k<p.numEqu;k++) {
214                                        m_t=0; /* mass of the element: m_t */                                        m_t=0; /* mass of the element: m_t */
215                                        for (q=0;q<p.numQuadTotal;q++) m_t+=Vol[q]*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)];
216                                                                                                        
217                                        diagS=0; /* diagonal sum: S */                                        diagS=0; /* diagonal sum: S */
218                                        for (s=0;s<p.row_numShapes;s++) {                                        for (s=0;s<p.row_numShapes;s++) {
219                                            rtmp=0;                                            rtmp=0;
220                                            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)];                                            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)];
221                                            EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;                                            EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
222                                            diagS+=rtmp;                                            diagS+=rtmp;
223                                         }                                         }
# Line 228  void Dudley_Assemble_LumpedSystem(Dudley Line 229  void Dudley_Assemble_LumpedSystem(Dudley
229                                    for (s=0;s<p.row_numShapes;s++) {                                    for (s=0;s<p.row_numShapes;s++) {
230                                        for (k=0;k<p.numEqu;k++) {                                        for (k=0;k<p.numEqu;k++) {
231                                           rtmp=0.;                                           rtmp=0.;
232                                            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)];                                            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)];
233                                             EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;                                             EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
234                                         }                                         }
235                    }                    }
# Line 244  void Dudley_Assemble_LumpedSystem(Dudley Line 245  void Dudley_Assemble_LumpedSystem(Dudley
245                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
246                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){
247                         if (elements->Color[e]==color) {                         if (elements->Color[e]==color) {
248                                Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);                  double vol=p.row_jac->absD[e]*p.row_jac->quadweight;                              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.numQuadTotal;q++) m_t+=Vol[q];                                        for (q=0;q<p.numQuadTotal;q++) m_t+=vol;
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.numQuadTotal;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*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.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];                                          for (q=0;q<p.numQuadTotal;q++) rtmp+=vol*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.3183  
changed lines
  Added in v.3184

  ViewVC Help
Powered by ViewVC 1.1.26