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

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