/[escript]/branches/trilinos_from_5897/dudley/src/Assemble_LumpedSystem.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/dudley/src/Assemble_LumpedSystem.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 6008 by caltinay, Mon Feb 22 06:59:27 2016 UTC revision 6009 by caltinay, Wed Mar 2 04:13:26 2016 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
17  /************************************************************************************/  /****************************************************************************
18    
19  /*    assembles the mass matrix in lumped form                */    Assembles the mass matrix in lumped form.
20    
21  /*    The coefficient D has to be defined on the integration points or not present. */    The coefficient D has to be defined on the integration points or not present.
22      lumpedMat has to be initialized before the routine is called.
23    
24  /*    lumpedMat has to be initialized before the routine is called. */  *****************************************************************************/
   
 /************************************************************************************/  
   
 #define ESNEEDPYTHON  
 #include "esysUtils/first.h"  
25    
26  #include "Assemble.h"  #include "Assemble.h"
 #include "Util.h"  
 #ifdef _OPENMP  
 #include <omp.h>  
 #endif  
   
27  #include "ShapeTable.h"  #include "ShapeTable.h"
28    #include "Util.h"
29    
30  /************************************************************************************/  namespace dudley {
31    
32  void Dudley_Assemble_LumpedSystem(Dudley_NodeFile * nodes, Dudley_ElementFile * elements, escript::Data * lumpedMat,  void Assemble_LumpedSystem(Dudley_NodeFile* nodes, Dudley_ElementFile* elements,
33                                    const escript::Data * D, const bool useHRZ)                             escript::Data& lumpedMat, const escript::Data& D,
34                               bool useHRZ)
35  {  {
36        if (!nodes || !elements || lumpedMat.isEmpty() || D.isEmpty())
     bool reducedIntegrationOrder = FALSE, expandedD;  
     Dudley_Assemble_Parameters p;  
     int dimensions[ESCRIPT_MAX_DATA_RANK];  
     dim_t k, e, len_EM_lumpedMat, q, s;  
     int funcspace;  
     index_t color, *row_index = NULL;  
     __const double *D_p = NULL;  
     const double *S = NULL;  
     double *EM_lumpedMat = NULL, *lumpedMat_p = NULL;  
     double rtmp;  
     double m_t = 0., diagS = 0.;  
   
     Dudley_resetError();  
   
     if (nodes == NULL || elements == NULL)  
         return;  
     if (lumpedMat->isEmpty() || D->isEmpty())  
37          return;          return;
38      if (lumpedMat->isEmpty() && !D->isEmpty())  
39      {      const int funcspace = D.getFunctionSpace().getTypeCode();
40          Dudley_setError(TYPE_ERROR, "Dudley_Assemble_LumpedSystem: coefficients are non-zero but no lumped matrix is given.");      bool reducedIntegrationOrder;
41          return;      // check function space of D
42      }      if (funcspace == DUDLEY_ELEMENTS) {
43      funcspace = D->getFunctionSpace().getTypeCode();          reducedIntegrationOrder = false;
44      /* check if all function spaces are the same */      } else if (funcspace == DUDLEY_FACE_ELEMENTS) {
45      if (funcspace == DUDLEY_ELEMENTS)          reducedIntegrationOrder = false;
46      {      } else if (funcspace == DUDLEY_REDUCED_ELEMENTS) {
47          reducedIntegrationOrder = FALSE;          reducedIntegrationOrder = true;
48      }      } else if (funcspace == DUDLEY_REDUCED_FACE_ELEMENTS) {
49      else if (funcspace == DUDLEY_FACE_ELEMENTS)          reducedIntegrationOrder = true;
50      {      } else {
51          reducedIntegrationOrder = FALSE;          throw DudleyException("Assemble_LumpedSystem: assemblage failed because of illegal function space.");
     }  
     else if (funcspace == DUDLEY_REDUCED_ELEMENTS)  
     {  
         reducedIntegrationOrder = TRUE;  
     }  
     else if (funcspace == DUDLEY_REDUCED_FACE_ELEMENTS)  
     {  
         reducedIntegrationOrder = TRUE;  
     }  
     else  
     {  
         Dudley_setError(TYPE_ERROR, "Dudley_Assemble_LumpedSystem: assemblage failed because of illegal function space.");  
52      }      }
     if (!Dudley_noError())  
         return;  
53    
54      /* set all parameters in p */      // initialize parameters
55      Dudley_Assemble_getAssembleParameters(nodes, elements, escript::ASM_ptr(), lumpedMat, reducedIntegrationOrder, &p);      Assemble_Parameters p;
56      if (!Dudley_noError())      Assemble_getAssembleParameters(nodes, elements, escript::ASM_ptr(),
57          return;                                        lumpedMat, reducedIntegrationOrder, &p);
58    
59      /* check if all function spaces are the same */      // check if all function spaces are the same
60      if (!numSamplesEqual(D, p.numQuad, elements->numElements))      if (!D.numSamplesEqual(p.numQuad, elements->numElements)) {
     {  
61          std::stringstream ss;          std::stringstream ss;
62          ss << "Dudley_Assemble_LumpedSystem: sample points of coefficient D "          ss << "Assemble_LumpedSystem: sample points of coefficient D "
63                "don't match (" << p.numQuad << ","                "don't match (" << p.numQuad << ","
64             << elements->numElements << ")";             << elements->numElements << ")";
65          std::string error_msg(ss.str());          const std::string msg(ss.str());
66          Dudley_setError(TYPE_ERROR, error_msg.c_str());          throw DudleyException(msg);
67      }      }
68    
69      /*  check the dimensions: */      // check the dimensions
70      if (p.numEqu == 1)      if (p.numEqu == 1) {
71      {          const escript::DataTypes::ShapeType dimensions; //dummy
72          if (!D->isEmpty())          if (D.getDataPointShape() != dimensions) {
73          {              throw DudleyException("Assemble_LumpedSystem: coefficient D, rank 0 expected.");
             if (!isDataPointShapeEqual(D, 0, dimensions))  
             {  
                 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_LumpedSystem: coefficient D, rank 0 expected.");  
             }  
   
74          }          }
75      }      } else {
76      else          const escript::DataTypes::ShapeType dimensions(1, p.numEqu);
77      {          if (D.getDataPointShape() != dimensions) {
78          if (!D->isEmpty())              std::stringstream ss;
79          {              ss << "Assemble_LumpedSystem: coefficient D, expected "
80              dimensions[0] = p.numEqu;                    "shape (" << p.numEqu << ",)";
81              if (!isDataPointShapeEqual(D, 1, dimensions))              const std::string msg(ss.str());
82              {              throw DudleyException(msg);
                 std::stringstream ss;  
                 ss << "Dudley_Assemble_LumpedSystem: coefficient D, expected "  
                       "shape (" << dimensions[0] << ",)";  
                 std::string error_msg(ss.str());  
                 Dudley_setError(TYPE_ERROR, error_msg.c_str());  
             }  
83          }          }
84      }      }
85      if (Dudley_noError())  
86      {      lumpedMat.requireWrite();
87          requireWrite(lumpedMat);      double* lumpedMat_p = lumpedMat.getSampleDataRW(0);
88          lumpedMat_p = getSampleDataRW(lumpedMat, 0);      
89                if (funcspace==DUDLEY_POINTS) {
90          if (funcspace==DUDLEY_POINTS) {  #pragma omp parallel
91                #pragma omp parallel private(color, D_p)          {
92                {              for (int color=elements->minColor; color<=elements->maxColor; color++) {
93                      for (color=elements->minColor;color<=elements->maxColor;color++) {                  // loop over all elements
94                        /*  open loop over all elements: */  #pragma omp for
95                        #pragma omp for private(e) schedule(static)                  for (index_t e=0; e<elements->numElements; e++) {
96                        for(e=0;e<elements->numElements;e++){                      if (elements->Color[e]==color) {
97                            if (elements->Color[e]==color) {                          const double* D_p = D.getSampleDataRO(e);
98                              D_p=getSampleDataRO(D, e);                          Dudley_Util_AddScatter(1,
99                              if (NULL!=D_p)  Dudley_Util_AddScatter(1,                                        &p.row_DOF[elements->Nodes[INDEX2(0,e,p.NN)]],
100                                                          &(p.row_DOF[elements->Nodes[INDEX2(0,e,p.NN)]]),                                        p.numEqu, D_p, lumpedMat_p,
101                                                          p.numEqu,                                        p.row_DOF_UpperBound);
102                                                          D_p,                      } // end color check
103                                                          lumpedMat_p,                  } // end element loop
104                                                          p.row_DOF_UpperBound);              } // end color loop
105                            } /* end color check */          } // end parallel region
106                        } /* end element loop */      } else {
107                    } /* end color loop */          bool expandedD = D.actsExpanded();
108              } /* end parallel region */          const double *S = NULL;
109          } else {            if (!getQuadShape(elements->numDim, reducedIntegrationOrder, &S)) {
110                              throw DudleyException("Assemble_LumpedSystem: Unable to locate shape function.");
               len_EM_lumpedMat = p.numShapes * p.numEqu;  
   
               expandedD = D->isExpanded();  
               if (!getQuadShape(elements->numDim, reducedIntegrationOrder, &S))  
               {  
                   Dudley_setError(TYPE_ERROR, "Dudley_Assemble_LumpedSystem: Unable to locate shape function.");  
               }  
               #pragma omp parallel private(color, EM_lumpedMat, row_index, D_p, s, q, k, rtmp, diagS, m_t)  
               {  
                   EM_lumpedMat = new double[len_EM_lumpedMat];  
                   row_index = new index_t[p.numShapes];  
                   if (!Dudley_checkPtr(EM_lumpedMat) && !Dudley_checkPtr(row_index))  
                   {  
                       if (p.numEqu == 1)  
                       {         /* single 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)   {  
                                             m_t = 0;    /* mass of the element: m_t */  
                                             for (q = 0; q < p.numQuad; q++)  
                                                 m_t += vol * D_p[INDEX2(q, 0, p.numQuad)];  
                                             diagS = 0;  /* diagonal sum: S */  
                                             for (s = 0; s < p.numShapes; s++)  
                                             {  
                                                 rtmp = 0;  
                                                 for (q = 0; q < p.numQuad; q++)  
                                                   rtmp +=  
                                                         vol * D_p[INDEX2(q, 0, p.numQuad)] * S[INDEX2(s, q, p.numShapes)] *  
                                                         S[INDEX2(s, q, p.numShapes)];  
                                                 EM_lumpedMat[INDEX2(0, 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++)  
                                                 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)] * D_p[INDEX2(q, 0, p.numQuad)];  
                                                 EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;  
                                             }  
                                           }  
                                           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++)  
                                                 {  
                                                   rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];  
                                                 }  
                                                 EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;  
                                                 diagS += rtmp;  
                                             }  
                                             /* 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];  
                                             }  
                                           }  
                                           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)];  
   
                                                 diagS = 0;      /* diagonal sum: S */  
                                                 for (s = 0; s < p.numShapes; s++)  
                                                 {  
                                                   rtmp = 0;  
                                                   for (q = 0; q < p.numQuad; q++)  
                                                         rtmp +=  
                                                             vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)] *  
                                                             S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];  
                                                   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++)  
                                                   EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp;  
                                             }  
                                           } else {                              /* row-sum lumping */  
                                             for (s = 0; s < p.numShapes; s++)  
                                             {  
                                                 for (k = 0; k < p.numEqu; k++)  
                                                 {  
                                                   rtmp = 0.;  
                                                   for (q = 0; q < p.numQuad; q++)  
                                                         rtmp +=  
                                                             vol * S[INDEX2(s, q, p.numShapes)] *  
                                                             D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];  
                                                   EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;  
                                                 }  
                                             }  
                                           }  
                                           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++)  
                                                   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.;  
                                                   for (q = 0; q < p.numQuad; q++)  
                                                       rtmp += vol * S[INDEX2(s, q, p.numShapes)];  
                                                   EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp * D_p[k];  
                                                 }  
                                             }  
                                           }  
                                           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 */  
                           }  
                       }  
                   }                     /* end of pointer check */  
                   delete[] EM_lumpedMat;  
                   delete[] row_index;  
               }                 /* end parallel region */  
111          }          }
112    #pragma omp parallel
113            {
114                std::vector<double> EM_lumpedMat(p.numShapes * p.numEqu);
115                std::vector<index_t> row_index(p.numShapes);
116    
117                if (p.numEqu == 1) { // single equation
118                    if (expandedD) { // with expanded D
119                        for (int color = elements->minColor; color <= elements->maxColor; color++) {
120                            // loop over all elements
121    #pragma omp for
122                            for (index_t e = 0; e < elements->numElements; e++) {
123                                if (elements->Color[e] == color) {
124                                    const double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
125                                    const double* D_p = D.getSampleDataRO(e);
126                                    if (useHRZ) {
127                                        double m_t = 0; // mass of the element
128                                        for (int q = 0; q < p.numQuad; q++)
129                                            m_t += vol * D_p[INDEX2(q, 0, p.numQuad)];
130                                        double diagS = 0;  // diagonal sum
131                                        double rtmp;
132                                        for (int s = 0; s < p.numShapes; s++) {
133                                            rtmp = 0.;
134                                            for (int q = 0; q < p.numQuad; q++)
135                                                rtmp +=
136                                                    vol * D_p[INDEX2(q, 0, p.numQuad)] * S[INDEX2(s, q, p.numShapes)] *
137                                                    S[INDEX2(s, q, p.numShapes)];
138                                            EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
139                                            diagS += rtmp;
140                                        }
141                                        // rescale diagonals by m_t/diagS to ensure
142                                        // consistent mass over element
143                                        rtmp = m_t / diagS;
144                                        for (int s = 0; s < p.numShapes; s++)
145                                            EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;
146                                    } else { // row-sum lumping
147                                        for (int s = 0; s < p.numShapes; s++) {
148                                            double rtmp = 0.;
149                                            for (int q = 0; q < p.numQuad; q++)
150                                                rtmp += vol * S[INDEX2(s, q, p.numShapes)] * D_p[INDEX2(q, 0, p.numQuad)];
151                                            EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
152                                        }
153                                    }
154                                    for (int q = 0; q < p.numShapes; q++)
155                                        row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
156                                    Dudley_Util_AddScatter(p.numShapes, &row_index[0],
157                                           p.numEqu, &EM_lumpedMat[0], lumpedMat_p,
158                                           p.row_DOF_UpperBound);
159                                } // end color check
160                            } // end element loop
161                        } // end color loop
162                    } else { // with constant D
163                        for (int color = elements->minColor; color <= elements->maxColor; color++) {
164                            // loop over all elements
165    #pragma omp for
166                            for (index_t e = 0; e < elements->numElements; e++) {
167                                if (elements->Color[e] == color) {
168                                    const double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
169                                    const double* D_p = D.getSampleDataRO(e);
170                                    if (useHRZ) { // HRZ lumping
171                                        // mass of the element
172                                        const double m_t = vol*p.numQuad;
173                                        double diagS = 0; // diagonal sum
174                                        double rtmp;
175                                        for (int s = 0; s < p.numShapes; s++) {
176                                            rtmp = 0.;
177                                            for (int q = 0; q < p.numQuad; q++) {
178                                                rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
179                                            }
180                                            EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
181                                            diagS += rtmp;
182                                        }
183                                        // rescale diagonals by m_t/diagS to ensure
184                                        // consistent mass over element
185                                        rtmp = m_t / diagS * D_p[0];
186                                        for (int s = 0; s < p.numShapes; s++)
187                                            EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;
188                                    } else { // row-sum lumping
189                                        for (int s = 0; s < p.numShapes; s++) {
190                                            double rtmp = 0.;
191                                            for (int q = 0; q < p.numQuad; q++)
192                                                rtmp += vol * S[INDEX2(s, q, p.numShapes)];
193                                            EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp * D_p[0];
194                                        }
195                                    }
196                                    for (int q = 0; q < p.numShapes; q++)
197                                        row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
198                                    Dudley_Util_AddScatter(p.numShapes, &row_index[0],
199                                           p.numEqu, &EM_lumpedMat[0], lumpedMat_p,
200                                           p.row_DOF_UpperBound);
201                                } // end color check
202                            } // end element loop
203                        } // end color loop
204                    }
205    
206                } else { // system of equations
207                    if (expandedD) { // with expanded D
208                        for (int color = elements->minColor; color <= elements->maxColor; color++) {
209                            // loop over all elements
210    #pragma omp for
211                            for (index_t e = 0; e < elements->numElements; e++) {
212                                if (elements->Color[e] == color) {
213                                    const double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
214                                    const double* D_p = D.getSampleDataRO(e);
215    
216                                    if (useHRZ) { // HRZ lumping
217                                        for (int k = 0; k < p.numEqu; k++) {
218                                            double m_t = 0; // mass of the element
219                                            for (int q = 0; q < p.numQuad; q++)
220                                                m_t += vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];
221    
222                                            double diagS = 0; // diagonal sum
223                                            double rtmp;
224                                            for (int s = 0; s < p.numShapes; s++) {
225                                                rtmp = 0.;
226                                                for (int q = 0; q < p.numQuad; q++)
227                                                    rtmp +=
228                                                        vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)] *
229                                                        S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
230                                                EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
231                                                diagS += rtmp;
232                                            }
233                                            // rescale diagonals by m_t/diagS to
234                                            // ensure consistent mass over element
235                                            rtmp = m_t / diagS;
236                                            for (int s = 0; s < p.numShapes; s++)
237                                                EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp;
238                                        }
239                                    } else { // row-sum lumping
240                                        for (int s = 0; s < p.numShapes; s++) {
241                                            for (int k = 0; k < p.numEqu; k++) {
242                                                double rtmp = 0.;
243                                                for (int q = 0; q < p.numQuad; q++)
244                                                    rtmp +=
245                                                        vol * S[INDEX2(s, q, p.numShapes)] *
246                                                        D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];
247                                                EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
248                                            }
249                                        }
250                                    }
251                                    for (int q = 0; q < p.numShapes; q++)
252                                        row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
253                                    Dudley_Util_AddScatter(p.numShapes, &row_index[0],
254                                           p.numEqu, &EM_lumpedMat[0], lumpedMat_p,
255                                           p.row_DOF_UpperBound);
256                                } // end color check
257                            } // end element loop
258                        } // end color loop
259                    } else { // with constant D
260                        for (int color = elements->minColor; color <= elements->maxColor; color++) {
261                            // loop over all elements
262    #pragma omp for
263                            for (index_t e = 0; e < elements->numElements; e++) {
264                                if (elements->Color[e] == color) {
265                                    const double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
266                                    const double* D_p = D.getSampleDataRO(e);
267    
268                                    if (useHRZ) { // HRZ lumping
269                                        double m_t = vol * p.numQuad; // mass of the element
270                                        double diagS = 0; // diagonal sum
271                                        double rtmp;
272                                        for (int s = 0; s < p.numShapes; s++) {
273                                            rtmp = 0.;
274                                            for (int q = 0; q < p.numQuad; q++)
275                                                rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
276                                            for (int k = 0; k < p.numEqu; k++)
277                                                EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
278                                            diagS += rtmp;
279                                        }
280                                        // rescale diagonals by m_t/diagS to ensure
281                                        // consistent mass over element
282                                        rtmp = m_t / diagS;
283                                        for (int s = 0; s < p.numShapes; s++) {
284                                            for (int k = 0; k < p.numEqu; k++)
285                                                EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp * D_p[k];
286                                        }
287                                    } else { // row-sum lumping
288                                        for (int s = 0; s < p.numShapes; s++) {
289                                            for (int k = 0; k < p.numEqu; k++) {
290                                                double rtmp = 0.;
291                                                for (int q = 0; q < p.numQuad; q++)
292                                                    rtmp += vol * S[INDEX2(s, q, p.numShapes)];
293                                                EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp * D_p[k];
294                                            }
295                                        }
296                                    }
297                                    for (int q = 0; q < p.numShapes; q++)
298                                        row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
299                                    Dudley_Util_AddScatter(p.numShapes, &row_index[0],
300                                           p.numEqu, &EM_lumpedMat[0], lumpedMat_p,
301                                           p.row_DOF_UpperBound);
302                                } // end color check
303                            } // end element loop
304                        } // end color loop
305                    }
306                }
307            } // end parallel region
308      }      }
309  }  }
310    
311    } // namespace dudley
312    

Legend:
Removed from v.6008  
changed lines
  Added in v.6009

  ViewVC Help
Powered by ViewVC 1.1.26