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

revision 6078 by caltinay, Wed Mar 2 04:13:26 2016 UTC revision 6079 by caltinay, Mon Mar 21 12:22:38 2016 UTC
# Line 14  Line 14
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16
/****************************************************************************

Assemblage of Jacobians: calculate the gradient of nodal data at quadrature
points

*****************************************************************************/

17  #include "Assemble.h"  #include "Assemble.h"
18  #include "Util.h"  #include "Util.h"
19
# Line 29  Line 22
22
23  namespace dudley {  namespace dudley {
24
25  void Assemble_gradient(Dudley_NodeFile* nodes, Dudley_ElementFile* elements,  void Assemble_gradient(const NodeFile* nodes, const ElementFile* elements,
26                                escript::Data* grad_data, const escript::Data* data)                         escript::Data& grad_data, const escript::Data& data)
27  {  {
28      if (!nodes || !elements)      if (!nodes || !elements)
29          return;          return;
30
31      const int numComps = data->getDataPointSize();      const int numComps = data.getDataPointSize();
32      const int NN = elements->numNodes;      const int NN = elements->numNodes;
33      const bool reducedIntegrationOrder = Assemble_reducedIntegrationOrder(grad_data);      const bool reducedIntegrationOrder = hasReducedIntegrationOrder(grad_data);
34      const int data_type = data->getFunctionSpace().getTypeCode();      const int data_type = data.getFunctionSpace().getTypeCode();
35
36      dim_t numNodes = 0;      dim_t numNodes = 0;
37      if (data_type == DUDLEY_NODES) {      if (data_type == DUDLEY_NODES) {
38          numNodes = nodes->nodesMapping->numTargets;          numNodes = nodes->getNumNodes();
} else if (data_type == DUDLEY_REDUCED_NODES) {
numNodes = nodes->reducedNodesMapping->numTargets;
39      } else if (data_type == DUDLEY_DEGREES_OF_FREEDOM) {      } else if (data_type == DUDLEY_DEGREES_OF_FREEDOM) {
40          if (elements->MPIInfo->size > 1) {          if (elements->MPIInfo->size > 1) {
41              throw DudleyException("Assemble_gradient: for more than one "              throw DudleyException("Assemble_gradient: for more than one "
42                  "processor DEGREES_OF_FREEDOM data are not accepted as input.");                  "processor DEGREES_OF_FREEDOM data are not accepted as input.");
43          }          }
44          numNodes = nodes->degreesOfFreedomMapping->numTargets;          numNodes = nodes->getNumDegreesOfFreedom();
} else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {
if (elements->MPIInfo->size > 1) {
throw DudleyException("Assemble_gradient: for more than one "
"processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
}
numNodes = nodes->reducedDegreesOfFreedomMapping->numTargets;
45      } else {      } else {
46          throw DudleyException("Assemble_gradient: Cannot calculate gradient "          throw DudleyException("Assemble_gradient: Cannot calculate gradient "
47                 "of data because of unsuitable input data representation.");                 "of data because of unsuitable input data representation.");
48      }      }
49
50      Dudley_ElementFile_Jacobians* jac = Dudley_ElementFile_borrowJacobians(      ElementFile_Jacobians* jac = elements->borrowJacobians(nodes,
51              elements, nodes, reducedIntegrationOrder);                                                       reducedIntegrationOrder);
52      const int numDim = jac->numDim;      const int numDim = jac->numDim;
53      const int numShapesTotal = jac->numShapes;      const int numShapesTotal = jac->numShapes;
54      const int numQuad = jac->numQuad;      const int numQuad = jac->numQuad;
55      const size_t localGradSize = sizeof(double) * numDim * numQuad * numComps;      const size_t localGradSize = sizeof(double) * numDim * numQuad * numComps;
56
57      // check the dimensions of data      // check the dimensions of data
58      if (!grad_data->numSamplesEqual(numQuad, elements->numElements)) {      if (!grad_data.numSamplesEqual(numQuad, elements->numElements)) {
59          throw DudleyException("Assemble_gradient: illegal number of samples in gradient Data object");          throw DudleyException("Assemble_gradient: illegal number of samples in gradient Data object");
60      } else if (!data->numSamplesEqual(1, numNodes)) {      } else if (!data.numSamplesEqual(1, numNodes)) {
61          throw DudleyException("Assemble_gradient: illegal number of samples of input Data object");          throw DudleyException("Assemble_gradient: illegal number of samples of input Data object");
62      } else if (numDim * numComps != grad_data->getDataPointSize()) {      } else if (numDim * numComps != grad_data.getDataPointSize()) {
63          throw DudleyException("Assemble_gradient: illegal number of components in gradient data object.");          throw DudleyException("Assemble_gradient: illegal number of components in gradient data object.");
64      } else if (!grad_data->actsExpanded()) {      } else if (!grad_data.actsExpanded()) {
65          throw DudleyException("Assemble_gradient: expanded Data object is expected for output data.");          throw DudleyException("Assemble_gradient: expanded Data object is expected for output data.");
66      }      }
67
68      grad_data->requireWrite();      grad_data.requireWrite();
69  #pragma omp parallel  #pragma omp parallel
70      {      {
71          if (data_type == DUDLEY_NODES) {          if (data_type == DUDLEY_NODES) {
72              if (numDim == 1) {              if (numDim == 1) {
73                  const dim_t numShapes = 2;                  const int numShapes = 2;
#pragma omp for
for (index_t e = 0; e < elements->numElements; e++) {
double* grad_data_e = grad_data->getSampleDataRW(e);
memset(grad_data_e, 0, localGradSize);
for (int s = 0; s < numShapes; s++) {
const index_t n = elements->Nodes[INDEX2(s, e, NN)];
const double* data_array = data->getSampleDataRO(n);
for (int q = 0; q < numQuad; q++) {
#pragma ivdep
for (int l = 0; l < numComps; l++) {
grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
data_array[l] *
jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
}
}
}
}
} else if (numDim == 2) {
const dim_t numShapes = 3;
#pragma omp for
for (index_t e = 0; e < elements->numElements; e++) {
double* grad_data_e = grad_data->getSampleDataRW(e);
memset(grad_data_e, 0, localGradSize);
for (int s = 0; s < numShapes; s++) {
const index_t n = elements->Nodes[INDEX2(s, e, NN)];
const double* data_array = data->getSampleDataRO(n);
for (int q = 0; q < numQuad; q++) {
#pragma ivdep
for (int l = 0; l < numComps; l++) {
grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
data_array[l] *
jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
data_array[l] *
jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
}
}
}
}
} else if (numDim == 3) {
const dim_t numShapes = 4;
#pragma omp for
for (index_t e = 0; e < elements->numElements; e++) {
double* grad_data_e = grad_data->getSampleDataRW(e);
memset(grad_data_e, 0, localGradSize);
for (int s = 0; s < numShapes; s++) {
const index_t n = elements->Nodes[INDEX2(s, e, NN)];
const double* data_array = data->getSampleDataRO(n);
for (int q = 0; q < numQuad; q++) {
#pragma ivdep
for (int l = 0; l < numComps; l++) {
grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
data_array[l] *
jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
data_array[l] *
jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
grad_data_e[INDEX4(l, 2, q, 0, numComps, numDim, numQuad)] +=
data_array[l] *
jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
}
}
}
}
}
} else if (data_type == DUDLEY_REDUCED_NODES) {
if (numDim == 1) {
const dim_t numShapes = 2;
74  #pragma omp for  #pragma omp for
75                  for (index_t e = 0; e < elements->numElements; e++) {                  for (index_t e = 0; e < elements->numElements; e++) {
76                      double* grad_data_e = grad_data->getSampleDataRW(e);                      double* grad_data_e = grad_data.getSampleDataRW(e);
77                      memset(grad_data_e, 0, localGradSize);                      memset(grad_data_e, 0, localGradSize);
78                      for (int s = 0; s < numShapes; s++) {                      for (int s = 0; s < numShapes; s++) {
79                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
80                          const double* data_array = data->getSampleDataRO(nodes->reducedNodesMapping->target[n]);                          const double* data_array = data.getSampleDataRO(n);
81                          for (int q = 0; q < numQuad; q++) {                          for (int q = 0; q < numQuad; q++) {
82  #pragma ivdep  #pragma ivdep
83                              for (int l = 0; l < numComps; l++) {                              for (int l = 0; l < numComps; l++) {
# Line 172  void Assemble_gradient(Dudley_NodeFile* Line 89  void Assemble_gradient(Dudley_NodeFile*
89                      }                      }
90                  }                  }
91              } else if (numDim == 2) {              } else if (numDim == 2) {
92                  const dim_t numShapes = 3;                  const int numShapes = 3;
93  #pragma omp for  #pragma omp for
94                  for (index_t e = 0; e < elements->numElements; e++) {                  for (index_t e = 0; e < elements->numElements; e++) {
95                      double* grad_data_e = grad_data->getSampleDataRW(e);                      double* grad_data_e = grad_data.getSampleDataRW(e);
96                      memset(grad_data_e, 0, localGradSize);                      memset(grad_data_e, 0, localGradSize);
97                      for (int s = 0; s < numShapes; s++) {                      for (int s = 0; s < numShapes; s++) {
98                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
99                          const double* data_array = data->getSampleDataRO(nodes->reducedNodesMapping->target[n]);                          const double* data_array = data.getSampleDataRO(n);
100                          for (int q = 0; q < numQuad; q++) {                          for (int q = 0; q < numQuad; q++) {
101  #pragma ivdep  #pragma ivdep
102                              for (int l = 0; l < numComps; l++) {                              for (int l = 0; l < numComps; l++) {
# Line 194  void Assemble_gradient(Dudley_NodeFile* Line 111  void Assemble_gradient(Dudley_NodeFile*
111                      }                      }
112                  }                  }
113              } else if (numDim == 3) {              } else if (numDim == 3) {
114                  const dim_t numShapes = 4;                  const int numShapes = 4;
115  #pragma omp for  #pragma omp for
116                  for (index_t e = 0; e < elements->numElements; e++) {                  for (index_t e = 0; e < elements->numElements; e++) {
117                      double* grad_data_e = grad_data->getSampleDataRW(e);                      double* grad_data_e = grad_data.getSampleDataRW(e);
118                      memset(grad_data_e, 0, localGradSize);                      memset(grad_data_e, 0, localGradSize);
119                      for (int s = 0; s < numShapes; s++) {                      for (int s = 0; s < numShapes; s++) {
120                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
121                          const double* data_array = data->getSampleDataRO(nodes->reducedNodesMapping->target[n]);                          const double* data_array = data.getSampleDataRO(n);
122                          for (int q = 0; q < numQuad; q++) {                          for (int q = 0; q < numQuad; q++) {
123  #pragma ivdep  #pragma ivdep
124                              for (int l = 0; l < numComps; l++) {                              for (int l = 0; l < numComps; l++) {
# Line 220  void Assemble_gradient(Dudley_NodeFile* Line 137  void Assemble_gradient(Dudley_NodeFile*
137                  }                  }
138              }              }
139          } else if (data_type == DUDLEY_DEGREES_OF_FREEDOM) {          } else if (data_type == DUDLEY_DEGREES_OF_FREEDOM) {
140              if (numDim == 1)              const index_t* target = nodes->borrowTargetDegreesOfFreedom();
{
const dim_t numShapes = 2;
#pragma omp for
for (index_t e = 0; e < elements->numElements; e++) {
double* grad_data_e = grad_data->getSampleDataRW(e);
memset(grad_data_e, 0, localGradSize);
for (int s = 0; s < numShapes; s++) {
const index_t n = elements->Nodes[INDEX2(s, e, NN)];
const double* data_array = data->getSampleDataRO(nodes->degreesOfFreedomMapping->target[n]);
for (int q = 0; q < numQuad; q++) {
#pragma ivdep
for (int l = 0; l < numComps; l++) {
grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
data_array[l] *
jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
}
}
}
}
} else if (numDim == 2) {
const dim_t numShapes = 3;
#pragma omp for
for (index_t e = 0; e < elements->numElements; e++) {
double* grad_data_e = grad_data->getSampleDataRW(e);
memset(grad_data_e, 0, localGradSize);
for (int s = 0; s < numShapes; s++) {
const index_t n = elements->Nodes[INDEX2(s, e, NN)];
const double* data_array = data->getSampleDataRO(nodes->degreesOfFreedomMapping->target[n]);
for (int q = 0; q < numQuad; q++) {
#pragma ivdep
for (int l = 0; l < numComps; l++) {
grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
data_array[l] *
jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
data_array[l] *
jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
}
}
}
}
} else if (numDim == 3) {
const dim_t numShapes = 4;
#pragma omp for
for (index_t e = 0; e < elements->numElements; e++) {
double* grad_data_e = grad_data->getSampleDataRW(e);
memset(grad_data_e, 0, localGradSize);
for (int s = 0; s < numShapes; s++) {
const index_t n = elements->Nodes[INDEX2(s, e, NN)];
const double* data_array = data->getSampleDataRO(nodes->degreesOfFreedomMapping->target[n]);
for (int q = 0; q < numQuad; q++) {
#pragma ivdep
for (int l = 0; l < numComps; l++) {
grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
data_array[l] *
jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
data_array[l] *
jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
grad_data_e[INDEX4(l, 2, q, 0, numComps, numDim, numQuad)] +=
data_array[l] *
jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
}
}
}
}
}
} else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {
141              if (numDim == 1) {              if (numDim == 1) {
142                  const dim_t numShapes = 2;                  const int numShapes = 2;
143  #pragma omp for  #pragma omp for
144                  for (index_t e = 0; e < elements->numElements; e++) {                  for (index_t e = 0; e < elements->numElements; e++) {
145                      double* grad_data_e = grad_data->getSampleDataRW(e);                      double* grad_data_e = grad_data.getSampleDataRW(e);
146                      memset(grad_data_e, 0, localGradSize);                      memset(grad_data_e, 0, localGradSize);
147                      for (int s = 0; s < numShapes; s++) {                      for (int s = 0; s < numShapes; s++) {
148                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
149                          const double* data_array = data->getSampleDataRO(nodes->reducedDegreesOfFreedomMapping->target[n]);                          const double* data_array = data.getSampleDataRO(target[n]);
150                          for (int q = 0; q < numQuad; q++) {                          for (int q = 0; q < numQuad; q++) {
151  #pragma ivdep  #pragma ivdep
152                              for (int l = 0; l < numComps; l++) {                              for (int l = 0; l < numComps; l++) {
# Line 309  void Assemble_gradient(Dudley_NodeFile* Line 158  void Assemble_gradient(Dudley_NodeFile*
158                      }                      }
159                  }                  }
160              } else if (numDim == 2) {              } else if (numDim == 2) {
161                  const dim_t numShapes = 3;                  const int numShapes = 3;
162  #pragma omp for  #pragma omp for
163                  for (index_t e = 0; e < elements->numElements; e++) {                  for (index_t e = 0; e < elements->numElements; e++) {
164                      double* grad_data_e = grad_data->getSampleDataRW(e);                      double* grad_data_e = grad_data.getSampleDataRW(e);
165                      memset(grad_data_e, 0, localGradSize);                      memset(grad_data_e, 0, localGradSize);
166                      for (int s = 0; s < numShapes; s++) {                      for (int s = 0; s < numShapes; s++) {
167                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
168                          const double* data_array = data->getSampleDataRO(nodes->reducedDegreesOfFreedomMapping->target[n]);                          const double* data_array = data.getSampleDataRO(target[n]);
169                          for (int q = 0; q < numQuad; q++) {                          for (int q = 0; q < numQuad; q++) {
170  #pragma ivdep  #pragma ivdep
171                              for (int l = 0; l < numComps; l++) {                              for (int l = 0; l < numComps; l++) {
# Line 331  void Assemble_gradient(Dudley_NodeFile* Line 180  void Assemble_gradient(Dudley_NodeFile*
180                      }                      }
181                  }                  }
182              } else if (numDim == 3) {              } else if (numDim == 3) {
183                  const dim_t numShapes = 4;                  const int numShapes = 4;
184  #pragma omp for  #pragma omp for
185                  for (index_t e = 0; e < elements->numElements; e++) {                  for (index_t e = 0; e < elements->numElements; e++) {
186                      double* grad_data_e = grad_data->getSampleDataRW(e);                      double* grad_data_e = grad_data.getSampleDataRW(e);
187                      memset(grad_data_e, 0, localGradSize);                      memset(grad_data_e, 0, localGradSize);
188                      for (int s = 0; s < numShapes; s++) {                      for (int s = 0; s < numShapes; s++) {
189                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];                          const index_t n = elements->Nodes[INDEX2(s, e, NN)];
190                          const double* data_array = data->getSampleDataRO(nodes->reducedDegreesOfFreedomMapping->target[n]);                          const double* data_array = data.getSampleDataRO(target[n]);
191                          for (int q = 0; q < numQuad; q++) {                          for (int q = 0; q < numQuad; q++) {
192  #pragma ivdep  #pragma ivdep
193                              for (int l = 0; l < numComps; l++) {                              for (int l = 0; l < numComps; l++) {

Legend:
 Removed from v.6078 changed lines Added in v.6079

 ViewVC Help Powered by ViewVC 1.1.26