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

revision 6008 by caltinay, Mon Feb 22 06:59:27 2016 UTC revision 6009 by caltinay, Wed Mar 2 04:13:26 2016 UTC
# Line 14  Line 14
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16
17  /************************************************************************************/  /****************************************************************************
18
19  /*        assemblage routines: interpolates nodal data in a data array onto elements (=integration points) */    Assemblage routines: interpolates nodal data in a data array onto elements
20      (=integration points)
21
22  /************************************************************************************/  *****************************************************************************/

#define ESNEEDPYTHON
#include "esysUtils/first.h"
23
24  #include "Assemble.h"  #include "Assemble.h"
25  #include "ShapeTable.h"  #include "ShapeTable.h"
26  #include "Util.h"  #include "Util.h"
27
28  /************************************************************************************/  namespace dudley {
29
30  void Dudley_Assemble_interpolate(Dudley_NodeFile* nodes,  void Assemble_interpolate(Dudley_NodeFile* nodes,
31                                   Dudley_ElementFile* elements,                            Dudley_ElementFile* elements,
32                                   const escript::Data* data,                            const escript::Data* data,
33                                   escript::Data* interpolated_data)                            escript::Data* interpolated_data)
34  {  {
35      const double *data_array;      if (!nodes || !elements)
double *local_data = NULL;
bool reduced_integration = FALSE;
dim_t q, i, NS_DOF, NN, numNodes = 0, e, numQuad = 0;
dim_t numComps = getDataPointSize(data);
index_t *map = NULL;
const double *shapeFns = 0;
int data_type = getFunctionSpaceType(data);
size_t numComps_size;
Dudley_resetError();
if (nodes == NULL || elements == NULL)
36          return;          return;
reduced_integration = Dudley_Assemble_reducedIntegrationOrder(interpolated_data);
NN = elements->numNodes;
37
38      /* set some parameter */      const int data_type = data->getFunctionSpace().getTypeCode();
39        const bool reduced_integration = Assemble_reducedIntegrationOrder(interpolated_data);
40
41      if (data_type == DUDLEY_NODES)      dim_t numNodes = 0;
42      {      index_t *map = NULL;
43
44        if (data_type == DUDLEY_NODES) {
45          numNodes = Dudley_NodeFile_getNumNodes(nodes);          numNodes = Dudley_NodeFile_getNumNodes(nodes);
46          map = Dudley_NodeFile_borrowTargetNodes(nodes);          map = Dudley_NodeFile_borrowTargetNodes(nodes);
47      }      } else if (data_type == DUDLEY_REDUCED_NODES) {
else if (data_type == DUDLEY_REDUCED_NODES)
{
48          numNodes = Dudley_NodeFile_getNumReducedNodes(nodes);          numNodes = Dudley_NodeFile_getNumReducedNodes(nodes);
49          map = Dudley_NodeFile_borrowTargetReducedNodes(nodes);          map = Dudley_NodeFile_borrowTargetReducedNodes(nodes);
50      }      } else if (data_type == DUDLEY_DEGREES_OF_FREEDOM) {
51      else if (data_type == DUDLEY_DEGREES_OF_FREEDOM)          if (elements->MPIInfo->size > 1) {
52      {              throw DudleyException("Assemble_interpolate: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");
if (elements->MPIInfo->size > 1)
{
Dudley_setError(TYPE_ERROR,
"Dudley_Assemble_interpolate: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");
return;
53          }          }
54          numNodes = Dudley_NodeFile_getNumDegreesOfFreedom(nodes);          numNodes = Dudley_NodeFile_getNumDegreesOfFreedom(nodes);
55          map = Dudley_NodeFile_borrowTargetDegreesOfFreedom(nodes);          map = Dudley_NodeFile_borrowTargetDegreesOfFreedom(nodes);
56      }      } else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {
57      else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM)          if (elements->MPIInfo->size > 1) {
58      {              throw DudleyException("Assemble_interpolate: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
if (elements->MPIInfo->size > 1)
{
Dudley_setError(TYPE_ERROR,
"Dudley_Assemble_interpolate: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
return;
59          }          }
60          numNodes = Dudley_NodeFile_getNumReducedDegreesOfFreedom(nodes);          numNodes = Dudley_NodeFile_getNumReducedDegreesOfFreedom(nodes);
61          map = Dudley_NodeFile_borrowTargetReducedDegreesOfFreedom(nodes);          map = Dudley_NodeFile_borrowTargetReducedDegreesOfFreedom(nodes);
62      }      } else {
63      else          throw DudleyException("Assemble_interpolate: Cannot interpolate data");
{
Dudley_setError(TYPE_ERROR, "Dudley_Assemble_interpolate: Cannot interpolate data");
return;
}

numQuad = reduced_integration ? 1 : (elements->numDim + 1);
NS_DOF = elements->numDim + 1;

/* check the dimensions of interpolated_data and data */

{
Dudley_setError(TYPE_ERROR, "Dudley_Assemble_interpolate: illegal number of samples of output Data object");
}
else if (!numSamplesEqual(data, 1, numNodes))
{
Dudley_setError(TYPE_ERROR, "Dudley_Assemble_interpolate: illegal number of samples of input Data object");
}
else if (numComps != getDataPointSize(interpolated_data))
{
Dudley_setError(TYPE_ERROR,
"Dudley_Assemble_interpolate: number of components of input and interpolated Data do not match.");
}
else if (!isExpanded(interpolated_data))
{
Dudley_setError(TYPE_ERROR, "Dudley_Assemble_interpolate: expanded Data object is expected for output data.");
64      }      }
65
66      if (Dudley_noError() && !getQuadShape(elements->numDim, reduced_integration, &shapeFns))      const dim_t numComps = data->getDataPointSize();
67      {      const int NN = elements->numNodes;
68          Dudley_setError(TYPE_ERROR, "Dudley_Assemble_interpolate: unable to locate shape function.");      const int numQuad = reduced_integration ? 1 : (elements->numDim + 1);
69      }      const int NS_DOF = elements->numDim + 1;
70        const double *shapeFns = NULL;
71      /* now we can start */
72        // check the dimensions of interpolated_data and data
73      if (Dudley_noError())      if (!interpolated_data->numSamplesEqual(numQuad, elements->numElements)) {
74      {          throw DudleyException("Assemble_interpolate: illegal number of samples of output Data object");
75          requireWrite(interpolated_data);      } else if (!data->numSamplesEqual(1, numNodes)) {
76  #pragma omp parallel private(local_data, numComps_size)          throw DudleyException("Assemble_interpolate: illegal number of samples of input Data object");
77          {      } else if (numComps != interpolated_data->getDataPointSize()) {
78              local_data = NULL;          throw DudleyException("Assemble_interpolate: number of components of input and interpolated Data do not match.");
79              /* allocation of work arrays */      } else if (!interpolated_data->actsExpanded()) {
80              local_data = new double[NS_DOF * numComps];          throw DudleyException("Assemble_interpolate: expanded Data object is expected for output data.");
81              if (!Dudley_checkPtr(local_data))      }
82              {
83                  numComps_size = (size_t) numComps *sizeof(double);      if (!getQuadShape(elements->numDim, reduced_integration, &shapeFns)) {
84                  /* open the element loop */          throw DudleyException("Assemble_interpolate: unable to locate shape function.");
85  #pragma omp for private(e,q,i,data_array) schedule(static)      }
86                  for (e = 0; e < elements->numElements; e++)
87                  {      interpolated_data->requireWrite();
88                      for (q = 0; q < NS_DOF; q++)  #pragma omp parallel
89                      {      {
90                          i = elements->Nodes[INDEX2(q, e, NN)];          std::vector<double> local_data(NS_DOF * numComps);
91                          data_array = getSampleDataRO(data, map[i]);          const size_t numComps_size = numComps *sizeof(double);
92                          memcpy(&(local_data[INDEX3(0, q, 0, numComps, NS_DOF)]), data_array, numComps_size);          /* open the element loop */
93                      }  #pragma omp for
94                      /*  calculate interpolated_data=local_data*S */          for (index_t e = 0; e < elements->numElements; e++) {
95                      Dudley_Util_SmallMatSetMult1(1, numComps, numQuad, getSampleDataRW(interpolated_data, e),              for (int q = 0; q < NS_DOF; q++) {
96                                                   NS_DOF, local_data, /*basis->S */ shapeFns);                  const index_t i = elements->Nodes[INDEX2(q, e, NN)];
97                  }               /* end of element loop */                  const double* data_array = data->getSampleDataRO(map[i]);
98                    memcpy(&local_data[INDEX3(0, q, 0, numComps, NS_DOF)],
99                           data_array, numComps_size);
100              }              }
101              delete[] local_data;              // calculate interpolated_data=local_data*S
102          }                       /* end of parallel region */              Dudley_Util_SmallMatSetMult1(1, numComps, numQuad,
103      }                              interpolated_data->getSampleDataRW(e), NS_DOF,
104                                &local_data[0], shapeFns);
105            } // end of element loop
106        } // end of parallel region
107  }  }
108
109    } // namespace dudley
110

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