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

Diff of /branches/trilinos_from_5897/dudley/src/Assemble_AverageElementData.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: averages data  
   
 *****************************************************************************/  
   
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_AverageElementData(Dudley_ElementFile* elements, escript::Data* out, const escript::Data* in)  void Assemble_AverageElementData(const ElementFile* elements,
24                                     escript::Data& out, const escript::Data& in)
25  {  {
26      if (!elements)      if (!elements)
27          return;          return;
28    
29      double wq;      double wq;
30      int numQuad_in, numQuad_out;      int numQuad_in, numQuad_out;
31      if (Assemble_reducedIntegrationOrder(in)) {      if (hasReducedIntegrationOrder(in)) {
32          numQuad_in = QuadNums[elements->numDim][0];          numQuad_in = QuadNums[elements->numDim][0];
33          wq = QuadWeight[elements->numDim][0];          wq = QuadWeight[elements->numDim][0];
34      } else {      } else {
35          numQuad_in = QuadNums[elements->numDim][1];          numQuad_in = QuadNums[elements->numDim][1];
36          wq = QuadWeight[elements->numDim][1];          wq = QuadWeight[elements->numDim][1];
37      }      }
38      if (Assemble_reducedIntegrationOrder(out)) {      if (hasReducedIntegrationOrder(out)) {
39          numQuad_out = QuadNums[elements->numDim][0];          numQuad_out = QuadNums[elements->numDim][0];
40      } else {      } else {
41          numQuad_out = QuadNums[elements->numDim][1];          numQuad_out = QuadNums[elements->numDim][1];
42      }      }
43    
44      const int numComps = out->getDataPointSize();      // check out and in
45      const dim_t numElements = elements->numElements;      const dim_t numElements = elements->numElements;
46        const int numComps = out.getDataPointSize();
47    
48      // check out and in      if (numComps != in.getDataPointSize()) {
     if (numComps != in->getDataPointSize()) {  
49          throw DudleyException("Assemble_AverageElementData: number of components of input and output Data do not match.");          throw DudleyException("Assemble_AverageElementData: number of components of input and output Data do not match.");
50      } else if (!in->numSamplesEqual(numQuad_in, numElements)) {      } else if (!in.numSamplesEqual(numQuad_in, numElements)) {
51          throw DudleyException("Assemble_AverageElementData: illegal number of samples of input Data object");          throw DudleyException("Assemble_AverageElementData: illegal number of samples of input Data object");
52      } else if (!out->numSamplesEqual(numQuad_out, numElements)) {      } else if (!out.numSamplesEqual(numQuad_out, numElements)) {
53          throw DudleyException("Assemble_AverageElementData: illegal number of samples of output Data object");          throw DudleyException("Assemble_AverageElementData: illegal number of samples of output Data object");
54      } else if (!out->actsExpanded()) {      } else if (!out.actsExpanded()) {
55          throw DudleyException("Assemble_AverageElementData: expanded Data object is expected for output data.");          throw DudleyException("Assemble_AverageElementData: expanded Data object is expected for output data.");
56      } else {      } else {
57          out->requireWrite();          out.requireWrite();
58          if (in->actsExpanded()) {          if (in.actsExpanded()) {
59              const double vol = wq * numQuad_in;              const double vol = wq * numQuad_in;
60              const double volinv = 1. / vol;              const double volinv = 1. / vol;
61  #pragma omp parallel for  #pragma omp parallel for
62              for (index_t n = 0; n < numElements; n++) {              for (index_t n = 0; n < numElements; n++) {
63                  const double* in_array = in->getSampleDataRO(n);                  const double* in_array = in.getSampleDataRO(n);
64                  double* out_array = out->getSampleDataRW(n);                  double* out_array = out.getSampleDataRW(n);
65                  for (int i = 0; i < numComps; ++i) {                  for (int i = 0; i < numComps; ++i) {
66                      double rtmp = 0.;                      double rtmp = 0.;
67                      for (int q = 0; q < numQuad_in; ++q)                      for (int q = 0; q < numQuad_in; ++q)
# Line 80  void Assemble_AverageElementData(Dudley_ Line 75  void Assemble_AverageElementData(Dudley_
75              const size_t numComps_size = numComps * sizeof(double);              const size_t numComps_size = numComps * sizeof(double);
76  #pragma omp parallel for  #pragma omp parallel for
77              for (index_t n = 0; n < numElements; n++) {              for (index_t n = 0; n < numElements; n++) {
78                  const double* in_array = in->getSampleDataRO(n);                  const double* in_array = in.getSampleDataRO(n);
79                  double* out_array = out->getSampleDataRW(n);                  double* out_array = out.getSampleDataRW(n);
80                  for (int q = 0; q < numQuad_out; q++)                  for (int q = 0; q < numQuad_out; q++)
81                      memcpy(out_array + q * numComps, in_array, numComps_size);                      memcpy(out_array + q * numComps, in_array, numComps_size);
82              }              }

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

  ViewVC Help
Powered by ViewVC 1.1.26