/[escript]/trunk/dudley/src/Assemble_LumpedSystem.c
ViewVC logotype

Diff of /trunk/dudley/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    
30  #include "ShapeTable.h"  #include "ShapeTable.h"
31    
 /* Disabled until the tests pass */  
 /* #define NEW_LUMPING */  
   
32  /**************************************************************/  /**************************************************************/
33    
34  void Dudley_Assemble_LumpedSystem(Dudley_NodeFile * nodes, Dudley_ElementFile * elements, escriptDataC * lumpedMat,  void Dudley_Assemble_LumpedSystem(Dudley_NodeFile * nodes, Dudley_ElementFile * elements, escriptDataC * lumpedMat,
35                    escriptDataC * D)                    escriptDataC * D, const bool_t useHRZ)
36  {  {
37    
38      bool_t reducedIntegrationOrder = FALSE, expandedD;      bool_t reducedIntegrationOrder = FALSE, expandedD;
# Line 49  void Dudley_Assemble_LumpedSystem(Dudley Line 46  void Dudley_Assemble_LumpedSystem(Dudley
46      double *EM_lumpedMat = NULL, *lumpedMat_p = NULL;      double *EM_lumpedMat = NULL, *lumpedMat_p = NULL;
47      register double rtmp;      register double rtmp;
48      size_t len_EM_lumpedMat_size;      size_t len_EM_lumpedMat_size;
 #ifdef NEW_LUMPING  
49      register double m_t = 0., diagS = 0.;      register double m_t = 0., diagS = 0.;
50  #endif  
51    
52      Dudley_resetError();      Dudley_resetError();
53    
# Line 138  void Dudley_Assemble_LumpedSystem(Dudley Line 134  void Dudley_Assemble_LumpedSystem(Dudley
134      {      {
135          Dudley_setError(TYPE_ERROR, "Dudley_Assemble_LumpedSystem: Unable to locate shape function.");          Dudley_setError(TYPE_ERROR, "Dudley_Assemble_LumpedSystem: Unable to locate shape function.");
136      }      }
137  #ifdef NEW_LUMPING      #pragma omp parallel private(color, EM_lumpedMat, row_index, D_p, s, q, k, rtmp, diagS, m_t)
 #pragma omp parallel private(color, EM_lumpedMat, row_index, D_p, s, q, k, rtmp, diagS, m_t)  
 #else  
 #pragma omp parallel private(color, EM_lumpedMat, row_index, D_p, s, q, k, rtmp)  
 #endif  
138      {      {
139          EM_lumpedMat = THREAD_MEMALLOC(len_EM_lumpedMat, double);          EM_lumpedMat = THREAD_MEMALLOC(len_EM_lumpedMat, double);
140          row_index = THREAD_MEMALLOC(p.numShapes, index_t);          row_index = THREAD_MEMALLOC(p.numShapes, index_t);
# Line 162  void Dudley_Assemble_LumpedSystem(Dudley Line 154  void Dudley_Assemble_LumpedSystem(Dudley
154                  {                  {
155                      double vol = p.row_jac->absD[e] * p.row_jac->quadweight;                      double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
156                      D_p = getSampleDataRO(D, e);                      D_p = getSampleDataRO(D, e);
157  #ifdef NEW_LUMPING      /* HRZ lumping */                      if (useHRZ) {
158                      m_t = 0;    /* mass of the element: m_t */                         m_t = 0; /* mass of the element: m_t */
159                      for (q = 0; q < p.numQuad; q++)                         for (q = 0; q < p.numQuad; q++)
160                      m_t += vol * D_p[INDEX2(q, 0, p.numQuad)];                        m_t += vol * D_p[INDEX2(q, 0, p.numQuad)];
161                      diagS = 0;  /* diagonal sum: S */                         diagS = 0;   /* diagonal sum: S */
162                      for (s = 0; s < p.numShapes; s++)                         for (s = 0; s < p.numShapes; s++)
163                      {                         {
164                      rtmp = 0;                        rtmp = 0;
165                      for (q = 0; q < p.numQuad; q++)                        for (q = 0; q < p.numQuad; q++)
166                          rtmp +=                           rtmp +=
167                          vol * D_p[INDEX2(q, 0, p.numQuad)] * S[INDEX2(s, q, p.numShapes)] *                             vol * D_p[INDEX2(q, 0, p.numQuad)] * S[INDEX2(s, q, p.numShapes)] *
168                          S[INDEX2(s, q, p.numShapes)];                             S[INDEX2(s, q, p.numShapes)];
169                      EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;                        EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
170                      diagS += rtmp;                        diagS += rtmp;
171                      }                         }
172                      /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                         /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
173                      rtmp = m_t / diagS;                         rtmp = m_t / diagS;
174                      for (s = 0; s < p.numShapes; s++)                         for (s = 0; s < p.numShapes; s++)
175                      EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;                        EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;
176    
177  #else               /* row-sum lumping */                      } else {/* row-sum lumping */
178                      for (s = 0; s < p.numShapes; s++)                         for (s = 0; s < p.numShapes; s++)
179                      {                         {
180                      rtmp = 0;                        rtmp = 0;
181                      for (q = 0; q < p.numQuad; q++)                        for (q = 0; q < p.numQuad; q++)
182                          rtmp += vol * S[INDEX2(s, q, p.numShapes)] * D_p[INDEX2(q, 0, p.numQuad)];                           rtmp += vol * S[INDEX2(s, q, p.numShapes)] * D_p[INDEX2(q, 0, p.numQuad)];
183                      EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;                        EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
184                           }
185                      }                      }
 #endif  
186                      for (q = 0; q < p.numShapes; q++)                      for (q = 0; q < p.numShapes; q++)
187                      {                      {
188                      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)]];
# Line 214  void Dudley_Assemble_LumpedSystem(Dudley Line 206  void Dudley_Assemble_LumpedSystem(Dudley
206                  {                  {
207                      double vol = p.row_jac->absD[e] * p.row_jac->quadweight;                      double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
208                      D_p = getSampleDataRO(D, e);                      D_p = getSampleDataRO(D, e);
209  #ifdef NEW_LUMPING      /* HRZ lumping */                      if (useHRZ) {   /* HRZ lumping */
210                      m_t = 0;    /* mass of the element: m_t */                         m_t = 0; /* mass of the element: m_t */
211                      for (q = 0; q < p.numQuad; q++)                         for (q = 0; q < p.numQuad; q++)
212                      m_t += vol;                        m_t += vol;
213                      diagS = 0;  /* diagonal sum: S */                         diagS = 0;   /* diagonal sum: S */
214                      for (s = 0; s < p.numShapes; s++)                         for (s = 0; s < p.numShapes; s++)
215                      {                         {
216                      rtmp = 0;                        rtmp = 0;
217                      for (q = 0; q < p.numQuad; q++)                        for (q = 0; q < p.numQuad; q++)
218                      {                        {
219                          rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];                           rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
220                      }                        }
221                      EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;                        EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
222                      diagS += rtmp;                        diagS += rtmp;
223                           }
224                           /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
225                           rtmp = m_t / diagS * D_p[0];
226                           for (s = 0; s < p.numShapes; s++)
227                          EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;
228                        } else {            /* row-sum lumping */
229                           for (s = 0; s < p.numShapes; s++)
230                           {
231                          rtmp = 0;
232                          for (q = 0; q < p.numQuad; q++)
233                             rtmp += vol * S[INDEX2(s, q, p.numShapes)];
234                          EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp * D_p[0];
235                           }
236                      }                      }
                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */  
                     rtmp = m_t / diagS * D_p[0];  
                     for (s = 0; s < p.numShapes; s++)  
                     EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;  
 #else               /* row-sum lumping */  
                     for (s = 0; s < p.numShapes; s++)  
                     {  
                     rtmp = 0;  
                     for (q = 0; q < p.numQuad; q++)  
                         rtmp += vol * S[INDEX2(s, q, p.numShapes)];  
                     EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp * D_p[0];  
                     }  
 #endif  
237                      for (q = 0; q < p.numShapes; q++)                      for (q = 0; q < p.numShapes; q++)
238                      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)]];
239                      Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,                      Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,
# Line 267  void Dudley_Assemble_LumpedSystem(Dudley Line 259  void Dudley_Assemble_LumpedSystem(Dudley
259                      double vol = p.row_jac->absD[e] * p.row_jac->quadweight;                      double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
260                      D_p = getSampleDataRO(D, e);                      D_p = getSampleDataRO(D, e);
261    
262  #ifdef NEW_LUMPING      /* HRZ lumping */                      if (useHRZ) {   /* HRZ lumping */
263                      for (k = 0; k < p.numEqu; k++)                         for (k = 0; k < p.numEqu; k++)
264                      {                         {
265                      m_t = 0;    /* mass of the element: m_t */                        m_t = 0;  /* mass of the element: m_t */
266                      for (q = 0; q < p.numQuad; q++)                        for (q = 0; q < p.numQuad; q++)
267                          m_t += vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];                           m_t += vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];
268    
269                      diagS = 0;  /* diagonal sum: S */                        diagS = 0;    /* diagonal sum: S */
270                      for (s = 0; s < p.numShapes; s++)                        for (s = 0; s < p.numShapes; s++)
271                      {                        {
272                          rtmp = 0;                           rtmp = 0;
273                          for (q = 0; q < p.numQuad; q++)                           for (q = 0; q < p.numQuad; q++)
274                          rtmp +=                             rtmp +=
275                              vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)] *                                vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)] *
276                              S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];                                S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
277                          EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;                           EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
278                          diagS += rtmp;                           diagS += rtmp;
279                      }                        }
280                      /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                        /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
281                      rtmp = m_t / diagS;                        rtmp = m_t / diagS;
282                      for (s = 0; s < p.numShapes; s++)                        for (s = 0; s < p.numShapes; s++)
283                          EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp;                           EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp;
284                      }                         }
285  #else               /* row-sum lumping */                      } else {                /* row-sum lumping */
286                      for (s = 0; s < p.numShapes; s++)                         for (s = 0; s < p.numShapes; s++)
287                      {                         {
288                      for (k = 0; k < p.numEqu; k++)                        for (k = 0; k < p.numEqu; k++)
289                      {                        {
290                          rtmp = 0.;                           rtmp = 0.;
291                          for (q = 0; q < p.numQuad; q++)                           for (q = 0; q < p.numQuad; q++)
292                          rtmp +=                             rtmp +=
293                              vol * S[INDEX2(s, q, p.numShapes)] *                                vol * S[INDEX2(s, q, p.numShapes)] *
294                              D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];                                D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];
295                          EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;                           EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
296                      }                        }
297                           }
298                      }                      }
 #endif  
299                      for (q = 0; q < p.numShapes; q++)                      for (q = 0; q < p.numShapes; q++)
300                      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)]];
301                      Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,                      Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,
# Line 325  void Dudley_Assemble_LumpedSystem(Dudley Line 317  void Dudley_Assemble_LumpedSystem(Dudley
317                      double vol = p.row_jac->absD[e] * p.row_jac->quadweight;                      double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
318                      D_p = getSampleDataRO(D, e);                      D_p = getSampleDataRO(D, e);
319    
320  #ifdef NEW_LUMPING      /* HRZ lumping */                      if (useHRZ)     { /* HRZ lumping */
321                      m_t = 0;    /* mass of the element: m_t */                         m_t = 0; /* mass of the element: m_t */
322                      for (q = 0; q < p.numQuad; q++)                         for (q = 0; q < p.numQuad; q++)
323                      m_t += vol;                        m_t += vol;
324                      diagS = 0;  /* diagonal sum: S */                         diagS = 0;   /* diagonal sum: S */
325                      for (s = 0; s < p.numShapes; s++)                         for (s = 0; s < p.numShapes; s++)
326                      {                         {
327                      rtmp = 0;                        rtmp = 0;
328                      for (q = 0; q < p.numQuad; q++)                        for (q = 0; q < p.numQuad; q++)
329                          rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];                           rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
330                      for (k = 0; k < p.numEqu; k++)                        for (k = 0; k < p.numEqu; k++)
331                          EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;                           EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
332                      diagS += rtmp;                        diagS += rtmp;
333                      }                         }
334                      /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                         /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
335                      rtmp = m_t / diagS;                         rtmp = m_t / diagS;
336                      for (s = 0; s < p.numShapes; s++)                         for (s = 0; s < p.numShapes; s++)
337                      {                         {
338                      for (k = 0; k < p.numEqu; k++)                        for (k = 0; k < p.numEqu; k++)
339                          EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp * D_p[k];                           EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp * D_p[k];
340                      }                         }
341  #else               /* row-sum lumping */                      } else {                /* row-sum lumping */
342                      for (s = 0; s < p.numShapes; s++)                         for (s = 0; s < p.numShapes; s++)
343                      {                         {
344                      for (k = 0; k < p.numEqu; k++)                        for (k = 0; k < p.numEqu; k++)
345                      {                        {
346                          rtmp = 0.;                          rtmp = 0.;
347                          for (q = 0; q < p.numQuad; q++)                          for (q = 0; q < p.numQuad; q++)
348                          rtmp += vol * S[INDEX2(s, q, p.numShapes)];                          rtmp += vol * S[INDEX2(s, q, p.numShapes)];
349                          EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp * D_p[k];                          EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp * D_p[k];
350                      }                        }
351                           }
352                      }                      }
 #endif  
353                      for (q = 0; q < p.numShapes; q++)                      for (q = 0; q < p.numShapes; q++)
354                      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)]];
355                      Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,                      Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,

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

  ViewVC Help
Powered by ViewVC 1.1.26