/[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 3521 by caltinay, Wed Mar 30 02:24:33 2011 UTC revision 3522 by gross, Tue May 24 00:57:58 2011 UTC
# Line 125  void Dudley_Assemble_LumpedSystem(Dudley Line 125  void Dudley_Assemble_LumpedSystem(Dudley
125      {      {
126      requireWrite(lumpedMat);      requireWrite(lumpedMat);
127      lumpedMat_p = getSampleDataRW(lumpedMat, 0);      lumpedMat_p = getSampleDataRW(lumpedMat, 0);
128      len_EM_lumpedMat = p.numShapes * p.numEqu;      
129        if (funcspace==DUDLEY_POINTS) {
130                  #pragma omp parallel private(color, D_p)
131              {
132                for (color=elements->minColor;color<=elements->maxColor;color++) {
133                  /*  open loop over all elements: */
134                  #pragma omp for private(e) schedule(static)
135                  for(e=0;e<elements->numElements;e++){
136                  if (elements->Color[e]==color) {
137                    D_p=getSampleDataRO(D, e);
138                    if (NULL!=D_p)  Dudley_Util_AddScatter(1,
139                                                            &(p.row_DOF[elements->Nodes[INDEX2(0,e,p.NN)]]),
140                                                            p.numEqu,
141                                                            D_p,
142                                                            lumpedMat_p,
143                                                            p.row_DOF_UpperBound);
144                  } /* end color check */
145                  } /* end element loop */
146              } /* end color loop */
147            } /* end parallel region */
148            } else {  
149              
150              len_EM_lumpedMat = p.numShapes * p.numEqu;
151    
152      expandedD = isExpanded(D);            expandedD = isExpanded(D);
153      if (!getQuadShape(elements->numDim, reducedIntegrationOrder, &S))            if (!getQuadShape(elements->numDim, reducedIntegrationOrder, &S))
154      {            {
155          Dudley_setError(TYPE_ERROR, "Dudley_Assemble_LumpedSystem: Unable to locate shape function.");            Dudley_setError(TYPE_ERROR, "Dudley_Assemble_LumpedSystem: Unable to locate shape function.");
156      }            }
157      #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)
158      {            {
159          EM_lumpedMat = THREAD_MEMALLOC(len_EM_lumpedMat, double);            EM_lumpedMat = THREAD_MEMALLOC(len_EM_lumpedMat, double);
160          row_index = THREAD_MEMALLOC(p.numShapes, index_t);            row_index = THREAD_MEMALLOC(p.numShapes, index_t);
161          if (!Dudley_checkPtr(EM_lumpedMat) && !Dudley_checkPtr(row_index))            if (!Dudley_checkPtr(EM_lumpedMat) && !Dudley_checkPtr(row_index))
162          {            {
163          if (p.numEqu == 1)                if (p.numEqu == 1)
164          {       /* single equation */                {     /* single equation */
165              if (expandedD)                if (expandedD)
166              {       /* with expanded D */                {     /* with expanded D */
167              for (color = elements->minColor; color <= elements->maxColor; color++)                    for (color = elements->minColor; color <= elements->maxColor; color++)
168              {                    {
169                  /*  open loop over all elements: */                    /*  open loop over all elements: */
170  #pragma omp for private(e) schedule(static)        #pragma omp for private(e) schedule(static)
171                  for (e = 0; e < elements->numElements; e++)                    for (e = 0; e < elements->numElements; e++)
172                  {                    {
173                  if (elements->Color[e] == color)                        if (elements->Color[e] == color)
174                  {                        {
175                      double vol = p.row_jac->absD[e] * p.row_jac->quadweight;                        double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
176                      D_p = getSampleDataRO(D, e);                        D_p = getSampleDataRO(D, e);
177                      if (useHRZ) {                        if (useHRZ)   {
178                         m_t = 0; /* mass of the element: m_t */                          m_t = 0;    /* mass of the element: m_t */
179                         for (q = 0; q < p.numQuad; q++)                          for (q = 0; q < p.numQuad; q++)
180                        m_t += vol * D_p[INDEX2(q, 0, p.numQuad)];                          m_t += vol * D_p[INDEX2(q, 0, p.numQuad)];
181                         diagS = 0;   /* diagonal sum: S */                          diagS = 0;  /* diagonal sum: S */
182                         for (s = 0; s < p.numShapes; s++)                          for (s = 0; s < p.numShapes; s++)
183                         {                          {
184                        rtmp = 0;                          rtmp = 0;
185                        for (q = 0; q < p.numQuad; q++)                          for (q = 0; q < p.numQuad; q++)
186                           rtmp +=                            rtmp +=
187                             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)] *
188                             S[INDEX2(s, q, p.numShapes)];                              S[INDEX2(s, q, p.numShapes)];
189                        EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;                          EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
190                        diagS += rtmp;                          diagS += rtmp;
191                         }                          }
192                         /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                          /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
193                         rtmp = m_t / diagS;                          rtmp = m_t / diagS;
194                         for (s = 0; s < p.numShapes; s++)                          for (s = 0; s < p.numShapes; s++)
195                        EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;                          EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;
196    
197                      } else {/* row-sum lumping */                        } else {/* row-sum lumping */
198                         for (s = 0; s < p.numShapes; s++)                          for (s = 0; s < p.numShapes; s++)
199                         {                          {
200                        rtmp = 0;                          rtmp = 0;
201                        for (q = 0; q < p.numQuad; q++)                          for (q = 0; q < p.numQuad; q++)
202                           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)];
203                        EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;                          EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
204                         }                          }
205                      }                        }
206                      for (q = 0; q < p.numShapes; q++)                        for (q = 0; q < p.numShapes; q++)
                     {  
                     row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];  
                     }  
                     Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,  
                                p.row_DOF_UpperBound);  
                 }   /* end color check */  
                 }   /* end element loop */  
             }   /* end color loop */  
             }  
             else  
             {       /* with constant D */  
   
             for (color = elements->minColor; color <= elements->maxColor; color++)  
             {  
                 /*  open loop over all elements: */  
 #pragma omp for private(e) schedule(static)  
                 for (e = 0; e < elements->numElements; e++)  
                 {  
                 if (elements->Color[e] == color)  
                 {  
                     double vol = p.row_jac->absD[e] * p.row_jac->quadweight;  
                     D_p = getSampleDataRO(D, e);  
                     if (useHRZ) {   /* HRZ lumping */  
                        m_t = 0; /* mass of the element: m_t */  
                        for (q = 0; q < p.numQuad; q++)  
                       m_t += vol;  
                        diagS = 0;   /* diagonal sum: S */  
                        for (s = 0; s < p.numShapes; s++)  
                        {  
                       rtmp = 0;  
                       for (q = 0; q < p.numQuad; q++)  
207                        {                        {
208                           rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];                            row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
209                        }                        }
210                        EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;                        Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,
211                        diagS += rtmp;                                  p.row_DOF_UpperBound);
212                         }                        } /* end color check */
213                         /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                    } /* end element loop */
214                         rtmp = m_t / diagS * D_p[0];                    } /* end color loop */
215                         for (s = 0; s < p.numShapes; s++)                }
216                        EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;                else
217                      } else {            /* row-sum lumping */                {     /* with constant D */
                        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];  
                        }  
                     }  
                     for (q = 0; q < p.numShapes; q++)  
                     row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];  
                     Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,  
                                p.row_DOF_UpperBound);  
                 }   /* end color check */  
                 }   /* end element loop */  
             }   /* end color loop */  
   
             }  
         }  
         else  
         {       /* system of  equation */  
             if (expandedD)  
             {       /* with expanded D */  
             for (color = elements->minColor; color <= elements->maxColor; color++)  
             {  
                 /*  open loop over all elements: */  
 #pragma omp for private(e) schedule(static)  
                 for (e = 0; e < elements->numElements; e++)  
                 {  
                 if (elements->Color[e] == color)  
                 {  
                     double vol = p.row_jac->absD[e] * p.row_jac->quadweight;  
                     D_p = getSampleDataRO(D, e);  
   
                     if (useHRZ) {   /* HRZ lumping */  
                        for (k = 0; k < p.numEqu; k++)  
                        {  
                       m_t = 0;  /* mass of the element: m_t */  
                       for (q = 0; q < p.numQuad; q++)  
                          m_t += vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];  
218    
219                        diagS = 0;    /* diagonal sum: S */                    for (color = elements->minColor; color <= elements->maxColor; color++)
220                        for (s = 0; s < p.numShapes; s++)                    {
221                        {                    /*  open loop over all elements: */
222                           rtmp = 0;        #pragma omp for private(e) schedule(static)
223                           for (q = 0; q < p.numQuad; q++)                    for (e = 0; e < elements->numElements; e++)
224                             rtmp +=                    {
225                                vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)] *                        if (elements->Color[e] == color)
226                                S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];                        {
227                           EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;                        double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
228                           diagS += rtmp;                        D_p = getSampleDataRO(D, e);
229                          if (useHRZ)   {   /* HRZ lumping */
230                            m_t = 0;    /* mass of the element: m_t */
231                            for (q = 0; q < p.numQuad; q++)
232                            m_t += vol;
233                            diagS = 0;  /* diagonal sum: S */
234                            for (s = 0; s < p.numShapes; s++)
235                            {
236                            rtmp = 0;
237                            for (q = 0; q < p.numQuad; q++)
238                            {
239                              rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
240                            }
241                            EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
242                            diagS += rtmp;
243                            }
244                            /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
245                            rtmp = m_t / diagS * D_p[0];
246                            for (s = 0; s < p.numShapes; s++)
247                            EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;
248                          } else {          /* row-sum lumping */
249                            for (s = 0; s < p.numShapes; s++)
250                            {
251                            rtmp = 0;
252                            for (q = 0; q < p.numQuad; q++)
253                              rtmp += vol * S[INDEX2(s, q, p.numShapes)];
254                            EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp * D_p[0];
255                            }
256                        }                        }
257                        /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                        for (q = 0; q < p.numShapes; q++)
258                        rtmp = m_t / diagS;                            row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
259                        for (s = 0; s < p.numShapes; s++)                        Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,
260                           EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp;                                  p.row_DOF_UpperBound);
261                         }                        } /* end color check */
262                      } else {                /* row-sum lumping */                    } /* end element loop */
263                         for (s = 0; s < p.numShapes; s++)                    } /* end color loop */
264                         {  
265                        for (k = 0; k < p.numEqu; k++)                }
266                        {                }
267                           rtmp = 0.;                else
268                           for (q = 0; q < p.numQuad; q++)                {     /* system of  equation */
269                             rtmp +=                if (expandedD)
270                                vol * S[INDEX2(s, q, p.numShapes)] *                {     /* with expanded D */
271                                D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];                    for (color = elements->minColor; color <= elements->maxColor; color++)
272                           EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;                    {
273                      /*  open loop over all elements: */
274          #pragma omp for private(e) schedule(static)
275                      for (e = 0; e < elements->numElements; e++)
276                      {
277                          if (elements->Color[e] == color)
278                          {
279                          double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
280                          D_p = getSampleDataRO(D, e);
281    
282                          if (useHRZ)   {   /* HRZ lumping */
283                            for (k = 0; k < p.numEqu; k++)
284                            {
285                            m_t = 0;    /* mass of the element: m_t */
286                            for (q = 0; q < p.numQuad; q++)
287                              m_t += vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];
288    
289                            diagS = 0;  /* diagonal sum: S */
290                            for (s = 0; s < p.numShapes; s++)
291                            {
292                              rtmp = 0;
293                              for (q = 0; q < p.numQuad; q++)
294                                rtmp +=
295                                    vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)] *
296                                    S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
297                              EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
298                              diagS += rtmp;
299                            }
300                            /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
301                            rtmp = m_t / diagS;
302                            for (s = 0; s < p.numShapes; s++)
303                              EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp;
304                            }
305                          } else {              /* row-sum lumping */
306                            for (s = 0; s < p.numShapes; s++)
307                            {
308                            for (k = 0; k < p.numEqu; k++)
309                            {
310                              rtmp = 0.;
311                              for (q = 0; q < p.numQuad; q++)
312                                rtmp +=
313                                    vol * S[INDEX2(s, q, p.numShapes)] *
314                                    D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];
315                              EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
316                            }
317                            }
318                        }                        }
319                         }                        for (q = 0; q < p.numShapes; q++)
320                      }                            row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
321                      for (q = 0; q < p.numShapes; q++)                        Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,
322                      row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];                                  p.row_DOF_UpperBound);
323                      Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,                        } /* end color check */
324                                 p.row_DOF_UpperBound);                    } /* end element loop */
325                  }   /* end color check */                    } /* end color loop */
326                  }   /* end element loop */                }
327              }   /* end color loop */                else
328              }                {     /* with constant D */
329              else                    for (color = elements->minColor; color <= elements->maxColor; color++)
330              {       /* with constant D */                    {
331              for (color = elements->minColor; color <= elements->maxColor; color++)                    /*  open loop over all elements: */
332              {        #pragma omp for private(e) schedule(static)
333                  /*  open loop over all elements: */                    for (e = 0; e < elements->numElements; e++)
334  #pragma omp for private(e) schedule(static)                    {
335                  for (e = 0; e < elements->numElements; e++)                        if (elements->Color[e] == color)
336                  {                        {
337                  if (elements->Color[e] == color)                        double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
338                  {                        D_p = getSampleDataRO(D, e);
339                      double vol = p.row_jac->absD[e] * p.row_jac->quadweight;  
340                      D_p = getSampleDataRO(D, e);                        if (useHRZ)       { /* HRZ lumping */
341                            m_t = 0;    /* mass of the element: m_t */
                     if (useHRZ)     { /* HRZ lumping */  
                        m_t = 0; /* mass of the element: m_t */  
                        for (q = 0; q < p.numQuad; q++)  
                       m_t += vol;  
                        diagS = 0;   /* diagonal sum: S */  
                        for (s = 0; s < p.numShapes; s++)  
                        {  
                       rtmp = 0;  
                       for (q = 0; q < p.numQuad; q++)  
                          rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];  
                       for (k = 0; k < p.numEqu; k++)  
                          EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;  
                       diagS += rtmp;  
                        }  
                        /* rescale diagonals by m_t/diagS to ensure consistent mass over element */  
                        rtmp = m_t / diagS;  
                        for (s = 0; s < p.numShapes; s++)  
                        {  
                       for (k = 0; k < p.numEqu; k++)  
                          EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp * D_p[k];  
                        }  
                     } else {                /* row-sum lumping */  
                        for (s = 0; s < p.numShapes; s++)  
                        {  
                       for (k = 0; k < p.numEqu; k++)  
                       {  
                         rtmp = 0.;  
342                          for (q = 0; q < p.numQuad; q++)                          for (q = 0; q < p.numQuad; q++)
343                          rtmp += vol * S[INDEX2(s, q, p.numShapes)];                          m_t += vol;
344                          EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp * D_p[k];                          diagS = 0;  /* diagonal sum: S */
345                            for (s = 0; s < p.numShapes; s++)
346                            {
347                            rtmp = 0;
348                            for (q = 0; q < p.numQuad; q++)
349                              rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
350                            for (k = 0; k < p.numEqu; k++)
351                              EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
352                            diagS += rtmp;
353                            }
354                            /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
355                            rtmp = m_t / diagS;
356                            for (s = 0; s < p.numShapes; s++)
357                            {
358                            for (k = 0; k < p.numEqu; k++)
359                              EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp * D_p[k];
360                            }
361                          } else {              /* row-sum lumping */
362                            for (s = 0; s < p.numShapes; s++)
363                            {
364                            for (k = 0; k < p.numEqu; k++)
365                            {
366                              rtmp = 0.;
367                              for (q = 0; q < p.numQuad; q++)
368                                  rtmp += vol * S[INDEX2(s, q, p.numShapes)];
369                              EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp * D_p[k];
370                            }
371                            }
372                        }                        }
373                         }                        for (q = 0; q < p.numShapes; q++)
374                      }                            row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
375                      for (q = 0; q < p.numShapes; q++)                        Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,
376                      row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];                                  p.row_DOF_UpperBound);
377                      Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,                        } /* end color check */
378                                 p.row_DOF_UpperBound);                    } /* end element loop */
379                  }   /* end color check */                    } /* end color loop */
380                  }   /* end element loop */                }
381              }   /* end color loop */                }
382              }            }         /* end of pointer check */
383          }            THREAD_MEMFREE(EM_lumpedMat);
384          }           /* end of pointer check */            THREAD_MEMFREE(row_index);
385          THREAD_MEMFREE(EM_lumpedMat);            }         /* end parallel region */
386          THREAD_MEMFREE(row_index);      }
     }           /* end parallel region */  
387      }      }
388  }  }

Legend:
Removed from v.3521  
changed lines
  Added in v.3522

  ViewVC Help
Powered by ViewVC 1.1.26