/[escript]/branches/diaplayground/ripley/src/Brick.cpp
ViewVC logotype

Diff of /branches/diaplayground/ripley/src/Brick.cpp

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

revision 3748 by caltinay, Thu Dec 15 07:36:19 2011 UTC revision 3753 by caltinay, Tue Jan 3 09:01:49 2012 UTC
# Line 47  Brick::Brick(int n0, int n1, int n2, dou Line 47  Brick::Brick(int n0, int n1, int n2, dou
47      if (m_NX*m_NY*m_NZ != m_mpiInfo->size)      if (m_NX*m_NY*m_NZ != m_mpiInfo->size)
48          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
49    
50      if (n0%m_NX > 0 || n1%m_NY > 0 || n2%m_NZ > 0)      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0 || (n2+1)%m_NZ > 0)
51          throw RipleyException("Number of elements must be separable into number of ranks in each dimension");          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");
52    
53      // local number of elements      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2) || (m_NZ > 1 && (n2+1)/m_NZ<2))
54      m_NE0 = n0/m_NX;          throw RipleyException("Too few elements for the number of ranks");
55      m_NE1 = n1/m_NY;  
56      m_NE2 = n2/m_NZ;      // local number of elements (including overlap)
57      // local number of nodes (not necessarily owned)      m_NE0 = (m_NX>1 ? (n0+1)/m_NX : n0);
58        if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)
59            m_NE0++;
60        m_NE1 = (m_NY>1 ? (n1+1)/m_NY : n1);
61        if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX>0 && m_mpiInfo->rank%(m_NX*m_NY)/m_NX<m_NY-1)
62            m_NE1++;
63        m_NE2 = (m_NZ>1 ? (n2+1)/m_NZ : n2);
64        if (m_mpiInfo->rank/(m_NX*m_NY)>0 && m_mpiInfo->rank/(m_NX*m_NY)<m_NZ-1)
65            m_NE2++;
66    
67        // local number of nodes
68      m_N0 = m_NE0+1;      m_N0 = m_NE0+1;
69      m_N1 = m_NE1+1;      m_N1 = m_NE1+1;
70      m_N2 = m_NE2+1;      m_N2 = m_NE2+1;
71    
72      // bottom-left-front node is at (offset0,offset1,offset2) in global mesh      // bottom-left-front node is at (offset0,offset1,offset2) in global mesh
73      m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);      m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);
74      m_offset1 = m_NE1*(m_mpiInfo->rank%(m_NX*m_NY)/m_NX);      if (m_offset0 > 0)
75      m_offset2 = m_NE2*(m_mpiInfo->rank/(m_NX*m_NY));          m_offset0--;
76        m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank%(m_NX*m_NY)/m_NX);
77        if (m_offset1 > 0)
78            m_offset1--;
79        m_offset2 = (n2+1)/m_NZ*(m_mpiInfo->rank/(m_NX*m_NY));
80        if (m_offset2 > 0)
81            m_offset2--;
82    
83      populateSampleIds();      populateSampleIds();
84  }  }
85    
# Line 239  const int* Brick::borrowSampleReferenceI Line 257  const int* Brick::borrowSampleReferenceI
257          case Nodes:          case Nodes:
258          case ReducedNodes: //FIXME: reduced          case ReducedNodes: //FIXME: reduced
259              return &m_nodeId[0];              return &m_nodeId[0];
260            case DegreesOfFreedom: //FIXME
261            case ReducedDegreesOfFreedom: //FIXME
262                return &m_dofId[0];
263          case Elements:          case Elements:
264          case ReducedElements:          case ReducedElements:
265              return &m_elementId[0];              return &m_elementId[0];
# Line 258  bool Brick::ownSample(int fsCode, index_ Line 279  bool Brick::ownSample(int fsCode, index_
279  {  {
280  #ifdef ESYS_MPI  #ifdef ESYS_MPI
281      if (fsCode == Nodes) {      if (fsCode == Nodes) {
282          const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];          const index_t left = (m_offset0==0 ? 0 : 1);
283          const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;          const index_t bottom = (m_offset1==0 ? 0 : 1);
284          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);          const index_t front = (m_offset2==0 ? 0 : 1);
285      } else          const index_t right = (m_mpiInfo->rank%m_NX==m_NX-1 ? m_N0 : m_N0-1);
286          throw RipleyException("ownSample() only implemented for Nodes");          const index_t top = (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1 ? m_N1 : m_N1-1);
287            const index_t back = (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1 ? m_N2 : m_N2-1);
288            const index_t x=id%m_N0;
289            const index_t y=id%(m_N0*m_N1)/m_N0;
290            const index_t z=id/(m_N0*m_N1);
291            return (x>=left && x<right && y>=bottom && y<top && z>=front && z<back);
292    
293        } else {
294            stringstream msg;
295            msg << "ownSample() not implemented for "
296                << functionSpaceTypeAsString(fsCode);
297            throw RipleyException(msg.str());
298        }
299  #else  #else
300      return true;      return true;
301  #endif  #endif
# Line 1226  pair<double,double> Brick::getFirstCoord Line 1259  pair<double,double> Brick::getFirstCoord
1259  }  }
1260    
1261  //protected  //protected
1262    dim_t Brick::getNumDOF() const
1263    {
1264        return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY*(m_gNE2+1)/m_NZ;
1265    }
1266    
1267    //protected
1268  dim_t Brick::getNumFaceElements() const  dim_t Brick::getNumFaceElements() const
1269  {  {
1270        const IndexVector faces = getNumFacesPerBoundary();
1271      dim_t n=0;      dim_t n=0;
1272      //left      for (size_t i=0; i<faces.size(); i++)
1273      if (m_offset0==0)          n+=faces[i];
         n+=m_NE1*m_NE2;  
     //right  
     if (m_mpiInfo->rank%m_NX==m_NX-1)  
         n+=m_NE1*m_NE2;  
     //bottom  
     if (m_offset1==0)  
         n+=m_NE0*m_NE2;  
     //top  
     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)  
         n+=m_NE0*m_NE2;  
     //front  
     if (m_offset2==0)  
         n+=m_NE0*m_NE1;  
     //back  
     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)  
         n+=m_NE0*m_NE1;  
   
1274      return n;      return n;
1275  }  }
1276    
# Line 1282  void Brick::assembleCoordinates(escript: Line 1305  void Brick::assembleCoordinates(escript:
1305  void Brick::populateSampleIds()  void Brick::populateSampleIds()
1306  {  {
1307      // identifiers are ordered from left to right, bottom to top, front to back      // identifiers are ordered from left to right, bottom to top, front to back
1308      // on each rank, except for the shared nodes which are owned by the rank      // globally
     // below / to the left / to the front of the current rank  
1309    
1310      // build node distribution vector first.      // build node distribution vector first.
     // m_nodeDistribution[i] is the first node id on rank i, that is  
1311      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
1312      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1313      m_nodeDistribution[1]=getNumNodes();      const dim_t numDOF=getNumDOF();
1314      for (dim_t k=1; k<m_mpiInfo->size-1; k++) {      for (dim_t k=1; k<m_mpiInfo->size; k++) {
1315          const index_t x = k%m_NX;          m_nodeDistribution[k]=k*numDOF;
         const index_t y = k%(m_NX*m_NY)/m_NX;  
         const index_t z = k/(m_NX*m_NY);  
         index_t numNodes=getNumNodes();  
         if (x>0)  
             numNodes-=m_N1*m_N2;  
         if (y>0)  
             numNodes-=m_N0*m_N2;  
         if (z>0)  
             numNodes-=m_N0*m_N1;  
         // if an edge was subtracted twice add it back  
         if (x>0 && y>0)  
             numNodes+=m_N2;  
         if (x>0 && z>0)  
             numNodes+=m_N1;  
         if (y>0 && z>0)  
             numNodes+=m_N0;  
         // the corner node was removed 3x and added back 3x, so subtract it  
         if (x>0 && y>0 && z>0)  
             numNodes--;  
         m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;  
1316      }      }
1317      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
   
1318      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
1319        m_dofId.resize(numDOF);
1320        m_elementId.resize(getNumElements());
1321        m_faceId.resize(getNumFaceElements());
1322    
1323      // the bottom, left and front planes are not owned by this rank so the  #pragma omp parallel
1324      // identifiers need to be computed accordingly      {
1325      const index_t left = (m_offset0==0 ? 0 : 1);  #pragma omp for nowait
1326      const index_t bottom = (m_offset1==0 ? 0 : 1);          // nodes
1327      const index_t front = (m_offset2==0 ? 0 : 1);          for (dim_t i2=0; i2<m_N2; i2++) {
1328                for (dim_t i1=0; i1<m_N1; i1++) {
1329      // case 1: all nodes on left plane are owned by rank on the left                  for (dim_t i0=0; i0<m_N0; i0++) {
1330      if (left>0) {                      m_nodeId[i0+i1*m_N0+i2*m_N0*m_N1] =
1331          const int neighbour=m_mpiInfo->rank-1;                          (m_offset2+i2)*(m_gNE0+1)*(m_gNE1+1)
1332          const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);                          +(m_offset1+i1)*(m_gNE0+1)
1333          const index_t leftN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);                          +m_offset0+i0;
1334  #pragma omp parallel for                  }
         for (dim_t i2=front; i2<m_N2; i2++) {  
             for (dim_t i1=bottom; i1<m_N1; i1++) {  
                 m_nodeId[i1*m_N0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]  
                     + (i1-bottom+1)*leftN0  
                     + (i2-front)*leftN0*leftN1 - 1;  
             }  
         }  
     }  
     // case 2: all nodes on bottom plane are owned by rank below  
     if (bottom>0) {  
         const int neighbour=m_mpiInfo->rank-m_NX;  
         const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t bottomN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);  
 #pragma omp parallel for  
         for (dim_t i2=front; i2<m_N2; i2++) {  
             for (dim_t i0=left; i0<m_N0; i0++) {  
                 m_nodeId[i0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]  
                     + bottomN0*(bottomN1-1)  
                     + (i2-front)*bottomN0*bottomN1 + i0-left;  
             }  
         }  
     }  
     // case 3: all nodes on front plane are owned by rank in front  
     if (front>0) {  
         const int neighbour=m_mpiInfo->rank-m_NX*m_NY;  
         const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);  
         const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);  
 #pragma omp parallel for  
         for (dim_t i1=bottom; i1<m_N1; i1++) {  
             for (dim_t i0=left; i0<m_N0; i0++) {  
                 m_nodeId[i0+i1*m_N0]=m_nodeDistribution[neighbour]  
                     + N0*N1*(N2-1)+(i1-bottom)*N0 + i0-left;  
1335              }              }
1336          }          }
     }  
     // case 4: nodes on front bottom edge are owned by the corresponding rank  
     if (front>0 && bottom>0) {  
         const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1);  
         const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);  
         const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);  
 #pragma omp parallel for  
         for (dim_t i0=left; i0<m_N0; i0++) {  
             m_nodeId[i0]=m_nodeDistribution[neighbour]  
                 + N0*N1*(N2-1)+(N1-1)*N0 + i0-left;  
         }  
     }  
     // case 5: nodes on left bottom edge are owned by the corresponding rank  
     if (left>0 && bottom>0) {  
         const int neighbour=m_mpiInfo->rank-m_NX-1;  
         const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);  
 #pragma omp parallel for  
         for (dim_t i2=front; i2<m_N2; i2++) {  
             m_nodeId[i2*m_N0*m_N1]=m_nodeDistribution[neighbour]  
                 + (1+i2-front)*N0*N1-1;  
         }  
     }  
     // case 6: nodes on left front edge are owned by the corresponding rank  
     if (left>0 && front>0) {  
         const int neighbour=m_mpiInfo->rank-m_NX*m_NY-1;  
         const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);  
         const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);  
 #pragma omp parallel for  
         for (dim_t i1=bottom; i1<m_N1; i1++) {  
             m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]  
                 + N0*N1*(N2-1)+N0-1+(i1-bottom)*N0;  
         }  
     }  
     // case 7: bottom-left-front corner node owned by corresponding rank  
     if (left>0 && bottom>0 && front>0) {  
         const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1)-1;  
         const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);  
         const index_t N2=(neighbour/(m_NX*m_NY) == 0 ? m_N2 : m_N2-1);  
         m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1*N2-1;  
     }  
1337    
1338      // the rest of the id's are contiguous          // degrees of freedom
1339      const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];  #pragma omp for nowait
1340  #pragma omp parallel for          for (dim_t k=0; k<numDOF; k++)
1341      for (dim_t i2=front; i2<m_N2; i2++) {              m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;
1342          for (dim_t i1=bottom; i1<m_N1; i1++) {  
1343              for (dim_t i0=left; i0<m_N0; i0++) {          // elements
1344                  m_nodeId[i0+i1*m_N0+i2*m_N0*m_N1] = firstId+i0-left  #pragma omp for nowait
1345                      +(i1-bottom)*(m_N0-left)          for (dim_t k=0; k<getNumElements(); k++)
1346                      +(i2-front)*(m_N0-left)*(m_N1-bottom);              m_elementId[k]=k;
1347              }  
1348          }          // face elements
1349      }  #pragma omp for
1350            for (dim_t k=0; k<getNumFaceElements(); k++)
1351                m_faceId[k]=k;
1352        } // end parallel section
1353    
1354      m_nodeTags.assign(getNumNodes(), 0);      m_nodeTags.assign(getNumNodes(), 0);
1355      updateTagsInUse(Nodes);      updateTagsInUse(Nodes);
1356    
     // elements  
     m_elementId.resize(getNumElements());  
 #pragma omp parallel for  
     for (dim_t k=0; k<getNumElements(); k++) {  
         m_elementId[k]=k;  
     }  
1357      m_elementTags.assign(getNumElements(), 0);      m_elementTags.assign(getNumElements(), 0);
1358      updateTagsInUse(Elements);      updateTagsInUse(Elements);
1359    
     // face elements  
     m_faceId.resize(getNumFaceElements());  
 #pragma omp parallel for  
     for (dim_t k=0; k<getNumFaceElements(); k++) {  
         m_faceId[k]=k;  
     }  
   
1360      // generate face offset vector and set face tags      // generate face offset vector and set face tags
1361      const IndexVector facesPerEdge = getNumFacesPerBoundary();      const IndexVector facesPerEdge = getNumFacesPerBoundary();
1362      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20, FRONT=100, BACK=200;      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20, FRONT=100, BACK=200;
# Line 1749  void Brick::interpolateNodesOnFaces(escr Line 1666  void Brick::interpolateNodesOnFaces(escr
1666      }      }
1667  }  }
1668    
1669    //protected
1670    void Brick::nodesToDOF(escript::Data& out, escript::Data& in) const
1671    {
1672        const dim_t numComp = in.getDataPointSize();
1673        out.requireWrite();
1674    
1675        const index_t left = (m_offset0==0 ? 0 : 1);
1676        const index_t bottom = (m_offset1==0 ? 0 : 1);
1677        const index_t front = (m_offset2==0 ? 0 : 1);
1678        const index_t right = (m_mpiInfo->rank%m_NX==m_NX-1 ? m_N0 : m_N0-1);
1679        const index_t top = (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1 ? m_N1 : m_N1-1);
1680        const index_t back = (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1 ? m_N2 : m_N2-1);
1681        index_t n=0;
1682        for (index_t i=front; i<back; i++) {
1683            for (index_t j=bottom; j<top; j++) {
1684                for (index_t k=left; k<right; k++, n++) {
1685                    const double* src=in.getSampleDataRO(k+j*m_N0+i*m_N0*m_N1);
1686                    copy(src, src+numComp, out.getSampleDataRW(n));
1687                }
1688            }
1689        }
1690    }
1691    
1692    
1693  } // end of namespace ripley  } // end of namespace ripley
1694    

Legend:
Removed from v.3748  
changed lines
  Added in v.3753

  ViewVC Help
Powered by ViewVC 1.1.26