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  /*    assemblage of Jacobians: calculate the gradient of nodal data at quadrature points */    Assemblage of Jacobians: calculate the gradient of nodal data at quadrature
20      points
21
22  /************************************************************************************/  *****************************************************************************/

#define ESNEEDPYTHON
#include "esysUtils/first.h"
23
24  #include "Assemble.h"  #include "Assemble.h"
25  #include "Util.h"  #include "Util.h"
#ifdef _OPENMP
#include <omp.h>
#endif
/***************************************************************************************/
26
27  /* Unless the loops in here get complicated again this file should be compiled for loop unrolling */  // Unless the loops in here get complicated again this file should be compiled
28    // with loop unrolling
29
30    namespace dudley {
31
32  void Dudley_Assemble_gradient(Dudley_NodeFile * nodes, Dudley_ElementFile * elements,  void Assemble_gradient(Dudley_NodeFile* nodes, Dudley_ElementFile* elements,
34  {  {
35      size_t localGradSize = 0;      if (!nodes || !elements)
dim_t e, q, l, s, n;
const double *data_array;
dim_t numNodes = 0, numShapesTotal = 0, numComps, NN = 0, numDim = 0, numQuad = 0;
int data_type = getFunctionSpaceType(data);
bool reducedIntegrationOrder = FALSE;
Dudley_ElementFile_Jacobeans *jac = NULL;

Dudley_resetError();
if (nodes == NULL || elements == NULL)
36          return;          return;
numComps = getDataPointSize(data);
NN = elements->numNodes;
37
38      if (data_type == DUDLEY_NODES)      const int numComps = data->getDataPointSize();
39      {      const int NN = elements->numNodes;
40        const bool reducedIntegrationOrder = Assemble_reducedIntegrationOrder(grad_data);
41        const int data_type = data->getFunctionSpace().getTypeCode();
42
43        dim_t numNodes = 0;
44        if (data_type == DUDLEY_NODES) {
45          numNodes = nodes->nodesMapping->numTargets;          numNodes = nodes->nodesMapping->numTargets;
46      }      } else if (data_type == DUDLEY_REDUCED_NODES) {
else if (data_type == DUDLEY_REDUCED_NODES)
{
47          numNodes = nodes->reducedNodesMapping->numTargets;          numNodes = nodes->reducedNodesMapping->numTargets;
48      }      } else if (data_type == DUDLEY_DEGREES_OF_FREEDOM) {
49      else if (data_type == DUDLEY_DEGREES_OF_FREEDOM)          if (elements->MPIInfo->size > 1) {
50      {              throw DudleyException("Assemble_gradient: for more than one "
51          if (elements->MPIInfo->size > 1)                  "processor DEGREES_OF_FREEDOM data are not accepted as input.");
{
Dudley_setError(TYPE_ERROR,
"Dudley_Assemble_gradient: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");
return;
52          }          }
53          numNodes = nodes->degreesOfFreedomMapping->numTargets;          numNodes = nodes->degreesOfFreedomMapping->numTargets;
54      }      } else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {
55      else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM)          if (elements->MPIInfo->size > 1) {
56      {              throw DudleyException("Assemble_gradient: for more than one "
57          if (elements->MPIInfo->size > 1)                  "processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
{
Dudley_setError(TYPE_ERROR,
"Dudley_Assemble_gradient: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
return;
58          }          }
59          numNodes = nodes->reducedDegreesOfFreedomMapping->numTargets;          numNodes = nodes->reducedDegreesOfFreedomMapping->numTargets;
60      }      } else {
62      {                 "of data because of unsuitable input data representation.");
Dudley_setError(TYPE_ERROR,
"Dudley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
63      }      }
64
65      jac = Dudley_ElementFile_borrowJacobeans(elements, nodes, reducedIntegrationOrder);      Dudley_ElementFile_Jacobians* jac = Dudley_ElementFile_borrowJacobians(
66      if (Dudley_noError())              elements, nodes, reducedIntegrationOrder);
67      {      const int numDim = jac->numDim;
68          numDim = jac->numDim;      const int numShapesTotal = jac->numShapes;
72          /* check the dimensions of data */      // check the dimensions of data
75          {      } else if (!data->numSamplesEqual(1, numNodes)) {
76              Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: illegal number of samples in gradient Data object");          throw DudleyException("Assemble_gradient: illegal number of samples of input Data object");
77          }      } else if (numDim * numComps != grad_data->getDataPointSize()) {
78          else if (!numSamplesEqual(data, 1, numNodes))          throw DudleyException("Assemble_gradient: illegal number of components in gradient data object.");
79          {      } else if (!grad_data->actsExpanded()) {
80              Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: illegal number of samples of input Data object");          throw DudleyException("Assemble_gradient: expanded Data object is expected for output data.");
}
else if (numDim * numComps != getDataPointSize(grad_data))
{
Dudley_setError(TYPE_ERROR,
}
{
Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: expanded Data object is expected for output data.");
}
81      }      }
/* now we can start */
82
84    #pragma omp parallel
85      {      {
86          requireWrite(grad_data);          if (data_type == DUDLEY_NODES) {
87  #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e)              if (numDim == 1) {
88          {                  const dim_t numShapes = 2;
89              if (data_type == DUDLEY_NODES)  #pragma omp for
90              {                  for (index_t e = 0; e < elements->numElements; e++) {
93                      const dim_t numShapes = 2;                      for (int s = 0; s < numShapes; s++) {
94  #define DIM 1                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
95  #pragma omp for schedule(static)                          const double* data_array = data->getSampleDataRO(n);
96                      for (e = 0; e < elements->numElements; e++)                          for (int q = 0; q < numQuad; q++) {
97                      {  #pragma ivdep
98                          grad_data_e = getSampleDataRW(grad_data, e);                              for (int l = 0; l < numComps; l++) {
100                          for (s = 0; s < numShapes; s++)                                      data_array[l] *
101                          {                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
102                              n = elements->Nodes[INDEX2(s, e, NN)];                              }
103                              data_array = getSampleDataRO(data, n);                          }
104                              for (q = 0; q < numQuad; q++)                      }
105                              {                  }
106  #pragma ivdep              } else if (numDim == 2) {
107                                  for (l = 0; l < numComps; l++)                  const dim_t numShapes = 3;
108                                  {  #pragma omp for
109                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                  for (index_t e = 0; e < elements->numElements; e++) {
112                                  }                      for (int s = 0; s < numShapes; s++) {
113                              }                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
114                          }                          const double* data_array = data->getSampleDataRO(n);
115                      }                          for (int q = 0; q < numQuad; q++) {
116  #undef DIM  #pragma ivdep
117                  }                              for (int l = 0; l < numComps; l++) {
118                  else if (numDim == 2)                                  grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
119                  {                                      data_array[l] *
120                      const dim_t numShapes = 3;                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
121  #define DIM 2                                  grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
122  #pragma omp for schedule(static)                                      data_array[l] *
123                      for (e = 0; e < elements->numElements; e++)                                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
124                      {                              }
127                          for (s = 0; s < numShapes; s++)                  }
128                          {              } else if (numDim == 3) {
129                              n = elements->Nodes[INDEX2(s, e, NN)];                  const dim_t numShapes = 4;
130                              data_array = getSampleDataRO(data, n);  #pragma omp for
131                              for (q = 0; q < numQuad; q++)                  for (index_t e = 0; e < elements->numElements; e++) {
134                                  for (l = 0; l < numComps; l++)                      for (int s = 0; s < numShapes; s++) {
135                                  {                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
136                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                          const double* data_array = data->getSampleDataRO(n);
137                                          data_array[l] *                          for (int q = 0; q < numQuad; q++) {
138                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  #pragma ivdep
139                                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=                              for (int l = 0; l < numComps; l++) {
141                                          jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                      data_array[l] *
142                                  }                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
144                          }                                      data_array[l] *
145                      }                                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
147                  }                                      data_array[l] *
148                  else if (numDim == 3)                                      jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
{
const dim_t numShapes = 4;
#define DIM 3
for (e = 0; e < elements->numElements; e++)
{
for (s = 0; s < numShapes; s++)
{
n = elements->Nodes[INDEX2(s, e, NN)];
data_array = getSampleDataRO(data, n);
for (q = 0; q < numQuad; q++)
{
#pragma ivdep
for (l = 0; l < numComps; l++)
{
data_array[l] *
jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
data_array[l] *
jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
data_array[l] *
jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
}
149                              }                              }
150                          }                          }
151                      }                      }
#undef DIM
152                  }                  }
153              }              }
154              else if (data_type == DUDLEY_REDUCED_NODES)          } else if (data_type == DUDLEY_REDUCED_NODES) {
155              {              if (numDim == 1) {
156                  if (numDim == 1)                  const dim_t numShapes = 2;
157                  {  #pragma omp for
158                      const dim_t numShapes = 2;                  for (index_t e = 0; e < elements->numElements; e++) {
161                      for (e = 0; e < elements->numElements; e++)                      for (int s = 0; s < numShapes; s++) {
162                      {                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
165                          for (s = 0; s < numShapes; s++)  #pragma ivdep
166                          {                              for (int l = 0; l < numComps; l++) {
167                              n = elements->Nodes[INDEX2(s, e, NN)];                                  grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
168                              data_array = getSampleDataRO(data, nodes->reducedNodesMapping->target[n]);                                      data_array[l] *
169                              for (q = 0; q < numQuad; q++)                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
170                              {                              }
171  #pragma ivdep                          }
172                                  for (l = 0; l < numComps; l++)                      }
173                                  {                  }
174                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=              } else if (numDim == 2) {
175                                          data_array[l] *                  const dim_t numShapes = 3;
176                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  #pragma omp for
177                                  }                  for (index_t e = 0; e < elements->numElements; e++) {
180                      }                      for (int s = 0; s < numShapes; s++) {
181  #undef DIM                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
182                  }                          const double* data_array = data->getSampleDataRO(nodes->reducedNodesMapping->target[n]);
183                  else if (numDim == 2)                          for (int q = 0; q < numQuad; q++) {
184                  {  #pragma ivdep
185                      const dim_t numShapes = 3;                              for (int l = 0; l < numComps; l++) {
186  #define DIM 2                                  grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
187  #pragma omp for schedule(static)                                      data_array[l] *
188                      for (e = 0; e < elements->numElements; e++)                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
192                          for (s = 0; s < numShapes; s++)                              }
193                          {                          }
194                              n = elements->Nodes[INDEX2(s, e, NN)];                      }
195                              data_array = getSampleDataRO(data, nodes->reducedNodesMapping->target[n]);                  }
196                              for (q = 0; q < numQuad; q++)              } else if (numDim == 3) {
197                              {                  const dim_t numShapes = 4;
198  #pragma ivdep  #pragma omp for
199                                  for (l = 0; l < numComps; l++)                  for (index_t e = 0; e < elements->numElements; e++) {
202                                          data_array[l] *                      for (int s = 0; s < numShapes; s++) {
203                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
204                                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=                          const double* data_array = data->getSampleDataRO(nodes->reducedNodesMapping->target[n]);
205                                          data_array[l] *                          for (int q = 0; q < numQuad; q++) {
206                                          jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  #pragma ivdep
207                                  }                              for (int l = 0; l < numComps; l++) {
209                          }                                      data_array[l] *
210                      }                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
212                  }                                      data_array[l] *
213                  else if (numDim == 3)                                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
215                      const dim_t numShapes = 4;                                      data_array[l] *
216  #define DIM 3                                      jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
#pragma omp for schedule(static)
for (e = 0; e < elements->numElements; e++)
{
for (s = 0; s < numShapes; s++)
{
n = elements->Nodes[INDEX2(s, e, NN)];
data_array = getSampleDataRO(data, nodes->reducedNodesMapping->target[n]);
for (q = 0; q < numQuad; q++)
{
#pragma ivdep
for (l = 0; l < numComps; l++)
{
data_array[l] *
jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
data_array[l] *
jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
data_array[l] *
jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
}
217                              }                              }
218                          }                          }
219                      }                      }
#undef DIM
220                  }                  }
221              }              }
222              else if (data_type == DUDLEY_DEGREES_OF_FREEDOM)          } else if (data_type == DUDLEY_DEGREES_OF_FREEDOM) {
223                if (numDim == 1)
224              {              {
225                  if (numDim == 1)                  const dim_t numShapes = 2;
226                  {  #pragma omp for
227                      const dim_t numShapes = 2;                  for (index_t e = 0; e < elements->numElements; e++) {
230                      for (e = 0; e < elements->numElements; e++)                      for (int s = 0; s < numShapes; s++) {
231                      {                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
234                          for (s = 0; s < numShapes; s++)  #pragma ivdep
235                          {                              for (int l = 0; l < numComps; l++) {
236                              n = elements->Nodes[INDEX2(s, e, NN)];                                  grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
237                              data_array = getSampleDataRO(data, nodes->degreesOfFreedomMapping->target[n]);                                      data_array[l] *
238                              for (q = 0; q < numQuad; q++)                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
239                              {                              }
240  #pragma ivdep                          }
241                                  for (l = 0; l < numComps; l++)                      }
242                                  {                  }
243                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=              } else if (numDim == 2) {
244                                          data_array[l] *                  const dim_t numShapes = 3;
245                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  #pragma omp for
246                                  }                  for (index_t e = 0; e < elements->numElements; e++) {
249                      }                      for (int s = 0; s < numShapes; s++) {
250  #undef DIM                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
251                  }                          const double* data_array = data->getSampleDataRO(nodes->degreesOfFreedomMapping->target[n]);
252                  else if (numDim == 2)                          for (int q = 0; q < numQuad; q++) {
253                  {  #pragma ivdep
254                      const dim_t numShapes = 3;                              for (int l = 0; l < numComps; l++) {
255  #define DIM 2                                  grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
256  #pragma omp for schedule(static)                                      data_array[l] *
257                      for (e = 0; e < elements->numElements; e++)                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
261                          for (s = 0; s < numShapes; s++)                              }
262                          {                          }
263                              n = elements->Nodes[INDEX2(s, e, NN)];                      }
264                              data_array = getSampleDataRO(data, nodes->degreesOfFreedomMapping->target[n]);                  }
265                              for (q = 0; q < numQuad; q++)              } else if (numDim == 3) {
266                              {                  const dim_t numShapes = 4;
267  #pragma ivdep  #pragma omp for
268                                  for (l = 0; l < numComps; l++)                  for (index_t e = 0; e < elements->numElements; e++) {
271                                          data_array[l] *                      for (int s = 0; s < numShapes; s++) {
272                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
273                                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=                          const double* data_array = data->getSampleDataRO(nodes->degreesOfFreedomMapping->target[n]);
274                                          data_array[l] *                          for (int q = 0; q < numQuad; q++) {
275                                          jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  #pragma ivdep
276                                  }                              for (int l = 0; l < numComps; l++) {
278                          }                                      data_array[l] *
279                      }                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
281                  }                                      data_array[l] *
282                  else if (numDim == 3)                                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
284                      const dim_t numShapes = 4;                                      data_array[l] *
285  #define DIM 3                                      jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
#pragma omp for schedule(static)
for (e = 0; e < elements->numElements; e++)
{
for (s = 0; s < numShapes; s++)
{
n = elements->Nodes[INDEX2(s, e, NN)];
data_array = getSampleDataRO(data, nodes->degreesOfFreedomMapping->target[n]);
for (q = 0; q < numQuad; q++)
{
#pragma ivdep
for (l = 0; l < numComps; l++)
{
data_array[l] *
jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
data_array[l] *
jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
data_array[l] *
jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
}
286                              }                              }
287                          }                          }
288                      }                      }
#undef DIM
289                  }                  }
290              }              }
291              else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM)          } else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {
292              {              if (numDim == 1) {
293                  if (numDim == 1)                  const dim_t numShapes = 2;
294                  {  #pragma omp for
295                      const dim_t numShapes = 2;                  for (index_t e = 0; e < elements->numElements; e++) {
298                      for (e = 0; e < elements->numElements; e++)                      for (int s = 0; s < numShapes; s++) {
299                      {                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
302                          for (s = 0; s < numShapes; s++)  #pragma ivdep
303                          {                              for (int l = 0; l < numComps; l++) {
304                              n = elements->Nodes[INDEX2(s, e, NN)];                                  grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
305                              data_array = getSampleDataRO(data, nodes->reducedDegreesOfFreedomMapping->target[n]);                                      data_array[l] *
306                              for (q = 0; q < numQuad; q++)                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
307                              {                              }
308  #pragma ivdep                          }
309                                  for (l = 0; l < numComps; l++)                      }
310                                  {                  }
311                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=              } else if (numDim == 2) {
312                                          data_array[l] *                  const dim_t numShapes = 3;
313                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  #pragma omp for
314                                  }                  for (index_t e = 0; e < elements->numElements; e++) {
317                      }                      for (int s = 0; s < numShapes; s++) {
318  #undef DIM                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
319                  }                          const double* data_array = data->getSampleDataRO(nodes->reducedDegreesOfFreedomMapping->target[n]);
320                  else if (numDim == 2)                          for (int q = 0; q < numQuad; q++) {
321                  {  #pragma ivdep
322                      const dim_t numShapes = 3;                              for (int l = 0; l < numComps; l++) {
323  #define DIM 2                                  grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
324  #pragma omp for schedule(static)                                      data_array[l] *
325                      for (e = 0; e < elements->numElements; e++)                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
329                          for (s = 0; s < numShapes; s++)                              }
330                          {                          }
331                              n = elements->Nodes[INDEX2(s, e, NN)];                      }
332                              data_array = getSampleDataRO(data, nodes->reducedDegreesOfFreedomMapping->target[n]);                  }
333                              for (q = 0; q < numQuad; q++)              } else if (numDim == 3) {
334                              {                  const dim_t numShapes = 4;
335  #pragma ivdep  #pragma omp for
336                                  for (l = 0; l < numComps; l++)                  for (index_t e = 0; e < elements->numElements; e++) {
339                                          data_array[l] *                      for (int s = 0; s < numShapes; s++) {
340                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
341                                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=                          const double* data_array = data->getSampleDataRO(nodes->reducedDegreesOfFreedomMapping->target[n]);
342                                          data_array[l] *                          for (int q = 0; q < numQuad; q++) {
343                                          jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  #pragma ivdep
344                                  }                              for (int l = 0; l < numComps; l++) {
346                          }                                      data_array[l] *
347                      }                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
349                  }                                      data_array[l] *
350                  else if (numDim == 3)                                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
352                      const dim_t numShapes = 4;                                      data_array[l] *
353  #define DIM 3                                      jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
#pragma omp for schedule(static)
for (e = 0; e < elements->numElements; e++)
{
for (s = 0; s < numShapes; s++)
{
n = elements->Nodes[INDEX2(s, e, NN)];
data_array = getSampleDataRO(data, nodes->reducedDegreesOfFreedomMapping->target[n]);
for (q = 0; q < numQuad; q++)
{
#pragma ivdep
for (l = 0; l < numComps; l++)
{
data_array[l] *
jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
data_array[l] *
jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
data_array[l] *
jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
}
354                              }                              }
355                          }                          }
356                      }                      }
#undef DIM
357                  }                  }
358              }              }
359          }                       /* end parallel region */          }
360      }      } // end parallel region
361  }  }
362
363    } // namespace dudley
364

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