/[escript]/trunk/speckley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/speckley/src/Rectangle.cpp

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

revision 5477 by jfenwick, Sat Feb 14 00:25:03 2015 UTC revision 5478 by sshaw, Wed Feb 18 01:57:49 2015 UTC
# Line 873  const dim_t* Rectangle::borrowSampleRefe Line 873  const dim_t* Rectangle::borrowSampleRefe
873          case Nodes:          case Nodes:
874              return &m_nodeId[0];              return &m_nodeId[0];
875          case Elements:          case Elements:
876            case ReducedElements:
877              return &m_elementId[0];              return &m_elementId[0];
878          case Points:          case Points:
879              return &m_diracPointNodeIDs[0];              return &m_diracPointNodeIDs[0];
# Line 1146  void Rectangle::addToMatrixAndRHS(escrip Line 1147  void Rectangle::addToMatrixAndRHS(escrip
1147      throw SpeckleyException("Rectangle::addToMatrixAndRHS, adding to matrix not supported");      throw SpeckleyException("Rectangle::addToMatrixAndRHS, adding to matrix not supported");
1148  }  }
1149    
1150    void Rectangle::reduceElements(escript::Data& out, const escript::Data& in) const
1151    {
1152        if (m_order == 2) {
1153            reduction_order2(in, out);
1154        } else if (m_order == 3) {
1155            reduction_order3(in, out);
1156        } else if (m_order == 4) {
1157            reduction_order4(in, out);
1158        } else if (m_order == 5) {
1159            reduction_order5(in, out);
1160        } else if (m_order == 6) {
1161            reduction_order6(in, out);
1162        } else if (m_order == 7) {
1163            reduction_order7(in, out);
1164        } else if (m_order == 8) {
1165            reduction_order8(in, out);
1166        } else if (m_order == 9) {
1167            reduction_order9(in, out);
1168        } else if (m_order == 10) {
1169            reduction_order10(in, out);
1170        }
1171    }
1172    
1173  //protected  //protected
1174  void Rectangle::interpolateNodesOnElements(escript::Data& out,  void Rectangle::interpolateNodesOnElements(escript::Data& out,
1175                                             const escript::Data& in) const                                             const escript::Data& in,
1176                                               bool reduced) const
1177  {  {
1178      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1179      const dim_t NE0 = m_NE[0];      const dim_t NE0 = m_NE[0];
# Line 1156  void Rectangle::interpolateNodesOnElemen Line 1181  void Rectangle::interpolateNodesOnElemen
1181      const int quads = m_order + 1;      const int quads = m_order + 1;
1182      const int max_x = m_NN[0];      const int max_x = m_NN[0];
1183      out.requireWrite();      out.requireWrite();
1184        if (reduced) { //going to ReducedElements
1185            escript::Data funcIn(in, escript::function(*this));
1186            reduceElements(out, funcIn);
1187            return;
1188        }
1189  #pragma omp parallel for  #pragma omp parallel for
1190      for (dim_t ey = 0; ey < NE1; ey++) {      for (dim_t ey = 0; ey < NE1; ey++) {
1191          for (dim_t ex = 0; ex < NE0; ex++) {          for (dim_t ex = 0; ex < NE0; ex++) {
# Line 1277  void Rectangle::interpolateElementsOnNod Line 1307  void Rectangle::interpolateElementsOnNod
1307      const dim_t max_x = (m_order*NE0) + 1;      const dim_t max_x = (m_order*NE0) + 1;
1308      const dim_t max_y = (m_order*NE1) + 1;      const dim_t max_y = (m_order*NE1) + 1;
1309      out.requireWrite();      out.requireWrite();
1310        const int inFS = in.getFunctionSpace().getTypeCode();
1311      // the summation portion      // the summation portion
1312      for (dim_t colouring = 0; colouring < 2; colouring++) {      if (inFS == ReducedElements) {
1313            for (dim_t colouring = 0; colouring < 2; colouring++) {
1314    #pragma omp parallel for
1315                for (dim_t ey = colouring; ey < NE1; ey += 2) {
1316                    for (dim_t ex = 0; ex < NE0; ex++) {
1317                        dim_t start = ex*m_order + ey*max_x*m_order;
1318                        const double *e_in = in.getSampleDataRO(ex + ey*NE0);
1319                        for (int qy = 0; qy < quads; qy++) {
1320                            for (int qx = 0; qx < quads; qx++) {
1321                                double *n_out = out.getSampleDataRW(start + max_x*qy + qx);
1322                                for (int comp = 0; comp < numComp; comp++) {
1323                                    n_out[comp] += e_in[comp];
1324                                }
1325                            }
1326                        }
1327                    }
1328                }
1329            }
1330        } else { //inFS == Elements
1331            for (dim_t colouring = 0; colouring < 2; colouring++) {
1332  #pragma omp parallel for  #pragma omp parallel for
1333          for (dim_t ey = colouring; ey < NE1; ey += 2) {              for (dim_t ey = colouring; ey < NE1; ey += 2) {
1334              for (dim_t ex = 0; ex < NE0; ex++) {                  for (dim_t ex = 0; ex < NE0; ex++) {
1335                  dim_t start = ex*m_order + ey*max_x*m_order;                      dim_t start = ex*m_order + ey*max_x*m_order;
1336                  const double *e_in = in.getSampleDataRO(ex + ey*NE0);                      const double *e_in = in.getSampleDataRO(ex + ey*NE0);
1337                  for (int qy = 0; qy < quads; qy++) {                      for (int qy = 0; qy < quads; qy++) {
1338                      for (int qx = 0; qx < quads; qx++) {                          for (int qx = 0; qx < quads; qx++) {
1339                          double *n_out = out.getSampleDataRW(start + max_x*qy + qx);                              double *n_out = out.getSampleDataRW(start + max_x*qy + qx);
1340                          for (int comp = 0; comp < numComp; comp++) {                              for (int comp = 0; comp < numComp; comp++) {
1341                              n_out[comp] += e_in[INDEX3(comp, qx, qy, numComp, quads)];                                  n_out[comp] += e_in[INDEX3(comp, qx, qy, numComp, quads)];
1342                                }
1343                          }                          }
1344                      }                      }
1345                  }                  }

Legend:
Removed from v.5477  
changed lines
  Added in v.5478

  ViewVC Help
Powered by ViewVC 1.1.26