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

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

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

revision 5962 by caltinay, Fri Feb 5 03:37:49 2016 UTC revision 5963 by caltinay, Mon Feb 22 06:59:27 2016 UTC
# Line 34  Line 34 
34  /* 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 for loop unrolling */
35    
36  void Dudley_Assemble_gradient(Dudley_NodeFile * nodes, Dudley_ElementFile * elements,  void Dudley_Assemble_gradient(Dudley_NodeFile * nodes, Dudley_ElementFile * elements,
37                    escript::Data* grad_data, const escript::Data* data)                                escript::Data* grad_data, const escript::Data* data)
38  {  {
39      size_t localGradSize = 0;      size_t localGradSize = 0;
40      register dim_t e, q, l, s, n;      dim_t e, q, l, s, n;
41      register __const double *data_array;      const double *data_array;
42      register double *grad_data_e;      double *grad_data_e;
43      dim_t numNodes = 0, numShapesTotal = 0, numComps, NN = 0, numDim = 0, numQuad = 0;      dim_t numNodes = 0, numShapesTotal = 0, numComps, NN = 0, numDim = 0, numQuad = 0;
44      type_t data_type = getFunctionSpaceType(data);      int data_type = getFunctionSpaceType(data);
45      bool reducedIntegrationOrder = FALSE;      bool reducedIntegrationOrder = FALSE;
46      Dudley_ElementFile_Jacobeans *jac = NULL;      Dudley_ElementFile_Jacobeans *jac = NULL;
47    
48      Dudley_resetError();      Dudley_resetError();
49      if (nodes == NULL || elements == NULL)      if (nodes == NULL || elements == NULL)
50      return;          return;
51      numComps = getDataPointSize(data);      numComps = getDataPointSize(data);
52      NN = elements->numNodes;      NN = elements->numNodes;
53      reducedIntegrationOrder = Dudley_Assemble_reducedIntegrationOrder(grad_data);      reducedIntegrationOrder = Dudley_Assemble_reducedIntegrationOrder(grad_data);
54    
55      if (data_type == DUDLEY_NODES)      if (data_type == DUDLEY_NODES)
56      {      {
57      numNodes = nodes->nodesMapping->numTargets;          numNodes = nodes->nodesMapping->numTargets;
58      }      }
59      else if (data_type == DUDLEY_REDUCED_NODES)      else if (data_type == DUDLEY_REDUCED_NODES)
60      {      {
61      numNodes = nodes->reducedNodesMapping->numTargets;          numNodes = nodes->reducedNodesMapping->numTargets;
62      }      }
63      else if (data_type == DUDLEY_DEGREES_OF_FREEDOM)      else if (data_type == DUDLEY_DEGREES_OF_FREEDOM)
64      {      {
65      if (elements->MPIInfo->size > 1)          if (elements->MPIInfo->size > 1)
66      {          {
67          Dudley_setError(TYPE_ERROR,              Dudley_setError(TYPE_ERROR,
68                  "Dudley_Assemble_gradient: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");                              "Dudley_Assemble_gradient: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");
69          return;              return;
70      }          }
71      numNodes = nodes->degreesOfFreedomMapping->numTargets;          numNodes = nodes->degreesOfFreedomMapping->numTargets;
72      }      }
73      else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM)      else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM)
74      {      {
75      if (elements->MPIInfo->size > 1)          if (elements->MPIInfo->size > 1)
76      {          {
77          Dudley_setError(TYPE_ERROR,              Dudley_setError(TYPE_ERROR,
78                  "Dudley_Assemble_gradient: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");                              "Dudley_Assemble_gradient: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
79          return;              return;
80      }          }
81      numNodes = nodes->reducedDegreesOfFreedomMapping->numTargets;          numNodes = nodes->reducedDegreesOfFreedomMapping->numTargets;
82      }      }
83      else      else
84      {      {
85      Dudley_setError(TYPE_ERROR,          Dudley_setError(TYPE_ERROR,
86              "Dudley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");                          "Dudley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
87      }      }
88    
89      jac = Dudley_ElementFile_borrowJacobeans(elements, nodes, reducedIntegrationOrder);      jac = Dudley_ElementFile_borrowJacobeans(elements, nodes, reducedIntegrationOrder);
90      if (Dudley_noError())      if (Dudley_noError())
91      {      {
92      numDim = jac->numDim;          numDim = jac->numDim;
93      numShapesTotal = jac->numShapes;          numShapesTotal = jac->numShapes;
94      numQuad = jac->numQuad;          numQuad = jac->numQuad;
95      localGradSize = sizeof(double) * numDim * numQuad * numComps;          localGradSize = sizeof(double) * numDim * numQuad * numComps;
96      /* check the dimensions of data */          /* check the dimensions of data */
97    
98      if (!numSamplesEqual(grad_data, numQuad, elements->numElements))          if (!numSamplesEqual(grad_data, numQuad, elements->numElements))
99      {          {
100          Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: illegal number of samples in gradient Data object");              Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: illegal number of samples in gradient Data object");
101      }          }
102      else if (!numSamplesEqual(data, 1, numNodes))          else if (!numSamplesEqual(data, 1, numNodes))
103      {          {
104          Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: illegal number of samples of input Data object");              Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: illegal number of samples of input Data object");
105      }          }
106      else if (numDim * numComps != getDataPointSize(grad_data))          else if (numDim * numComps != getDataPointSize(grad_data))
107      {          {
108          Dudley_setError(TYPE_ERROR,              Dudley_setError(TYPE_ERROR,
109                  "Dudley_Assemble_gradient: illegal number of components in gradient data object.");                              "Dudley_Assemble_gradient: illegal number of components in gradient data object.");
110      }          }
111      else if (!isExpanded(grad_data))          else if (!isExpanded(grad_data))
112      {          {
113          Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: expanded Data object is expected for output data.");              Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: expanded Data object is expected for output data.");
114      }          }
115      }      }
116      /* now we can start */      /* now we can start */
117    
118      if (Dudley_noError())      if (Dudley_noError())
119      {      {
120      requireWrite(grad_data);          requireWrite(grad_data);
121  #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e)  #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e)
122      {          {
123          if (data_type == DUDLEY_NODES)              if (data_type == DUDLEY_NODES)
124          {              {
125          if (numDim == 1)                  if (numDim == 1)
126          {                  {
127              const dim_t numShapes = 2;                      const dim_t numShapes = 2;
128  #define DIM 1  #define DIM 1
129  #pragma omp for schedule(static)  #pragma omp for schedule(static)
130              for (e = 0; e < elements->numElements; e++)                      for (e = 0; e < elements->numElements; e++)
131              {                      {
132              grad_data_e = getSampleDataRW(grad_data, e);                          grad_data_e = getSampleDataRW(grad_data, e);
133              memset(grad_data_e, 0, localGradSize);                          memset(grad_data_e, 0, localGradSize);
134              for (s = 0; s < numShapes; s++)                          for (s = 0; s < numShapes; s++)
135              {                          {
136                  n = elements->Nodes[INDEX2(s, e, NN)];                              n = elements->Nodes[INDEX2(s, e, NN)];
137                  data_array = getSampleDataRO(data, n);                              data_array = getSampleDataRO(data, n);
138                  for (q = 0; q < numQuad; q++)                              for (q = 0; q < numQuad; q++)
139                  {                              {
140  #pragma ivdep  #pragma ivdep
141                  for (l = 0; l < numComps; l++)                                  for (l = 0; l < numComps; l++)
142                  {                                  {
143                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
144                      data_array[l] *                                          data_array[l] *
145                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
146                  }                                  }
147                  }                              }
148              }                          }
149              }                      }
150  #undef DIM  #undef DIM
151          }                  }
152          else if (numDim == 2)                  else if (numDim == 2)
153          {                  {
154              const dim_t numShapes = 3;                      const dim_t numShapes = 3;
155  #define DIM 2  #define DIM 2
156  #pragma omp for schedule(static)  #pragma omp for schedule(static)
157              for (e = 0; e < elements->numElements; e++)                      for (e = 0; e < elements->numElements; e++)
158              {                      {
159              grad_data_e = getSampleDataRW(grad_data, e);                          grad_data_e = getSampleDataRW(grad_data, e);
160              memset(grad_data_e, 0, localGradSize);                          memset(grad_data_e, 0, localGradSize);
161              for (s = 0; s < numShapes; s++)                          for (s = 0; s < numShapes; s++)
162              {                          {
163                  n = elements->Nodes[INDEX2(s, e, NN)];                              n = elements->Nodes[INDEX2(s, e, NN)];
164                  data_array = getSampleDataRO(data, n);                              data_array = getSampleDataRO(data, n);
165                  for (q = 0; q < numQuad; q++)                              for (q = 0; q < numQuad; q++)
166                  {                              {
167  #pragma ivdep  #pragma ivdep
168                  for (l = 0; l < numComps; l++)                                  for (l = 0; l < numComps; l++)
169                  {                                  {
170                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
171                      data_array[l] *                                          data_array[l] *
172                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
173                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
174                      data_array[l] *                                          data_array[l] *
175                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
176                  }                                  }
177                  }                              }
178              }                          }
179              }                      }
180  #undef DIM  #undef DIM
181          }                  }
182          else if (numDim == 3)                  else if (numDim == 3)
183          {                  {
184              const dim_t numShapes = 4;                      const dim_t numShapes = 4;
185  #define DIM 3  #define DIM 3
186  #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)  #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
187              for (e = 0; e < elements->numElements; e++)                      for (e = 0; e < elements->numElements; e++)
188              {                      {
189              grad_data_e = getSampleDataRW(grad_data, e);                          grad_data_e = getSampleDataRW(grad_data, e);
190              memset(grad_data_e, 0, localGradSize);                          memset(grad_data_e, 0, localGradSize);
191              for (s = 0; s < numShapes; s++)                          for (s = 0; s < numShapes; s++)
192              {                          {
193                  n = elements->Nodes[INDEX2(s, e, NN)];                              n = elements->Nodes[INDEX2(s, e, NN)];
194                  data_array = getSampleDataRO(data, n);                              data_array = getSampleDataRO(data, n);
195                  for (q = 0; q < numQuad; q++)                              for (q = 0; q < numQuad; q++)
196                  {                              {
197  #pragma ivdep  #pragma ivdep
198                  for (l = 0; l < numComps; l++)                                  for (l = 0; l < numComps; l++)
199                  {                                  {
200                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
201                      data_array[l] *                                          data_array[l] *
202                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
203                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
204                      data_array[l] *                                          data_array[l] *
205                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
206                      grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=
207                      data_array[l] *                                          data_array[l] *
208                      jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
209                  }                                  }
210                  }                              }
211              }                          }
212              }                      }
213  #undef DIM  #undef DIM
214          }                  }
215          }              }
216          else if (data_type == DUDLEY_REDUCED_NODES)              else if (data_type == DUDLEY_REDUCED_NODES)
217          {              {
218          if (numDim == 1)                  if (numDim == 1)
219          {                  {
220              const dim_t numShapes = 2;                      const dim_t numShapes = 2;
221  #define DIM 1  #define DIM 1
222  #pragma omp for schedule(static)  #pragma omp for schedule(static)
223              for (e = 0; e < elements->numElements; e++)                      for (e = 0; e < elements->numElements; e++)
224              {                      {
225              grad_data_e = getSampleDataRW(grad_data, e);                          grad_data_e = getSampleDataRW(grad_data, e);
226              memset(grad_data_e, 0, localGradSize);                          memset(grad_data_e, 0, localGradSize);
227              for (s = 0; s < numShapes; s++)                          for (s = 0; s < numShapes; s++)
228              {                          {
229                  n = elements->Nodes[INDEX2(s, e, NN)];                              n = elements->Nodes[INDEX2(s, e, NN)];
230                  data_array = getSampleDataRO(data, nodes->reducedNodesMapping->target[n]);                              data_array = getSampleDataRO(data, nodes->reducedNodesMapping->target[n]);
231                  for (q = 0; q < numQuad; q++)                              for (q = 0; q < numQuad; q++)
232                  {                              {
233  #pragma ivdep  #pragma ivdep
234                  for (l = 0; l < numComps; l++)                                  for (l = 0; l < numComps; l++)
235                  {                                  {
236                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
237                      data_array[l] *                                          data_array[l] *
238                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
239                  }                                  }
240                  }                              }
241              }                          }
242              }                      }
243  #undef DIM  #undef DIM
244          }                  }
245          else if (numDim == 2)                  else if (numDim == 2)
246          {                  {
247              const dim_t numShapes = 3;                      const dim_t numShapes = 3;
248  #define DIM 2  #define DIM 2
249  #pragma omp for schedule(static)  #pragma omp for schedule(static)
250              for (e = 0; e < elements->numElements; e++)                      for (e = 0; e < elements->numElements; e++)
251              {                      {
252              grad_data_e = getSampleDataRW(grad_data, e);                          grad_data_e = getSampleDataRW(grad_data, e);
253              memset(grad_data_e, 0, localGradSize);                          memset(grad_data_e, 0, localGradSize);
254              for (s = 0; s < numShapes; s++)                          for (s = 0; s < numShapes; s++)
255              {                          {
256                  n = elements->Nodes[INDEX2(s, e, NN)];                              n = elements->Nodes[INDEX2(s, e, NN)];
257                  data_array = getSampleDataRO(data, nodes->reducedNodesMapping->target[n]);                              data_array = getSampleDataRO(data, nodes->reducedNodesMapping->target[n]);
258                  for (q = 0; q < numQuad; q++)                              for (q = 0; q < numQuad; q++)
259                  {                              {
260  #pragma ivdep  #pragma ivdep
261                  for (l = 0; l < numComps; l++)                                  for (l = 0; l < numComps; l++)
262                  {                                  {
263                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
264                      data_array[l] *                                          data_array[l] *
265                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
266                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
267                      data_array[l] *                                          data_array[l] *
268                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
269                  }                                  }
270                  }                              }
271              }                          }
272              }                      }
273  #undef DIM  #undef DIM
274          }                  }
275          else if (numDim == 3)                  else if (numDim == 3)
276          {                  {
277              const dim_t numShapes = 4;                      const dim_t numShapes = 4;
278  #define DIM 3  #define DIM 3
279  #pragma omp for schedule(static)  #pragma omp for schedule(static)
280              for (e = 0; e < elements->numElements; e++)                      for (e = 0; e < elements->numElements; e++)
281              {                      {
282              grad_data_e = getSampleDataRW(grad_data, e);                          grad_data_e = getSampleDataRW(grad_data, e);
283              memset(grad_data_e, 0, localGradSize);                          memset(grad_data_e, 0, localGradSize);
284              for (s = 0; s < numShapes; s++)                          for (s = 0; s < numShapes; s++)
285              {                          {
286                  n = elements->Nodes[INDEX2(s, e, NN)];                              n = elements->Nodes[INDEX2(s, e, NN)];
287                  data_array = getSampleDataRO(data, nodes->reducedNodesMapping->target[n]);                              data_array = getSampleDataRO(data, nodes->reducedNodesMapping->target[n]);
288                  for (q = 0; q < numQuad; q++)                              for (q = 0; q < numQuad; q++)
289                  {                              {
290  #pragma ivdep  #pragma ivdep
291                  for (l = 0; l < numComps; l++)                                  for (l = 0; l < numComps; l++)
292                  {                                  {
293                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
294                      data_array[l] *                                          data_array[l] *
295                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
296                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
297                      data_array[l] *                                          data_array[l] *
298                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
299                      grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=
300                      data_array[l] *                                          data_array[l] *
301                      jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
302                  }                                  }
303                  }                              }
304              }                          }
305              }                      }
306  #undef DIM  #undef DIM
307          }                  }
308          }              }
309          else if (data_type == DUDLEY_DEGREES_OF_FREEDOM)              else if (data_type == DUDLEY_DEGREES_OF_FREEDOM)
310          {              {
311          if (numDim == 1)                  if (numDim == 1)
312          {                  {
313              const dim_t numShapes = 2;                      const dim_t numShapes = 2;
314  #define DIM 1  #define DIM 1
315  #pragma omp for schedule(static)  #pragma omp for schedule(static)
316              for (e = 0; e < elements->numElements; e++)                      for (e = 0; e < elements->numElements; e++)
317              {                      {
318              grad_data_e = getSampleDataRW(grad_data, e);                          grad_data_e = getSampleDataRW(grad_data, e);
319              memset(grad_data_e, 0, localGradSize);                          memset(grad_data_e, 0, localGradSize);
320              for (s = 0; s < numShapes; s++)                          for (s = 0; s < numShapes; s++)
321              {                          {
322                  n = elements->Nodes[INDEX2(s, e, NN)];                              n = elements->Nodes[INDEX2(s, e, NN)];
323                  data_array = getSampleDataRO(data, nodes->degreesOfFreedomMapping->target[n]);                              data_array = getSampleDataRO(data, nodes->degreesOfFreedomMapping->target[n]);
324                  for (q = 0; q < numQuad; q++)                              for (q = 0; q < numQuad; q++)
325                  {                              {
326  #pragma ivdep  #pragma ivdep
327                  for (l = 0; l < numComps; l++)                                  for (l = 0; l < numComps; l++)
328                  {                                  {
329                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
330                      data_array[l] *                                          data_array[l] *
331                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
332                  }                                  }
333                  }                              }
334              }                          }
335              }                      }
336  #undef DIM  #undef DIM
337          }                  }
338          else if (numDim == 2)                  else if (numDim == 2)
339          {                  {
340              const dim_t numShapes = 3;                      const dim_t numShapes = 3;
341  #define DIM 2  #define DIM 2
342  #pragma omp for schedule(static)  #pragma omp for schedule(static)
343              for (e = 0; e < elements->numElements; e++)                      for (e = 0; e < elements->numElements; e++)
344              {                      {
345              grad_data_e = getSampleDataRW(grad_data, e);                          grad_data_e = getSampleDataRW(grad_data, e);
346              memset(grad_data_e, 0, localGradSize);                          memset(grad_data_e, 0, localGradSize);
347              for (s = 0; s < numShapes; s++)                          for (s = 0; s < numShapes; s++)
348              {                          {
349                  n = elements->Nodes[INDEX2(s, e, NN)];                              n = elements->Nodes[INDEX2(s, e, NN)];
350                  data_array = getSampleDataRO(data, nodes->degreesOfFreedomMapping->target[n]);                              data_array = getSampleDataRO(data, nodes->degreesOfFreedomMapping->target[n]);
351                  for (q = 0; q < numQuad; q++)                              for (q = 0; q < numQuad; q++)
352                  {                              {
353  #pragma ivdep  #pragma ivdep
354                  for (l = 0; l < numComps; l++)                                  for (l = 0; l < numComps; l++)
355                  {                                  {
356                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
357                      data_array[l] *                                          data_array[l] *
358                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
359                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
360                      data_array[l] *                                          data_array[l] *
361                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
362                  }                                  }
363                  }                              }
364              }                          }
365              }                      }
366  #undef DIM  #undef DIM
367          }                  }
368          else if (numDim == 3)                  else if (numDim == 3)
369          {                  {
370              const dim_t numShapes = 4;                      const dim_t numShapes = 4;
371  #define DIM 3  #define DIM 3
372  #pragma omp for schedule(static)  #pragma omp for schedule(static)
373              for (e = 0; e < elements->numElements; e++)                      for (e = 0; e < elements->numElements; e++)
374              {                      {
375              grad_data_e = getSampleDataRW(grad_data, e);                          grad_data_e = getSampleDataRW(grad_data, e);
376              memset(grad_data_e, 0, localGradSize);                          memset(grad_data_e, 0, localGradSize);
377              for (s = 0; s < numShapes; s++)                          for (s = 0; s < numShapes; s++)
378              {                          {
379                  n = elements->Nodes[INDEX2(s, e, NN)];                              n = elements->Nodes[INDEX2(s, e, NN)];
380                  data_array = getSampleDataRO(data, nodes->degreesOfFreedomMapping->target[n]);                              data_array = getSampleDataRO(data, nodes->degreesOfFreedomMapping->target[n]);
381                  for (q = 0; q < numQuad; q++)                              for (q = 0; q < numQuad; q++)
382                  {                              {
383  #pragma ivdep  #pragma ivdep
384                  for (l = 0; l < numComps; l++)                                  for (l = 0; l < numComps; l++)
385                  {                                  {
386                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
387                      data_array[l] *                                          data_array[l] *
388                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
389                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
390                      data_array[l] *                                          data_array[l] *
391                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
392                      grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=
393                      data_array[l] *                                          data_array[l] *
394                      jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
395                  }                                  }
396                  }                              }
397              }                          }
398              }                      }
399  #undef DIM  #undef DIM
400          }                  }
401          }              }
402          else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM)              else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM)
403          {              {
404          if (numDim == 1)                  if (numDim == 1)
405          {                  {
406              const dim_t numShapes = 2;                      const dim_t numShapes = 2;
407  #define DIM 1  #define DIM 1
408  #pragma omp for schedule(static)  #pragma omp for schedule(static)
409              for (e = 0; e < elements->numElements; e++)                      for (e = 0; e < elements->numElements; e++)
410              {                      {
411              grad_data_e = getSampleDataRW(grad_data, e);                          grad_data_e = getSampleDataRW(grad_data, e);
412              memset(grad_data_e, 0, localGradSize);                          memset(grad_data_e, 0, localGradSize);
413              for (s = 0; s < numShapes; s++)                          for (s = 0; s < numShapes; s++)
414              {                          {
415                  n = elements->Nodes[INDEX2(s, e, NN)];                              n = elements->Nodes[INDEX2(s, e, NN)];
416                  data_array = getSampleDataRO(data, nodes->reducedDegreesOfFreedomMapping->target[n]);                              data_array = getSampleDataRO(data, nodes->reducedDegreesOfFreedomMapping->target[n]);
417                  for (q = 0; q < numQuad; q++)                              for (q = 0; q < numQuad; q++)
418                  {                              {
419  #pragma ivdep  #pragma ivdep
420                  for (l = 0; l < numComps; l++)                                  for (l = 0; l < numComps; l++)
421                  {                                  {
422                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
423                      data_array[l] *                                          data_array[l] *
424                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
425                  }                                  }
426                  }                              }
427              }                          }
428              }                      }
429  #undef DIM  #undef DIM
430          }                  }
431          else if (numDim == 2)                  else if (numDim == 2)
432          {                  {
433              const dim_t numShapes = 3;                      const dim_t numShapes = 3;
434  #define DIM 2  #define DIM 2
435  #pragma omp for schedule(static)  #pragma omp for schedule(static)
436              for (e = 0; e < elements->numElements; e++)                      for (e = 0; e < elements->numElements; e++)
437              {                      {
438              grad_data_e = getSampleDataRW(grad_data, e);                          grad_data_e = getSampleDataRW(grad_data, e);
439              memset(grad_data_e, 0, localGradSize);                          memset(grad_data_e, 0, localGradSize);
440              for (s = 0; s < numShapes; s++)                          for (s = 0; s < numShapes; s++)
441              {                          {
442                  n = elements->Nodes[INDEX2(s, e, NN)];                              n = elements->Nodes[INDEX2(s, e, NN)];
443                  data_array = getSampleDataRO(data, nodes->reducedDegreesOfFreedomMapping->target[n]);                              data_array = getSampleDataRO(data, nodes->reducedDegreesOfFreedomMapping->target[n]);
444                  for (q = 0; q < numQuad; q++)                              for (q = 0; q < numQuad; q++)
445                  {                              {
446  #pragma ivdep  #pragma ivdep
447                  for (l = 0; l < numComps; l++)                                  for (l = 0; l < numComps; l++)
448                  {                                  {
449                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
450                      data_array[l] *                                          data_array[l] *
451                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
452                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
453                      data_array[l] *                                          data_array[l] *
454                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
455                  }                                  }
456                  }                              }
457              }                          }
458              }                      }
459  #undef DIM  #undef DIM
460          }                  }
461          else if (numDim == 3)                  else if (numDim == 3)
462          {                  {
463              const dim_t numShapes = 4;                      const dim_t numShapes = 4;
464  #define DIM 3  #define DIM 3
465  #pragma omp for schedule(static)  #pragma omp for schedule(static)
466              for (e = 0; e < elements->numElements; e++)                      for (e = 0; e < elements->numElements; e++)
467              {                      {
468              grad_data_e = getSampleDataRW(grad_data, e);                          grad_data_e = getSampleDataRW(grad_data, e);
469              memset(grad_data_e, 0, localGradSize);                          memset(grad_data_e, 0, localGradSize);
470              for (s = 0; s < numShapes; s++)                          for (s = 0; s < numShapes; s++)
471              {                          {
472                  n = elements->Nodes[INDEX2(s, e, NN)];                              n = elements->Nodes[INDEX2(s, e, NN)];
473                  data_array = getSampleDataRO(data, nodes->reducedDegreesOfFreedomMapping->target[n]);                              data_array = getSampleDataRO(data, nodes->reducedDegreesOfFreedomMapping->target[n]);
474                  for (q = 0; q < numQuad; q++)                              for (q = 0; q < numQuad; q++)
475                  {                              {
476  #pragma ivdep  #pragma ivdep
477                  for (l = 0; l < numComps; l++)                                  for (l = 0; l < numComps; l++)
478                  {                                  {
479                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
480                      data_array[l] *                                          data_array[l] *
481                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
482                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
483                      data_array[l] *                                          data_array[l] *
484                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
485                      grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=                                      grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=
486                      data_array[l] *                                          data_array[l] *
487                      jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                                          jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
488                  }                                  }
489                  }                              }
490              }                          }
491              }                      }
492  #undef DIM  #undef DIM
493          }                  }
494          }              }
495      }           /* end parallel region */          }                       /* end parallel region */
496      }      }
497  }  }

Legend:
Removed from v.5962  
changed lines
  Added in v.5963

  ViewVC Help
Powered by ViewVC 1.1.26