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

Diff of /trunk/speckley/src/Brick.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 1037  const dim_t* Brick::borrowSampleReferenc Line 1037  const dim_t* Brick::borrowSampleReferenc
1037          case Nodes:          case Nodes:
1038              return &m_nodeId[0];              return &m_nodeId[0];
1039          case Elements:          case Elements:
1040            case ReducedElements:
1041              return &m_elementId[0];              return &m_elementId[0];
1042          case Points:          case Points:
1043              return &m_diracPointNodeIDs[0];              return &m_diracPointNodeIDs[0];
# Line 1369  void Brick::interpolateElementsOnNodes(e Line 1370  void Brick::interpolateElementsOnNodes(e
1370      const dim_t max_x = m_NN[0];      const dim_t max_x = m_NN[0];
1371      const dim_t max_y = m_NN[1];      const dim_t max_y = m_NN[1];
1372      const dim_t max_z = m_NN[2];      const dim_t max_z = m_NN[2];
1373        const int inFS = in.getFunctionSpace().getTypeCode();
1374      out.requireWrite();      out.requireWrite();
1375      //init to zero so we can do some sums without undefined, may not be required      //init to zero so we can do some sums without undefined, may not be required
1376      memset(out.getSampleDataRW(0), 0, sizeof(double)*quads*quads*numComp);      memset(out.getSampleDataRW(0), 0, sizeof(double)*quads*quads*numComp);
1377      // the summation portion      // the summation portion
1378      for (dim_t colouring = 0; colouring < 2; colouring++) {      if (inFS == ReducedElements) {
1379  #pragma omp parallel for          for (dim_t colouring = 0; colouring < 2; colouring++) {
1380          for (dim_t ez = colouring; ez < NE2; ez += 2) {      #pragma omp parallel for
1381              for (dim_t ey = 0; ey < NE1; ey++) {              for (dim_t ez = colouring; ez < NE2; ez += 2) {
1382                  for (dim_t ex = 0; ex < NE0; ex++) {                  for (dim_t ey = 0; ey < NE1; ey++) {
1383                      dim_t start = m_order * (INDEX3(ex, ey, ez, max_x, max_y));                      for (dim_t ex = 0; ex < NE0; ex++) {
1384                      const double *e_in = in.getSampleDataRO(INDEX3(ex,ey,ez,NE0,NE1));                          dim_t start = m_order * (INDEX3(ex, ey, ez, max_x, max_y));
1385                      for (int qz = 0; qz < quads; qz++) {                          const double *e_in = in.getSampleDataRO(INDEX3(ex,ey,ez,NE0,NE1));
1386                          for (int qy = 0; qy < quads; qy++) {                          for (int qz = 0; qz < quads; qz++) {
1387                              for (int qx = 0; qx < quads; qx++) {                              for (int qy = 0; qy < quads; qy++) {
1388                                  double *n_out = out.getSampleDataRW(start + INDEX3(qx, qy, qz, max_x, max_y));                                  for (int qx = 0; qx < quads; qx++) {
1389                                  for (dim_t comp = 0; comp < numComp; comp++) {                                      double *n_out = out.getSampleDataRW(start + INDEX3(qx, qy, qz, max_x, max_y));
1390                                      n_out[comp] += e_in[INDEX4(comp, qx, qy, qz, numComp, quads, quads)];                                      for (dim_t comp = 0; comp < numComp; comp++) {
1391                                            n_out[comp] += e_in[comp];
1392                                        }
1393                                    }
1394                                }
1395                            }
1396                        }
1397                    }
1398                }
1399            }        
1400        } else {
1401            for (dim_t colouring = 0; colouring < 2; colouring++) {
1402        #pragma omp parallel for
1403                for (dim_t ez = colouring; ez < NE2; ez += 2) {
1404                    for (dim_t ey = 0; ey < NE1; ey++) {
1405                        for (dim_t ex = 0; ex < NE0; ex++) {
1406                            dim_t start = m_order * (INDEX3(ex, ey, ez, max_x, max_y));
1407                            const double *e_in = in.getSampleDataRO(INDEX3(ex,ey,ez,NE0,NE1));
1408                            for (int qz = 0; qz < quads; qz++) {
1409                                for (int qy = 0; qy < quads; qy++) {
1410                                    for (int qx = 0; qx < quads; qx++) {
1411                                        double *n_out = out.getSampleDataRW(start + INDEX3(qx, qy, qz, max_x, max_y));
1412                                        for (dim_t comp = 0; comp < numComp; comp++) {
1413                                            n_out[comp] += e_in[INDEX4(comp, qx, qy, qz, numComp, quads, quads)];
1414                                        }
1415                                  }                                  }
1416                              }                              }
1417                          }                          }
# Line 1446  void Brick::interpolateElementsOnNodes(e Line 1472  void Brick::interpolateElementsOnNodes(e
1472      }      }
1473  }  }
1474    
1475    void Brick::reduceElements(escript::Data& out, const escript::Data& in) const
1476    {
1477        if (m_order == 2) {
1478            reduction_order2(in, out);
1479        } else if (m_order == 3) {
1480            reduction_order3(in, out);
1481        } else if (m_order == 4) {
1482            reduction_order4(in, out);
1483        } else if (m_order == 5) {
1484            reduction_order5(in, out);
1485        } else if (m_order == 6) {
1486            reduction_order6(in, out);
1487        } else if (m_order == 7) {
1488            reduction_order7(in, out);
1489        } else if (m_order == 8) {
1490            reduction_order8(in, out);
1491        } else if (m_order == 9) {
1492            reduction_order9(in, out);
1493        } else if (m_order == 10) {
1494            reduction_order10(in, out);
1495        }
1496    }
1497    
1498  //protected  //protected
1499  void Brick::interpolateNodesOnElements(escript::Data& out,  void Brick::interpolateNodesOnElements(escript::Data& out,
1500                                         const escript::Data& in) const                                         const escript::Data& in,
1501                                           bool reduced) const
1502  {  {
1503        if (reduced) { //going to ReducedElements
1504            escript::Data funcIn(in, escript::function(*this));
1505            reduceElements(out, funcIn);
1506            return;
1507        }
1508      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1509      const dim_t NE0 = m_NE[0];      const dim_t NE0 = m_NE[0];
1510      const dim_t NE1 = m_NE[1];      const dim_t NE1 = m_NE[1];

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

  ViewVC Help
Powered by ViewVC 1.1.26