/[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 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,
33                                escript::Data* grad_data, const escript::Data* data)                                escript::Data* grad_data, const escript::Data* data)
34  {  {
35      size_t localGradSize = 0;      if (!nodes || !elements)
     dim_t e, q, l, s, n;  
     const double *data_array;  
     double *grad_data_e;  
     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;  
     reducedIntegrationOrder = Dudley_Assemble_reducedIntegrationOrder(grad_data);  
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 {
61      else          throw DudleyException("Assemble_gradient: Cannot calculate gradient "
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;
69          numShapesTotal = jac->numShapes;      const int numQuad = jac->numQuad;
70          numQuad = jac->numQuad;      const size_t localGradSize = sizeof(double) * numDim * numQuad * numComps;
71          localGradSize = sizeof(double) * numDim * numQuad * numComps;  
72          /* check the dimensions of data */      // check the dimensions of data
73        if (!grad_data->numSamplesEqual(numQuad, elements->numElements)) {
74          if (!numSamplesEqual(grad_data, numQuad, elements->numElements))          throw DudleyException("Assemble_gradient: illegal number of samples in gradient Data object");
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_Assemble_gradient: illegal number of components in gradient data object.");  
         }  
         else if (!isExpanded(grad_data))  
         {  
             Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: expanded Data object is expected for output data.");  
         }  
81      }      }
     /* now we can start */  
82    
83      if (Dudley_noError())      grad_data->requireWrite();
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++) {
91                  if (numDim == 1)                      double* grad_data_e = grad_data->getSampleDataRW(e);
92                  {                      memset(grad_data_e, 0, localGradSize);
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++) {
99                          memset(grad_data_e, 0, localGradSize);                                  grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
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++) {
110                                          data_array[l] *                      double* grad_data_e = grad_data->getSampleDataRW(e);
111                                          jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];                      memset(grad_data_e, 0, localGradSize);
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                      {                              }
125                          grad_data_e = getSampleDataRW(grad_data, e);                          }
126                          memset(grad_data_e, 0, localGradSize);                      }
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++) {
132                              {                      double* grad_data_e = grad_data->getSampleDataRW(e);
133  #pragma ivdep                      memset(grad_data_e, 0, localGradSize);
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++) {
140                                          data_array[l] *                                  grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
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)];
143                              }                                  grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
144                          }                                      data_array[l] *
145                      }                                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
146  #undef DIM                                  grad_data_e[INDEX4(l, 2, q, 0, numComps, numDim, numQuad)] +=
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  
 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)  
                     for (e = 0; e < elements->numElements; e++)  
                     {  
                         grad_data_e = getSampleDataRW(grad_data, e);  
                         memset(grad_data_e, 0, localGradSize);  
                         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++)  
                                 {  
                                     grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=  
                                         data_array[l] *  
                                         jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  
                                     grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=  
                                         data_array[l] *  
                                         jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  
                                     grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=  
                                         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++) {
159  #define DIM 1                      double* grad_data_e = grad_data->getSampleDataRW(e);
160  #pragma omp for schedule(static)                      memset(grad_data_e, 0, localGradSize);
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)];
163                          grad_data_e = getSampleDataRW(grad_data, e);                          const double* data_array = data->getSampleDataRO(nodes->reducedNodesMapping->target[n]);
164                          memset(grad_data_e, 0, localGradSize);                          for (int q = 0; q < numQuad; q++) {
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++) {
178                              }                      double* grad_data_e = grad_data->getSampleDataRW(e);
179                          }                      memset(grad_data_e, 0, localGradSize);
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)];
189                      {                                  grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
190                          grad_data_e = getSampleDataRW(grad_data, e);                                      data_array[l] *
191                          memset(grad_data_e, 0, localGradSize);                                      jac->DSDX[INDEX5(s, 1, 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++) {
200                                  {                      double* grad_data_e = grad_data->getSampleDataRW(e);
201                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                      memset(grad_data_e, 0, localGradSize);
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++) {
208                              }                                  grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
209                          }                                      data_array[l] *
210                      }                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
211  #undef DIM                                  grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
212                  }                                      data_array[l] *
213                  else if (numDim == 3)                                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
214                  {                                  grad_data_e[INDEX4(l, 2, q, 0, numComps, numDim, numQuad)] +=
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++)  
                     {  
                         grad_data_e = getSampleDataRW(grad_data, e);  
                         memset(grad_data_e, 0, localGradSize);  
                         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++)  
                                 {  
                                     grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=  
                                         data_array[l] *  
                                         jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  
                                     grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=  
                                         data_array[l] *  
                                         jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  
                                     grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=  
                                         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++) {
228  #define DIM 1                      double* grad_data_e = grad_data->getSampleDataRW(e);
229  #pragma omp for schedule(static)                      memset(grad_data_e, 0, localGradSize);
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)];
232                          grad_data_e = getSampleDataRW(grad_data, e);                          const double* data_array = data->getSampleDataRO(nodes->degreesOfFreedomMapping->target[n]);
233                          memset(grad_data_e, 0, localGradSize);                          for (int q = 0; q < numQuad; q++) {
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++) {
247                              }                      double* grad_data_e = grad_data->getSampleDataRW(e);
248                          }                      memset(grad_data_e, 0, localGradSize);
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)];
258                      {                                  grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
259                          grad_data_e = getSampleDataRW(grad_data, e);                                      data_array[l] *
260                          memset(grad_data_e, 0, localGradSize);                                      jac->DSDX[INDEX5(s, 1, 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++) {
269                                  {                      double* grad_data_e = grad_data->getSampleDataRW(e);
270                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                      memset(grad_data_e, 0, localGradSize);
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++) {
277                              }                                  grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
278                          }                                      data_array[l] *
279                      }                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
280  #undef DIM                                  grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
281                  }                                      data_array[l] *
282                  else if (numDim == 3)                                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
283                  {                                  grad_data_e[INDEX4(l, 2, q, 0, numComps, numDim, numQuad)] +=
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++)  
                     {  
                         grad_data_e = getSampleDataRW(grad_data, e);  
                         memset(grad_data_e, 0, localGradSize);  
                         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++)  
                                 {  
                                     grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=  
                                         data_array[l] *  
                                         jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  
                                     grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=  
                                         data_array[l] *  
                                         jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  
                                     grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=  
                                         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++) {
296  #define DIM 1                      double* grad_data_e = grad_data->getSampleDataRW(e);
297  #pragma omp for schedule(static)                      memset(grad_data_e, 0, localGradSize);
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)];
300                          grad_data_e = getSampleDataRW(grad_data, e);                          const double* data_array = data->getSampleDataRO(nodes->reducedDegreesOfFreedomMapping->target[n]);
301                          memset(grad_data_e, 0, localGradSize);                          for (int q = 0; q < numQuad; q++) {
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++) {
315                              }                      double* grad_data_e = grad_data->getSampleDataRW(e);
316                          }                      memset(grad_data_e, 0, localGradSize);
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)];
326                      {                                  grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
327                          grad_data_e = getSampleDataRW(grad_data, e);                                      data_array[l] *
328                          memset(grad_data_e, 0, localGradSize);                                      jac->DSDX[INDEX5(s, 1, 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++) {
337                                  {                      double* grad_data_e = grad_data->getSampleDataRW(e);
338                                      grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=                      memset(grad_data_e, 0, localGradSize);
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++) {
345                              }                                  grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
346                          }                                      data_array[l] *
347                      }                                      jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
348  #undef DIM                                  grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
349                  }                                      data_array[l] *
350                  else if (numDim == 3)                                      jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
351                  {                                  grad_data_e[INDEX4(l, 2, q, 0, numComps, numDim, numQuad)] +=
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++)  
                     {  
                         grad_data_e = getSampleDataRW(grad_data, e);  
                         memset(grad_data_e, 0, localGradSize);  
                         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++)  
                                 {  
                                     grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=  
                                         data_array[l] *  
                                         jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  
                                     grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=  
                                         data_array[l] *  
                                         jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];  
                                     grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=  
                                         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

  ViewVC Help
Powered by ViewVC 1.1.26