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

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

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

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 routines: interpolates nodal data in a data array onto elements  
   (=integration points)  
   
 *****************************************************************************/  
   
17  #include "Assemble.h"  #include "Assemble.h"
18  #include "ShapeTable.h"  #include "ShapeTable.h"
19  #include "Util.h"  #include "Util.h"
20    
21  namespace dudley {  namespace dudley {
22    
23  void Assemble_interpolate(Dudley_NodeFile* nodes,  void Assemble_interpolate(const NodeFile* nodes, const ElementFile* elements,
24                            Dudley_ElementFile* elements,                            const escript::Data& data,
25                            const escript::Data* data,                            escript::Data& interpolated_data)
                           escript::Data* interpolated_data)  
26  {  {
27      if (!nodes || !elements)      if (!nodes || !elements)
28          return;          return;
29    
30      const int data_type = data->getFunctionSpace().getTypeCode();      const int data_type = data.getFunctionSpace().getTypeCode();
31      const bool reduced_integration = Assemble_reducedIntegrationOrder(interpolated_data);      const bool reduced_integration = hasReducedIntegrationOrder(interpolated_data);
32    
33      dim_t numNodes = 0;      dim_t numNodes = 0;
34      index_t *map = NULL;      const index_t* map = NULL;
35    
36      if (data_type == DUDLEY_NODES) {      if (data_type == DUDLEY_NODES) {
37          numNodes = Dudley_NodeFile_getNumNodes(nodes);          numNodes = nodes->getNumNodes();
38          map = Dudley_NodeFile_borrowTargetNodes(nodes);          map = nodes->borrowTargetNodes();
     } else if (data_type == DUDLEY_REDUCED_NODES) {  
         numNodes = Dudley_NodeFile_getNumReducedNodes(nodes);  
         map = Dudley_NodeFile_borrowTargetReducedNodes(nodes);  
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_interpolate: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");              throw DudleyException("Assemble_interpolate: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");
42          }          }
43          numNodes = Dudley_NodeFile_getNumDegreesOfFreedom(nodes);          numNodes = nodes->getNumDegreesOfFreedom();
44          map = Dudley_NodeFile_borrowTargetDegreesOfFreedom(nodes);          map = nodes->borrowTargetDegreesOfFreedom();
     } else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {  
         if (elements->MPIInfo->size > 1) {  
             throw DudleyException("Assemble_interpolate: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");  
         }  
         numNodes = Dudley_NodeFile_getNumReducedDegreesOfFreedom(nodes);  
         map = Dudley_NodeFile_borrowTargetReducedDegreesOfFreedom(nodes);  
45      } else {      } else {
46          throw DudleyException("Assemble_interpolate: Cannot interpolate data");          throw DudleyException("Assemble_interpolate: Cannot interpolate data");
47      }      }
48    
49      const dim_t numComps = data->getDataPointSize();      const dim_t numComps = data.getDataPointSize();
50      const int NN = elements->numNodes;      const int NN = elements->numNodes;
51      const int numQuad = reduced_integration ? 1 : (elements->numDim + 1);      const int numQuad = reduced_integration ? 1 : elements->numNodes;
52      const int NS_DOF = elements->numDim + 1;      const int NS_DOF = elements->numDim + 1;
53      const double *shapeFns = NULL;      const double *shapeFns = NULL;
54    
55      // check the dimensions of interpolated_data and data      // check the dimensions of interpolated_data and data
56      if (!interpolated_data->numSamplesEqual(numQuad, elements->numElements)) {      if (!interpolated_data.numSamplesEqual(numQuad, elements->numElements)) {
57          throw DudleyException("Assemble_interpolate: illegal number of samples of output Data object");          throw DudleyException("Assemble_interpolate: illegal number of samples of output Data object");
58      } else if (!data->numSamplesEqual(1, numNodes)) {      } else if (!data.numSamplesEqual(1, numNodes)) {
59          throw DudleyException("Assemble_interpolate: illegal number of samples of input Data object");          throw DudleyException("Assemble_interpolate: illegal number of samples of input Data object");
60      } else if (numComps != interpolated_data->getDataPointSize()) {      } else if (numComps != interpolated_data.getDataPointSize()) {
61          throw DudleyException("Assemble_interpolate: number of components of input and interpolated Data do not match.");          throw DudleyException("Assemble_interpolate: number of components of input and interpolated Data do not match.");
62      } else if (!interpolated_data->actsExpanded()) {      } else if (!interpolated_data.actsExpanded()) {
63          throw DudleyException("Assemble_interpolate: expanded Data object is expected for output data.");          throw DudleyException("Assemble_interpolate: expanded Data object is expected for output data.");
64      }      }
65    
# Line 84  void Assemble_interpolate(Dudley_NodeFil Line 67  void Assemble_interpolate(Dudley_NodeFil
67          throw DudleyException("Assemble_interpolate: unable to locate shape function.");          throw DudleyException("Assemble_interpolate: unable to locate shape function.");
68      }      }
69    
70      interpolated_data->requireWrite();      interpolated_data.requireWrite();
71  #pragma omp parallel  #pragma omp parallel
72      {      {
73          std::vector<double> local_data(NS_DOF * numComps);          std::vector<double> local_data(NS_DOF * numComps);
74          const size_t numComps_size = numComps *sizeof(double);          const size_t numComps_size = numComps *sizeof(double);
75          /* open the element loop */          // open the element loop
76  #pragma omp for  #pragma omp for
77          for (index_t e = 0; e < elements->numElements; e++) {          for (index_t e = 0; e < elements->numElements; e++) {
78              for (int q = 0; q < NS_DOF; q++) {              for (int q = 0; q < NS_DOF; q++) {
79                  const index_t i = elements->Nodes[INDEX2(q, e, NN)];                  const index_t i = elements->Nodes[INDEX2(q, e, NN)];
80                  const double* data_array = data->getSampleDataRO(map[i]);                  const double* data_array = data.getSampleDataRO(map[i]);
81                  memcpy(&local_data[INDEX3(0, q, 0, numComps, NS_DOF)],                  memcpy(&local_data[INDEX3(0, q, 0, numComps, NS_DOF)],
82                         data_array, numComps_size);                         data_array, numComps_size);
83              }              }
84              // calculate interpolated_data=local_data*S              // calculate interpolated_data=local_data*S
85              Dudley_Util_SmallMatSetMult1(1, numComps, numQuad,              util::smallMatSetMult1(1, numComps, numQuad,
86                              interpolated_data->getSampleDataRW(e), NS_DOF,                              interpolated_data.getSampleDataRW(e), NS_DOF,
87                              &local_data[0], shapeFns);                              &local_data[0], shapeFns);
88          } // end of element loop          } // end of element loop
89      } // end of parallel region      } // end of parallel region

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

  ViewVC Help
Powered by ViewVC 1.1.26