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

Diff of /trunk/ripley/src/Brick.cpp

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

revision 3755 by caltinay, Thu Jan 5 06:51:31 2012 UTC revision 3756 by caltinay, Fri Jan 6 02:35:19 2012 UTC
# Line 81  Brick::Brick(int n0, int n1, int n2, dou Line 81  Brick::Brick(int n0, int n1, int n2, dou
81          m_offset2--;          m_offset2--;
82    
83      populateSampleIds();      populateSampleIds();
84        createPattern();
85  }  }
86    
87    
# Line 279  bool Brick::ownSample(int fsCode, index_ Line 280  bool Brick::ownSample(int fsCode, index_
280  {  {
281  #ifdef ESYS_MPI  #ifdef ESYS_MPI
282      if (fsCode == Nodes) {      if (fsCode == Nodes) {
283          const index_t left = (m_offset0==0 ? 0 : 1);          return (m_dofMap[id] < getNumDOF());
         const index_t bottom = (m_offset1==0 ? 0 : 1);  
         const index_t front = (m_offset2==0 ? 0 : 1);  
         const index_t right = (m_mpiInfo->rank%m_NX==m_NX-1 ? m_N0 : m_N0-1);  
         const index_t top = (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1 ? m_N1 : m_N1-1);  
         const index_t back = (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1 ? m_N2 : m_N2-1);  
         const index_t x=id%m_N0;  
         const index_t y=id%(m_N0*m_N1)/m_N0;  
         const index_t z=id/(m_N0*m_N1);  
         return (x>=left && x<right && y>=bottom && y<top && z>=front && z<back);  
   
284      } else {      } else {
285          stringstream msg;          stringstream msg;
286          msg << "ownSample() not implemented for "          msg << "ownSample() not implemented for "
# Line 1182  Paso_SystemMatrixPattern* Brick::getPatt Line 1173  Paso_SystemMatrixPattern* Brick::getPatt
1173      if (reducedRowOrder || reducedColOrder)      if (reducedRowOrder || reducedColOrder)
1174          throw RipleyException("getPattern() not implemented for reduced order");          throw RipleyException("getPattern() not implemented for reduced order");
1175    
1176      // connector      return m_pattern;
1177      RankVector neighbour;  }
1178      IndexVector offsetInShared(1,0);  
1179      IndexVector sendShared, recvShared;  void Brick::Print_Mesh_Info(const bool full) const
1180      const IndexVector faces=getNumFacesPerBoundary();  {
1181      const index_t nDOF0 = (m_gNE0+1)/m_NX;      RipleyDomain::Print_Mesh_Info(full);
1182      const index_t nDOF1 = (m_gNE1+1)/m_NY;      if (full) {
1183      const index_t nDOF2 = (m_gNE2+1)/m_NZ;          cout << "     Id  Coordinates" << endl;
1184      const int numDOF=nDOF0*nDOF1*nDOF2;          cout.precision(15);
1185            cout.setf(ios::scientific, ios::floatfield);
1186            pair<double,double> xdx = getFirstCoordAndSpacing(0);
1187            pair<double,double> ydy = getFirstCoordAndSpacing(1);
1188            pair<double,double> zdz = getFirstCoordAndSpacing(2);
1189            for (index_t i=0; i < getNumNodes(); i++) {
1190                cout << "  " << setw(5) << m_nodeId[i]
1191                    << "  " << xdx.first+(i%m_N0)*xdx.second
1192                    << "  " << ydy.first+(i%(m_N0*m_N1)/m_N0)*ydy.second
1193                    << "  " << zdz.first+(i/(m_N0*m_N1))*zdz.second << endl;
1194            }
1195        }
1196    }
1197    
1198    IndexVector Brick::getNumNodesPerDim() const
1199    {
1200        IndexVector ret;
1201        ret.push_back(m_N0);
1202        ret.push_back(m_N1);
1203        ret.push_back(m_N2);
1204        return ret;
1205    }
1206    
1207    IndexVector Brick::getNumElementsPerDim() const
1208    {
1209        IndexVector ret;
1210        ret.push_back(m_NE0);
1211        ret.push_back(m_NE1);
1212        ret.push_back(m_NE2);
1213        return ret;
1214    }
1215    
1216    IndexVector Brick::getNumFacesPerBoundary() const
1217    {
1218        IndexVector ret(6, 0);
1219        //left
1220        if (m_offset0==0)
1221            ret[0]=m_NE1*m_NE2;
1222        //right
1223        if (m_mpiInfo->rank%m_NX==m_NX-1)
1224            ret[1]=m_NE1*m_NE2;
1225        //bottom
1226        if (m_offset1==0)
1227            ret[2]=m_NE0*m_NE2;
1228        //top
1229        if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
1230            ret[3]=m_NE0*m_NE2;
1231        //front
1232        if (m_offset2==0)
1233            ret[4]=m_NE0*m_NE1;
1234        //back
1235        if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
1236            ret[5]=m_NE0*m_NE1;
1237        return ret;
1238    }
1239    
1240    pair<double,double> Brick::getFirstCoordAndSpacing(dim_t dim) const
1241    {
1242        if (dim==0)
1243            return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
1244        else if (dim==1)
1245            return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
1246        else if (dim==2)
1247            return pair<double,double>((m_l2*m_offset2)/m_gNE2, m_l2/m_gNE2);
1248    
1249        throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
1250    }
1251    
1252    //protected
1253    dim_t Brick::getNumDOF() const
1254    {
1255        return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY*(m_gNE2+1)/m_NZ;
1256    }
1257    
1258    //protected
1259    dim_t Brick::getNumFaceElements() const
1260    {
1261        const IndexVector faces = getNumFacesPerBoundary();
1262        dim_t n=0;
1263        for (size_t i=0; i<faces.size(); i++)
1264            n+=faces[i];
1265        return n;
1266    }
1267    
1268    //protected
1269    void Brick::assembleCoordinates(escript::Data& arg) const
1270    {
1271        escriptDataC x = arg.getDataC();
1272        int numDim = m_numDim;
1273        if (!isDataPointShapeEqual(&x, 1, &numDim))
1274            throw RipleyException("setToX: Invalid Data object shape");
1275        if (!numSamplesEqual(&x, 1, getNumNodes()))
1276            throw RipleyException("setToX: Illegal number of samples in Data object");
1277    
1278        pair<double,double> xdx = getFirstCoordAndSpacing(0);
1279        pair<double,double> ydy = getFirstCoordAndSpacing(1);
1280        pair<double,double> zdz = getFirstCoordAndSpacing(2);
1281        arg.requireWrite();
1282    #pragma omp parallel for
1283        for (dim_t i2 = 0; i2 < m_N2; i2++) {
1284            for (dim_t i1 = 0; i1 < m_N1; i1++) {
1285                for (dim_t i0 = 0; i0 < m_N0; i0++) {
1286                    double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);
1287                    point[0] = xdx.first+i0*xdx.second;
1288                    point[1] = ydy.first+i1*ydy.second;
1289                    point[2] = zdz.first+i2*zdz.second;
1290                }
1291            }
1292        }
1293    }
1294    
1295    //protected
1296    dim_t Brick::insertNeighbourNodes(IndexVector& index, index_t node) const
1297    {
1298        const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1299        const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1300        const dim_t nDOF2 = (m_gNE2+1)/m_NZ;
1301        const int x=node%nDOF0;
1302        const int y=node%(nDOF0*nDOF1)/nDOF0;
1303        const int z=node/(nDOF0*nDOF1);
1304        int num=0;
1305        // loop through potential neighbours and add to index if positions are
1306        // within bounds
1307        for (int i2=-1; i2<2; i2++) {
1308            for (int i1=-1; i1<2; i1++) {
1309                for (int i0=-1; i0<2; i0++) {
1310                    // skip node itself
1311                    if (i0==0 && i1==0 && i2==0)
1312                        continue;
1313                    // location of neighbour node
1314                    const int nx=x+i0;
1315                    const int ny=y+i1;
1316                    const int nz=z+i2;
1317                    if (nx>=0 && ny>=0 && nz>=0
1318                            && nx<nDOF0 && ny<nDOF1 && nz<nDOF2) {
1319                        index.push_back(nz*nDOF0*nDOF1+ny*nDOF0+nx);
1320                        num++;
1321                    }
1322                }
1323            }
1324        }
1325    
1326        return num;
1327    }
1328    
1329    //protected
1330    void Brick::nodesToDOF(escript::Data& out, escript::Data& in) const
1331    {
1332        const dim_t numComp = in.getDataPointSize();
1333        out.requireWrite();
1334    
1335        const index_t left = (m_offset0==0 ? 0 : 1);
1336        const index_t bottom = (m_offset1==0 ? 0 : 1);
1337        const index_t front = (m_offset2==0 ? 0 : 1);
1338        const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1339        const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1340        const dim_t nDOF2 = (m_gNE2+1)/m_NZ;
1341    #pragma omp parallel for
1342        for (index_t i=0; i<nDOF2; i++) {
1343            for (index_t j=0; j<nDOF1; j++) {
1344                for (index_t k=0; k<nDOF0; k++) {
1345                    const index_t n=k+left+(j+bottom)*m_N0+(i+front)*m_N0*m_N1;
1346                    const double* src=in.getSampleDataRO(n);
1347                    copy(src, src+numComp, out.getSampleDataRW(k+j*nDOF0+i*nDOF0*nDOF1));
1348                }
1349            }
1350        }
1351    }
1352    
1353    //protected
1354    void Brick::dofToNodes(escript::Data& out, escript::Data& in) const
1355    {
1356        const dim_t numComp = in.getDataPointSize();
1357        Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);
1358        in.requireWrite();
1359        Paso_Coupler_startCollect(coupler, in.getSampleDataRW(0));
1360    
1361        const dim_t numDOF = getNumDOF();
1362        out.requireWrite();
1363        const double* buffer = Paso_Coupler_finishCollect(coupler);
1364    
1365    #pragma omp parallel for
1366        for (index_t i=0; i<getNumNodes(); i++) {
1367            const double* src=(m_dofMap[i]<numDOF ?
1368                    in.getSampleDataRO(m_dofMap[i])
1369                    : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1370            copy(src, src+numComp, out.getSampleDataRW(i));
1371        }
1372    }
1373    
1374    //private
1375    void Brick::populateSampleIds()
1376    {
1377        // identifiers are ordered from left to right, bottom to top, front to back
1378        // globally
1379    
1380        // build node distribution vector first.
1381        // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
1382        m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1383        const dim_t numDOF=getNumDOF();
1384        for (dim_t k=1; k<m_mpiInfo->size; k++) {
1385            m_nodeDistribution[k]=k*numDOF;
1386        }
1387        m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
1388        m_nodeId.resize(getNumNodes());
1389        m_dofId.resize(numDOF);
1390        m_elementId.resize(getNumElements());
1391        m_faceId.resize(getNumFaceElements());
1392    
1393    #pragma omp parallel
1394        {
1395    #pragma omp for nowait
1396            // nodes
1397            for (dim_t i2=0; i2<m_N2; i2++) {
1398                for (dim_t i1=0; i1<m_N1; i1++) {
1399                    for (dim_t i0=0; i0<m_N0; i0++) {
1400                        m_nodeId[i0+i1*m_N0+i2*m_N0*m_N1] =
1401                            (m_offset2+i2)*(m_gNE0+1)*(m_gNE1+1)
1402                            +(m_offset1+i1)*(m_gNE0+1)
1403                            +m_offset0+i0;
1404                    }
1405                }
1406            }
1407    
1408            // degrees of freedom
1409    #pragma omp for nowait
1410            for (dim_t k=0; k<numDOF; k++)
1411                m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;
1412    
1413            // elements
1414    #pragma omp for nowait
1415            for (dim_t i2=0; i2<m_NE2; i2++) {
1416                for (dim_t i1=0; i1<m_NE1; i1++) {
1417                    for (dim_t i0=0; i0<m_NE0; i0++) {
1418                        m_elementId[i0+i1*m_NE0+i2*m_NE0*m_NE1] =
1419                            (m_offset2+i2)*m_gNE0*m_gNE1
1420                            +(m_offset1+i1)*m_gNE0
1421                            +m_offset0+i0;
1422                    }
1423                }
1424            }
1425    
1426            // face elements
1427    #pragma omp for
1428            for (dim_t k=0; k<getNumFaceElements(); k++)
1429                m_faceId[k]=k;
1430        } // end parallel section
1431    
1432        m_nodeTags.assign(getNumNodes(), 0);
1433        updateTagsInUse(Nodes);
1434    
1435        m_elementTags.assign(getNumElements(), 0);
1436        updateTagsInUse(Elements);
1437    
1438        // generate face offset vector and set face tags
1439        const IndexVector facesPerEdge = getNumFacesPerBoundary();
1440        const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20, FRONT=100, BACK=200;
1441        const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP, FRONT, BACK };
1442        m_faceOffset.assign(facesPerEdge.size(), -1);
1443        m_faceTags.clear();
1444        index_t offset=0;
1445        for (size_t i=0; i<facesPerEdge.size(); i++) {
1446            if (facesPerEdge[i]>0) {
1447                m_faceOffset[i]=offset;
1448                offset+=facesPerEdge[i];
1449                m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);
1450            }
1451        }
1452        setTagMap("left", LEFT);
1453        setTagMap("right", RIGHT);
1454        setTagMap("bottom", BOTTOM);
1455        setTagMap("top", TOP);
1456        setTagMap("front", FRONT);
1457        setTagMap("back", BACK);
1458        updateTagsInUse(FaceElements);
1459    }
1460    
1461    //private
1462    void Brick::createPattern()
1463    {
1464        const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1465        const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1466        const dim_t nDOF2 = (m_gNE2+1)/m_NZ;
1467      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
1468      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
1469      const index_t front = (m_offset2==0 ? 0 : 1);      const index_t front = (m_offset2==0 ? 0 : 1);
     vector<IndexVector> colIndices(numDOF); // for the couple blocks  
     int numShared=0;  
1470    
1471        // populate node->DOF mapping with own degrees of freedom.
1472        // The rest is assigned in the loop further down
1473      m_dofMap.assign(getNumNodes(), 0);      m_dofMap.assign(getNumNodes(), 0);
1474  #pragma omp parallel for  #pragma omp parallel for
1475      for (index_t i=front; i<m_N2; i++) {      for (index_t i=front; i<m_N2; i++) {
# Line 1210  Paso_SystemMatrixPattern* Brick::getPatt Line 1483  Paso_SystemMatrixPattern* Brick::getPatt
1483      // build list of shared components and neighbours by looping through      // build list of shared components and neighbours by looping through
1484      // all potential neighbouring ranks and checking if positions are      // all potential neighbouring ranks and checking if positions are
1485      // within bounds      // within bounds
1486        const dim_t numDOF=nDOF0*nDOF1*nDOF2;
1487        vector<IndexVector> colIndices(numDOF); // for the couple blocks
1488        RankVector neighbour;
1489        IndexVector offsetInShared(1,0);
1490        IndexVector sendShared, recvShared;
1491        int numShared=0;
1492      const int x=m_mpiInfo->rank%m_NX;      const int x=m_mpiInfo->rank%m_NX;
1493      const int y=m_mpiInfo->rank%(m_NX*m_NY)/m_NX;      const int y=m_mpiInfo->rank%(m_NX*m_NY)/m_NX;
1494      const int z=m_mpiInfo->rank/(m_NX*m_NY);      const int z=m_mpiInfo->rank/(m_NX*m_NY);
1495      for (int i2=-1; i2<2; i2++) {      for (int i2=-1; i2<2; i2++) {
1496          for (int i1=-1; i1<2; i1++) {          for (int i1=-1; i1<2; i1++) {
1497              for (int i0=-1; i0<2; i0++) {              for (int i0=-1; i0<2; i0++) {
1498                  // skip rank itself                  // skip this rank
1499                  if (i0==0 && i1==0 && i2==0)                  if (i0==0 && i1==0 && i2==0)
1500                      continue;                      continue;
1501                  // location of neighbour rank                  // location of neighbour rank
# Line 1397  Paso_SystemMatrixPattern* Brick::getPatt Line 1676  Paso_SystemMatrixPattern* Brick::getPatt
1676          }          }
1677      }      }
1678    
1679        // create connector
1680        Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
1681                numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
1682                &offsetInShared[0], 1, 0, m_mpiInfo);
1683        Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
1684                numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
1685                &offsetInShared[0], 1, 0, m_mpiInfo);
1686        m_connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
1687        Paso_SharedComponents_free(snd_shcomp);
1688        Paso_SharedComponents_free(rcv_shcomp);
1689    
1690        // create main and couple blocks
1691        Paso_Pattern *mainPattern = createMainPattern();
1692        Paso_Pattern *colPattern, *rowPattern;
1693        createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);
1694    
1695        // allocate paso distribution
1696        Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
1697                const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);
1698    
1699        // finally create the system matrix
1700        m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,
1701                distribution, distribution, mainPattern, colPattern, rowPattern,
1702                m_connector, m_connector);
1703    
1704        Paso_Distribution_free(distribution);
1705    
1706        // useful debug output
1707      /*      /*
1708      cout << "--- rcv_shcomp ---" << endl;      cout << "--- rcv_shcomp ---" << endl;
1709      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
# Line 1421  Paso_SystemMatrixPattern* Brick::getPatt Line 1728  Paso_SystemMatrixPattern* Brick::getPatt
1728      }      }
1729      */      */
1730    
     Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(  
             numDOF, neighbour.size(), &neighbour[0], &sendShared[0],  
             &offsetInShared[0], 1, 0, m_mpiInfo);  
     Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(  
             numDOF, neighbour.size(), &neighbour[0], &recvShared[0],  
             &offsetInShared[0], 1, 0, m_mpiInfo);  
     Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);  
     Paso_SharedComponents_free(snd_shcomp);  
     Paso_SharedComponents_free(rcv_shcomp);  
   
     // create patterns  
     dim_t M, N;  
     IndexVector ptr(1,0);  
     IndexVector index;  
   
     // main pattern  
     for (index_t i=0; i<numDOF; i++) {  
         // always add the node itself  
         index.push_back(i);  
         const int num=insertNeighbours(index, i);  
         ptr.push_back(ptr.back()+num+1);  
     }  
     M=N=ptr.size()-1;  
     // paso will manage the memory  
     index_t* indexC = MEMALLOC(index.size(),index_t);  
     index_t* ptrC = MEMALLOC(ptr.size(), index_t);  
     copy(index.begin(), index.end(), indexC);  
     copy(ptr.begin(), ptr.end(), ptrC);  
     Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,  
             M, N, ptrC, indexC);  
   
1731      /*      /*
1732      cout << "--- main_pattern ---" << endl;      cout << "--- main_pattern ---" << endl;
1733      cout << "M=" << M << ", N=" << N << endl;      cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
1734      for (size_t i=0; i<ptr.size(); i++) {      for (size_t i=0; i<mainPattern->numOutput+1; i++) {
1735          cout << "ptr[" << i << "]=" << ptr[i] << endl;          cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;
1736      }      }
1737      for (size_t i=0; i<index.size(); i++) {      for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {
1738          cout << "index[" << i << "]=" << index[i] << endl;          cout << "index[" << i << "]=" << mainPattern->index[i] << endl;
1739      }      }
1740      */      */
1741    
     // column & row couple patterns  
     ptr.assign(1, 0);  
     index.clear();  
   
     for (index_t i=0; i<numDOF; i++) {  
         index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());  
         ptr.push_back(ptr.back()+colIndices[i].size());  
     }  
   
     // paso will manage the memory  
     indexC = MEMALLOC(index.size(), index_t);  
     ptrC = MEMALLOC(ptr.size(), index_t);  
     copy(index.begin(), index.end(), indexC);  
     copy(ptr.begin(), ptr.end(), ptrC);  
     M=ptr.size()-1;  
     N=numShared;  
     Paso_Pattern *colCouplePattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,  
             M, N, ptrC, indexC);  
   
1742      /*      /*
1743      cout << "--- colCouple_pattern ---" << endl;      cout << "--- colCouple_pattern ---" << endl;
1744      cout << "M=" << M << ", N=" << N << endl;      cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
1745      for (size_t i=0; i<ptr.size(); i++) {      for (size_t i=0; i<colPattern->numOutput+1; i++) {
1746          cout << "ptr[" << i << "]=" << ptr[i] << endl;          cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
1747      }      }
1748      for (size_t i=0; i<index.size(); i++) {      for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
1749          cout << "index[" << i << "]=" << index[i] << endl;          cout << "index[" << i << "]=" << colPattern->index[i] << endl;
1750      }      }
1751      */      */
1752    
     // now build the row couple pattern  
     IndexVector ptr2(1,0);  
     IndexVector index2;  
     for (dim_t id=0; id<N; id++) {  
         dim_t n=0;  
         for (dim_t i=0; i<M; i++) {  
             for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {  
                 if (index[j]==id) {  
                     index2.push_back(i);  
                     n++;  
                     break;  
                 }  
             }  
         }  
         ptr2.push_back(ptr2.back()+n);  
     }  
   
     // paso will manage the memory  
     indexC = MEMALLOC(index2.size(), index_t);  
     ptrC = MEMALLOC(ptr2.size(), index_t);  
     copy(index2.begin(), index2.end(), indexC);  
     copy(ptr2.begin(), ptr2.end(), ptrC);  
     Paso_Pattern *rowCouplePattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,  
             N, M, ptrC, indexC);  
   
1753      /*      /*
1754      cout << "--- rowCouple_pattern ---" << endl;      cout << "--- rowCouple_pattern ---" << endl;
1755      cout << "M=" << N << ", N=" << M << endl;      cout << "M=" << rowPattern->numOutput << ", N=" << rowPattern->numInput << endl;
1756      for (size_t i=0; i<ptr2.size(); i++) {      for (size_t i=0; i<rowPattern->numOutput+1; i++) {
1757          cout << "ptr[" << i << "]=" << ptr2[i] << endl;          cout << "ptr[" << i << "]=" << rowPattern->ptr[i] << endl;
1758      }      }
1759      for (size_t i=0; i<index2.size(); i++) {      for (size_t i=0; i<rowPattern->ptr[rowPattern->numOutput]; i++) {
1760          cout << "index[" << i << "]=" << index2[i] << endl;          cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
1761      }      }
1762      */      */
1763    
     // allocate paso distribution  
     Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,  
             const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);  
   
     Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(  
             MATRIX_FORMAT_DEFAULT, distribution, distribution,  
             mainPattern, colCouplePattern, rowCouplePattern,  
             connector, connector);  
1764      Paso_Pattern_free(mainPattern);      Paso_Pattern_free(mainPattern);
1765      Paso_Pattern_free(colCouplePattern);      Paso_Pattern_free(colPattern);
1766      Paso_Pattern_free(rowCouplePattern);      Paso_Pattern_free(rowPattern);
     Paso_Distribution_free(distribution);  
     return pattern;  
 }  
   
 void Brick::Print_Mesh_Info(const bool full) const  
 {  
     RipleyDomain::Print_Mesh_Info(full);  
     if (full) {  
         cout << "     Id  Coordinates" << endl;  
         cout.precision(15);  
         cout.setf(ios::scientific, ios::floatfield);  
         pair<double,double> xdx = getFirstCoordAndSpacing(0);  
         pair<double,double> ydy = getFirstCoordAndSpacing(1);  
         pair<double,double> zdz = getFirstCoordAndSpacing(2);  
         for (index_t i=0; i < getNumNodes(); i++) {  
             cout << "  " << setw(5) << m_nodeId[i]  
                 << "  " << xdx.first+(i%m_N0)*xdx.second  
                 << "  " << ydy.first+(i%(m_N0*m_N1)/m_N0)*ydy.second  
                 << "  " << zdz.first+(i/(m_N0*m_N1))*zdz.second << endl;  
         }  
     }  
 }  
   
 IndexVector Brick::getNumNodesPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_N0);  
     ret.push_back(m_N1);  
     ret.push_back(m_N2);  
     return ret;  
 }  
   
 IndexVector Brick::getNumElementsPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_NE0);  
     ret.push_back(m_NE1);  
     ret.push_back(m_NE2);  
     return ret;  
 }  
   
 IndexVector Brick::getNumFacesPerBoundary() const  
 {  
     IndexVector ret(6, 0);  
     //left  
     if (m_offset0==0)  
         ret[0]=m_NE1*m_NE2;  
     //right  
     if (m_mpiInfo->rank%m_NX==m_NX-1)  
         ret[1]=m_NE1*m_NE2;  
     //bottom  
     if (m_offset1==0)  
         ret[2]=m_NE0*m_NE2;  
     //top  
     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)  
         ret[3]=m_NE0*m_NE2;  
     //front  
     if (m_offset2==0)  
         ret[4]=m_NE0*m_NE1;  
     //back  
     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)  
         ret[5]=m_NE0*m_NE1;  
     return ret;  
 }  
   
 pair<double,double> Brick::getFirstCoordAndSpacing(dim_t dim) const  
 {  
     if (dim==0)  
         return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);  
     else if (dim==1)  
         return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);  
     else if (dim==2)  
         return pair<double,double>((m_l2*m_offset2)/m_gNE2, m_l2/m_gNE2);  
   
     throw RipleyException("getFirstCoordAndSpacing(): invalid argument");  
 }  
   
 //protected  
 dim_t Brick::getNumDOF() const  
 {  
     return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY*(m_gNE2+1)/m_NZ;  
 }  
   
 //protected  
 dim_t Brick::getNumFaceElements() const  
 {  
     const IndexVector faces = getNumFacesPerBoundary();  
     dim_t n=0;  
     for (size_t i=0; i<faces.size(); i++)  
         n+=faces[i];  
     return n;  
 }  
   
 //protected  
 void Brick::assembleCoordinates(escript::Data& arg) const  
 {  
     escriptDataC x = arg.getDataC();  
     int numDim = m_numDim;  
     if (!isDataPointShapeEqual(&x, 1, &numDim))  
         throw RipleyException("setToX: Invalid Data object shape");  
     if (!numSamplesEqual(&x, 1, getNumNodes()))  
         throw RipleyException("setToX: Illegal number of samples in Data object");  
   
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
     pair<double,double> zdz = getFirstCoordAndSpacing(2);  
     arg.requireWrite();  
 #pragma omp parallel for  
     for (dim_t i2 = 0; i2 < m_N2; i2++) {  
         for (dim_t i1 = 0; i1 < m_N1; i1++) {  
             for (dim_t i0 = 0; i0 < m_N0; i0++) {  
                 double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);  
                 point[0] = xdx.first+i0*xdx.second;  
                 point[1] = ydy.first+i1*ydy.second;  
                 point[2] = zdz.first+i2*zdz.second;  
             }  
         }  
     }  
 }  
   
 //private  
 void Brick::populateSampleIds()  
 {  
     // identifiers are ordered from left to right, bottom to top, front to back  
     // globally  
   
     // build node distribution vector first.  
     // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes  
     m_nodeDistribution.assign(m_mpiInfo->size+1, 0);  
     const dim_t numDOF=getNumDOF();  
     for (dim_t k=1; k<m_mpiInfo->size; k++) {  
         m_nodeDistribution[k]=k*numDOF;  
     }  
     m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();  
     m_nodeId.resize(getNumNodes());  
     m_dofId.resize(numDOF);  
     m_elementId.resize(getNumElements());  
     m_faceId.resize(getNumFaceElements());  
   
 #pragma omp parallel  
     {  
 #pragma omp for nowait  
         // nodes  
         for (dim_t i2=0; i2<m_N2; i2++) {  
             for (dim_t i1=0; i1<m_N1; i1++) {  
                 for (dim_t i0=0; i0<m_N0; i0++) {  
                     m_nodeId[i0+i1*m_N0+i2*m_N0*m_N1] =  
                         (m_offset2+i2)*(m_gNE0+1)*(m_gNE1+1)  
                         +(m_offset1+i1)*(m_gNE0+1)  
                         +m_offset0+i0;  
                 }  
             }  
         }  
   
         // degrees of freedom  
 #pragma omp for nowait  
         for (dim_t k=0; k<numDOF; k++)  
             m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;  
   
         // elements  
 #pragma omp for nowait  
         for (dim_t i2=0; i2<m_NE2; i2++) {  
             for (dim_t i1=0; i1<m_NE1; i1++) {  
                 for (dim_t i0=0; i0<m_NE0; i0++) {  
                     m_elementId[i0+i1*m_NE0+i2*m_NE0*m_NE1] =  
                         (m_offset2+i2)*m_gNE0*m_gNE1  
                         +(m_offset1+i1)*m_gNE0  
                         +m_offset0+i0;  
                 }  
             }  
         }  
   
         // face elements  
 #pragma omp for  
         for (dim_t k=0; k<getNumFaceElements(); k++)  
             m_faceId[k]=k;  
     } // end parallel section  
   
     m_nodeTags.assign(getNumNodes(), 0);  
     updateTagsInUse(Nodes);  
   
     m_elementTags.assign(getNumElements(), 0);  
     updateTagsInUse(Elements);  
   
     // generate face offset vector and set face tags  
     const IndexVector facesPerEdge = getNumFacesPerBoundary();  
     const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20, FRONT=100, BACK=200;  
     const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP, FRONT, BACK };  
     m_faceOffset.assign(facesPerEdge.size(), -1);  
     m_faceTags.clear();  
     index_t offset=0;  
     for (size_t i=0; i<facesPerEdge.size(); i++) {  
         if (facesPerEdge[i]>0) {  
             m_faceOffset[i]=offset;  
             offset+=facesPerEdge[i];  
             m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);  
         }  
     }  
     setTagMap("left", LEFT);  
     setTagMap("right", RIGHT);  
     setTagMap("bottom", BOTTOM);  
     setTagMap("top", TOP);  
     setTagMap("front", FRONT);  
     setTagMap("back", BACK);  
     updateTagsInUse(FaceElements);  
 }  
   
 //private  
 int Brick::insertNeighbours(IndexVector& index, index_t node) const  
 {  
     const index_t nDOF0 = (m_gNE0+1)/m_NX;  
     const index_t nDOF1 = (m_gNE1+1)/m_NY;  
     const index_t nDOF2 = (m_gNE2+1)/m_NZ;  
     const int x=node%nDOF0;  
     const int y=node%(nDOF0*nDOF1)/nDOF0;  
     const int z=node/(nDOF0*nDOF1);  
     int num=0;  
     // loop through potential neighbours and add to index if positions are  
     // within bounds  
     for (int i2=-1; i2<2; i2++) {  
         for (int i1=-1; i1<2; i1++) {  
             for (int i0=-1; i0<2; i0++) {  
                 // skip node itself  
                 if (i0==0 && i1==0 && i2==0)  
                     continue;  
                 // location of neighbour node  
                 const int nx=x+i0;  
                 const int ny=y+i1;  
                 const int nz=z+i2;  
                 if (nx>=0 && ny>=0 && nz>=0  
                         && nx<nDOF0 && ny<nDOF1 && nz<nDOF2) {  
                     index.push_back(nz*nDOF0*nDOF1+ny*nDOF0+nx);  
                     num++;  
                 }  
             }  
         }  
     }  
   
     return num;  
 }  
   
 //protected  
 void Brick::addToSystemMatrix(Paso_SystemMatrix* mat,  
        const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,  
        dim_t num_Sol, const vector<double>& array) const  
 {  
     const dim_t numMyCols = mat->pattern->mainPattern->numInput;  
     const dim_t numMyRows = mat->pattern->mainPattern->numOutput;  
   
     const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;  
     const index_t* mainBlock_index = mat->mainBlock->pattern->index;  
     double* mainBlock_val = mat->mainBlock->val;  
     const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;  
     const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;  
     double* col_coupleBlock_val = mat->col_coupleBlock->val;  
     const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;  
     const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;  
     double* row_coupleBlock_val = mat->row_coupleBlock->val;  
   
     for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {  
         // down columns of array  
         const dim_t j_Eq = nodes_Eq[k_Eq];  
         const dim_t i_row = j_Eq;  
 //printf("row:%d\n", i_row);  
         // only look at the matrix rows stored on this processor  
         if (i_row < numMyRows) {  
             for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {  
                 const dim_t i_col = nodes_Sol[k_Sol];  
                 if (i_col < numMyCols) {  
                     for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {  
                         if (mainBlock_index[k] == i_col) {  
                             mainBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];  
                             break;  
                         }  
                     }  
                 } else {  
                     for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {  
 //cout << "col:" << i_col-numMyCols << " colIdx:" << col_coupleBlock_index[k] << endl;  
                         if (col_coupleBlock_index[k] == i_col - numMyCols) {  
                             col_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];  
                             break;  
                         }  
                     }  
                 }  
             }  
         } else {  
             for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {  
                 // across rows of array  
                 const dim_t i_col = nodes_Sol[k_Sol];  
                 if (i_col < numMyCols) {  
                     for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];  
                          k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)  
                     {  
 //cout << "col:" << i_col << " rowIdx:" << row_coupleBlock_index[k] << endl;  
                         if (row_coupleBlock_index[k] == i_col) {  
                             row_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];  
                             break;  
                         }  
                     }  
                 }  
             }  
         }  
     }  
1767  }  }
1768    
1769  //protected  //protected
# Line 2131  void Brick::interpolateNodesOnFaces(escr Line 2052  void Brick::interpolateNodesOnFaces(escr
2052      }      }
2053  }  }
2054    
 //protected  
 void Brick::nodesToDOF(escript::Data& out, escript::Data& in) const  
 {  
     const dim_t numComp = in.getDataPointSize();  
     out.requireWrite();  
   
     const index_t left = (m_offset0==0 ? 0 : 1);  
     const index_t bottom = (m_offset1==0 ? 0 : 1);  
     const index_t front = (m_offset2==0 ? 0 : 1);  
     const index_t nDOF0 = (m_gNE0+1)/m_NX;  
     const index_t nDOF1 = (m_gNE1+1)/m_NY;  
     const index_t nDOF2 = (m_gNE2+1)/m_NZ;  
 #pragma omp parallel for  
     for (index_t i=0; i<nDOF2; i++) {  
         for (index_t j=0; j<nDOF1; j++) {  
             for (index_t k=0; k<nDOF0; k++) {  
                 const index_t n=k+left+(j+bottom)*m_N0+(i+front)*m_N0*m_N1;  
                 const double* src=in.getSampleDataRO(n);  
                 copy(src, src+numComp, out.getSampleDataRW(k+j*nDOF0+i*nDOF0*nDOF1));  
             }  
         }  
     }  
 }  
   
 //protected  
 void Brick::dofToNodes(escript::Data& out, escript::Data& in) const  
 {  
     const dim_t numComp = in.getDataPointSize();  
     out.requireWrite();  
   
     //TODO: use coupler to get the rest of the values  
   
     const index_t left = (m_offset0==0 ? 0 : 1);  
     const index_t bottom = (m_offset1==0 ? 0 : 1);  
     const index_t front = (m_offset2==0 ? 0 : 1);  
     const index_t nDOF0 = (m_gNE0+1)/m_NX;  
     const index_t nDOF1 = (m_gNE1+1)/m_NY;  
     const index_t nDOF2 = (m_gNE2+1)/m_NZ;  
 #pragma omp parallel for  
     for (index_t i=0; i<nDOF2; i++) {  
         for (index_t j=0; j<nDOF1; j++) {  
             for (index_t k=0; k<nDOF0; k++) {  
                 const double* src=in.getSampleDataRO(k+j*nDOF0+i*nDOF0*nDOF1);  
                 const index_t n=k+left+(j+bottom)*m_N0+(i+front)*m_N0*m_N1;  
                 copy(src, src+numComp, out.getSampleDataRW(n));  
             }  
         }  
     }  
 }  
   
2055  //protected  //protected
2056  void Brick::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,  void Brick::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
2057          const escript::Data& A, const escript::Data& B,          const escript::Data& A, const escript::Data& B,

Legend:
Removed from v.3755  
changed lines
  Added in v.3756

  ViewVC Help
Powered by ViewVC 1.1.26