/[escript]/trunk/finley/src/NodeFile.cpp
ViewVC logotype

Diff of /trunk/finley/src/NodeFile.cpp

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

revision 4495 by caltinay, Wed Jul 3 06:32:11 2013 UTC revision 4496 by caltinay, Mon Jul 15 06:53:44 2013 UTC
# Line 54  static void scatterEntries(int n, int* i Line 54  static void scatterEntries(int n, int* i
54  }  }
55    
56  // helper function  // helper function
57  static void gatherEntries(int n, int* index, int min_index, int max_index,  static void gatherEntries(int n, const int* index, int min_index, int max_index,
58                            int* Id_out, int* Id_in, int* Tag_out, int* Tag_in,                            int* Id_out, int* Id_in, int* Tag_out, int* Tag_in,
59                            int* globalDegreesOfFreedom_out,                            int* globalDegreesOfFreedom_out,
60                            int* globalDegreesOfFreedom_in,                            int* globalDegreesOfFreedom_in,
# Line 207  void NodeFile::setCoordinates(const escr Line 207  void NodeFile::setCoordinates(const escr
207          ss << "NodeFile::setCoordinates: number of dimensions of new "          ss << "NodeFile::setCoordinates: number of dimensions of new "
208              "coordinates has to be " << numDim;              "coordinates has to be " << numDim;
209          const std::string errorMsg(ss.str());          const std::string errorMsg(ss.str());
210          Finley_setError(VALUE_ERROR, errorMsg.c_str());          setError(VALUE_ERROR, errorMsg.c_str());
211      } else if (cNewX.getNumDataPointsPerSample() != 1 ||      } else if (cNewX.getNumDataPointsPerSample() != 1 ||
212              cNewX.getNumSamples() != numNodes) {              cNewX.getNumSamples() != numNodes) {
213          std::stringstream ss;          std::stringstream ss;
214          ss << "NodeFile::setCoordinates: number of given nodes must be "          ss << "NodeFile::setCoordinates: number of given nodes must be "
215              << numNodes;              << numNodes;
216          const std::string errorMsg(ss.str());          const std::string errorMsg(ss.str());
217          Finley_setError(VALUE_ERROR, errorMsg.c_str());          setError(VALUE_ERROR, errorMsg.c_str());
218      } else {      } else {
219          const size_t numDim_size=numDim*sizeof(double);          const size_t numDim_size=numDim*sizeof(double);
220          Finley_increaseStatus(this);          ++status;
221          escript::Data& newX = *const_cast<escript::Data*>(&cNewX);          escript::Data& newX = *const_cast<escript::Data*>(&cNewX);
222  #pragma omp parallel for  #pragma omp parallel for
223          for (int n=0; n<numNodes; n++) {          for (int n=0; n<numNodes; n++) {
# Line 229  void NodeFile::setCoordinates(const escr Line 229  void NodeFile::setCoordinates(const escr
229  /// sets tags to newTag where mask>0  /// sets tags to newTag where mask>0
230  void NodeFile::setTags(const int newTag, const escript::Data& cMask)  void NodeFile::setTags(const int newTag, const escript::Data& cMask)
231  {  {
232      Finley_resetError();      resetError();
233    
234      if (1 != cMask.getDataPointSize()) {      if (1 != cMask.getDataPointSize()) {
235         Finley_setError(TYPE_ERROR, "NodeFile::setTags: number of components of mask must be 1.");         setError(TYPE_ERROR, "NodeFile::setTags: number of components of mask must be 1.");
236         return;         return;
237      } else if (cMask.getNumDataPointsPerSample() != 1 ||      } else if (cMask.getNumDataPointsPerSample() != 1 ||
238              cMask.getNumSamples() != numNodes) {              cMask.getNumSamples() != numNodes) {
239         Finley_setError(TYPE_ERROR, "NodeFile::setTags: illegal number of samples of mask Data object");         setError(TYPE_ERROR, "NodeFile::setTags: illegal number of samples of mask Data object");
240         return;         return;
241      }      }
242    
# Line 320  void NodeFile::copyTable(int offset, int Line 320  void NodeFile::copyTable(int offset, int
320  {  {
321      // check number of dimensions and table size      // check number of dimensions and table size
322      if (numDim != in->numDim) {      if (numDim != in->numDim) {
323          Finley_setError(TYPE_ERROR, "NodeFile::copyTable: dimensions of node files don't match");          setError(TYPE_ERROR, "NodeFile::copyTable: dimensions of node files don't match");
324          return;          return;
325      }      }
326      if (numNodes < in->numNodes+offset) {      if (numNodes < in->numNodes+offset) {
327          Finley_setError(MEMORY_ERROR, "NodeFile::copyTable: node table is too small.");          setError(MEMORY_ERROR, "NodeFile::copyTable: node table is too small.");
328          return;          return;
329      }      }
330    
# Line 360  void NodeFile::gather(int* index, const Line 360  void NodeFile::gather(int* index, const
360              numDim, Coordinates, in->Coordinates);              numDim, Coordinates, in->Coordinates);
361  }  }
362    
363  void NodeFile::gather_global(int* index, const NodeFile* in)  void NodeFile::gather_global(const std::vector<int>& index, const NodeFile* in)
364  {  {
365      // get the global range of node ids      // get the global range of node ids
366      const std::pair<int,int> id_range(in->getGlobalIdRange());      const std::pair<int,int> id_range(in->getGlobalIdRange());
# Line 423  void NodeFile::gather_global(int* index, Line 423  void NodeFile::gather_global(int* index,
423  #endif  #endif
424      buffer_rank=in->MPIInfo->rank;      buffer_rank=in->MPIInfo->rank;
425      for (int p=0; p<in->MPIInfo->size; ++p) {      for (int p=0; p<in->MPIInfo->size; ++p) {
426          gatherEntries(numNodes, index, distribution[buffer_rank],          gatherEntries(numNodes, &index[0], distribution[buffer_rank],
427                  distribution[buffer_rank+1], Id, Id_buffer, Tag, Tag_buffer,                  distribution[buffer_rank+1], Id, Id_buffer, Tag, Tag_buffer,
428                  globalDegreesOfFreedom, globalDegreesOfFreedom_buffer, numDim,                  globalDegreesOfFreedom, globalDegreesOfFreedom_buffer, numDim,
429                  Coordinates, Coordinates_buffer);                  Coordinates, Coordinates_buffer);
# Line 454  void NodeFile::gather_global(int* index, Line 454  void NodeFile::gather_global(int* index,
454              ss << "NodeFile::gather_global: Node id " << Id[n]              ss << "NodeFile::gather_global: Node id " << Id[n]
455                  << " at position " << n << " is referenced but not defined.";                  << " at position " << n << " is referenced but not defined.";
456              const std::string errorMsg(ss.str());              const std::string errorMsg(ss.str());
457              Finley_setError(VALUE_ERROR, errorMsg.c_str());              setError(VALUE_ERROR, errorMsg.c_str());
458          }          }
459      }      }
460      delete[] Id_buffer;      delete[] Id_buffer;
# Line 465  void NodeFile::gather_global(int* index, Line 465  void NodeFile::gather_global(int* index,
465      Esys_MPIInfo_noError(in->MPIInfo);      Esys_MPIInfo_noError(in->MPIInfo);
466  }  }
467    
468  void NodeFile::assignMPIRankToDOFs(Esys_MPI_rank* mpiRankOfDOF,  void NodeFile::assignMPIRankToDOFs(std::vector<int>& mpiRankOfDOF,
469                                     int *distribution)                                     const std::vector<int>& distribution)
470  {  {
471      Esys_MPI_rank p_min=MPIInfo->size, p_max=-1;      Esys_MPI_rank p_min=MPIInfo->size, p_max=-1;
472      // first we retrieve the min and max DOF on this processor to reduce      // first we retrieve the min and max DOF on this processor to reduce
# Line 489  void NodeFile::assignMPIRankToDOFs(Esys_ Line 489  void NodeFile::assignMPIRankToDOFs(Esys_
489      }      }
490  }  }
491    
492  int NodeFile::prepareLabeling(int* mask, std::vector<int>& buffer,  int NodeFile::prepareLabeling(const std::vector<short>& mask,
493                                  std::vector<int>& buffer,
494                                std::vector<int>& distribution, bool useNodes)                                std::vector<int>& distribution, bool useNodes)
495  {  {
496      const int UNSET_ID=-1,SET_ID=1;      const int UNSET_ID=-1,SET_ID=1;
# Line 528  int NodeFile::prepareLabeling(int* mask, Line 529  int NodeFile::prepareLabeling(int* mask,
529          const int id1=distribution[buffer_rank+1];          const int id1=distribution[buffer_rank+1];
530  #pragma omp parallel for  #pragma omp parallel for
531          for (int n=0; n<numNodes; n++) {          for (int n=0; n<numNodes; n++) {
532              if (!mask || mask[n]>-1) {              if (mask.size()<numNodes || mask[n]>-1) {
533                  const int k=indexArray[n];                  const int k=indexArray[n];
534                  if (id0<=k && k<id1) {                  if (id0<=k && k<id1) {
535                      buffer[k-id0] = SET_ID;                      buffer[k-id0] = SET_ID;
# Line 557  int NodeFile::createDenseDOFLabeling() Line 558  int NodeFile::createDenseDOFLabeling()
558      int new_numGlobalDOFs=0;      int new_numGlobalDOFs=0;
559    
560      // retrieve the number of own DOFs and fill buffer      // retrieve the number of own DOFs and fill buffer
561      loc_offsets[MPIInfo->rank]=prepareLabeling(NULL, DOF_buffer, distribution,      loc_offsets[MPIInfo->rank]=prepareLabeling(std::vector<short>(),
562                                                 false);              DOF_buffer, distribution, false);
563  #ifdef ESYS_MPI  #ifdef ESYS_MPI
564      MPI_Allreduce(&loc_offsets[0], &offsets[0], MPIInfo->size, MPI_INT,      MPI_Allreduce(&loc_offsets[0], &offsets[0], MPIInfo->size, MPI_INT,
565                    MPI_SUM, MPIInfo->comm);                    MPI_SUM, MPIInfo->comm);
# Line 611  int NodeFile::createDenseDOFLabeling() Line 612  int NodeFile::createDenseDOFLabeling()
612      return new_numGlobalDOFs;      return new_numGlobalDOFs;
613  }  }
614    
615  int NodeFile::createDenseNodeLabeling(int* node_distribution,  int NodeFile::createDenseNodeLabeling(std::vector<int>& nodeDistribution,
616                                        const int* dof_distribution)                                        const std::vector<int>& dofDistribution)
617  {  {
618      const int UNSET_ID=-1, SET_ID=1;      const int UNSET_ID=-1, SET_ID=1;
619      const int myFirstDOF=dof_distribution[MPIInfo->rank];      const int myFirstDOF=dofDistribution[MPIInfo->rank];
620      const int myLastDOF=dof_distribution[MPIInfo->rank+1];      const int myLastDOF=dofDistribution[MPIInfo->rank+1];
621    
622      // find the range of node ids controlled by me      // find the range of node ids controlled by me
623      int min_id=std::numeric_limits<int>::max();      int min_id=std::numeric_limits<int>::max();
# Line 672  int NodeFile::createDenseNodeLabeling(in Line 673  int NodeFile::createDenseNodeLabeling(in
673      }      }
674      // make the local number of nodes globally available      // make the local number of nodes globally available
675  #ifdef ESYS_MPI  #ifdef ESYS_MPI
676      MPI_Allgather(&myNewNumNodes, 1, MPI_INT, node_distribution, 1, MPI_INT,      MPI_Allgather(&myNewNumNodes, 1, MPI_INT, &nodeDistribution[0], 1, MPI_INT,
677                    MPIInfo->comm);                    MPIInfo->comm);
678  #else  #else
679      node_distribution[0]=myNewNumNodes;      nodeDistribution[0]=myNewNumNodes;
680  #endif  #endif
681    
682      int globalNumNodes=0;      int globalNumNodes=0;
683      for (int p=0; p<MPIInfo->size; ++p) {      for (int p=0; p<MPIInfo->size; ++p) {
684          const int itmp=node_distribution[p];          const int itmp=nodeDistribution[p];
685          node_distribution[p]=globalNumNodes;          nodeDistribution[p]=globalNumNodes;
686          globalNumNodes+=itmp;          globalNumNodes+=itmp;
687      }      }
688      node_distribution[MPIInfo->size]=globalNumNodes;      nodeDistribution[MPIInfo->size]=globalNumNodes;
689    
690      // offset node buffer      // offset node buffer
691  #pragma omp parallel for  #pragma omp parallel for
692      for (int n=0; n<my_buffer_len; n++)      for (int n=0; n<my_buffer_len; n++)
693          Node_buffer[n+header_len]+=node_distribution[MPIInfo->rank];          Node_buffer[n+header_len]+=nodeDistribution[MPIInfo->rank];
694    
695      // now we send this buffer around to assign global node index      // now we send this buffer around to assign global node index
696  #ifdef ESYS_MPI  #ifdef ESYS_MPI
# Line 700  int NodeFile::createDenseNodeLabeling(in Line 701  int NodeFile::createDenseNodeLabeling(in
701      for (int p=0; p<MPIInfo->size; ++p) {      for (int p=0; p<MPIInfo->size; ++p) {
702          const int nodeID_0=Node_buffer[0];          const int nodeID_0=Node_buffer[0];
703          const int nodeID_1=Node_buffer[1];          const int nodeID_1=Node_buffer[1];
704          const int dof0=dof_distribution[buffer_rank];          const int dof0=dofDistribution[buffer_rank];
705          const int dof1=dof_distribution[buffer_rank+1];          const int dof1=dofDistribution[buffer_rank+1];
706          if (nodeID_0 <= nodeID_1) {          if (nodeID_0 <= nodeID_1) {
707  #pragma omp parallel for  #pragma omp parallel for
708              for (int n=0; n<numNodes; n++) {              for (int n=0; n<numNodes; n++) {
# Line 725  int NodeFile::createDenseNodeLabeling(in Line 726  int NodeFile::createDenseNodeLabeling(in
726      return globalNumNodes;      return globalNumNodes;
727  }  }
728    
729  int NodeFile::createDenseReducedLabeling(int* reducedMask, bool useNodes)  int NodeFile::createDenseReducedLabeling(const std::vector<short>& reducedMask,
730                                             bool useNodes)
731  {  {
732      std::vector<int> buffer;      std::vector<int> buffer;
733      std::vector<int> distribution;      std::vector<int> distribution;
# Line 832  void NodeFile::createDOFMappingAndCoupli Line 834  void NodeFile::createDOFMappingAndCoupli
834      }      }
835    
836      if (!((min_DOF<=myFirstDOF) && (myLastDOF-1<=max_DOF))) {      if (!((min_DOF<=myFirstDOF) && (myLastDOF-1<=max_DOF))) {
837          Finley_setError(SYSTEM_ERROR, "Local elements do not span local degrees of freedom.");          setError(SYSTEM_ERROR, "Local elements do not span local degrees of freedom.");
838          return;          return;
839      }      }
840      const int UNUSED = -1;      const int UNUSED = -1;
# Line 934  void NodeFile::createDOFMappingAndCoupli Line 936  void NodeFile::createDOFMappingAndCoupli
936    
937      // define how to get DOF values for controlled but other processors      // define how to get DOF values for controlled but other processors
938  #ifdef BOUNDS_CHECK  #ifdef BOUNDS_CHECK
939      if (offsetInShared[numNeighbours] >= numNodes*(p_max-p_min+1)) {      if (offsetInShared[numNeighbors] >= numNodes*(p_max-p_min+1)) {
940          printf("BOUNDS_CHECK %s %d\n", __FILE__, __LINE__);          printf("BOUNDS_CHECK %s %d\n", __FILE__, __LINE__);
941          exit(1);          exit(1);
942      }      }
# Line 1000  void NodeFile::createDOFMappingAndCoupli Line 1002  void NodeFile::createDOFMappingAndCoupli
1002              myLastDOF-myFirstDOF, numNeighbors, &neighbor[0], &shared[0],              myLastDOF-myFirstDOF, numNeighbors, &neighbor[0], &shared[0],
1003              &offsetInShared[0], 1, 0, MPIInfo);              &offsetInShared[0], 1, 0, MPIInfo);
1004    
1005      if (Finley_noError()) {      if (noError()) {
1006          if (use_reduced_elements) {          if (use_reduced_elements) {
1007              reducedDegreesOfFreedomConnector=Paso_Connector_alloc(snd_shcomp, rcv_shcomp);              reducedDegreesOfFreedomConnector=Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
1008          } else {          } else {
# Line 1012  void NodeFile::createDOFMappingAndCoupli Line 1014  void NodeFile::createDOFMappingAndCoupli
1014      Paso_SharedComponents_free(snd_shcomp);      Paso_SharedComponents_free(snd_shcomp);
1015  }  }
1016    
1017  void NodeFile::createNodeMappings(int numReducedNodes,  void NodeFile::createNodeMappings(const std::vector<int>& indexReducedNodes,
1018                                    const std::vector<int>& indexReducedNodes,                                    const std::vector<int>& dofDist,
1019                                    const int* dofDist, const int* nodeDist)                                    const std::vector<int>& nodeDist)
1020  {  {
1021      const int mpiSize=MPIInfo->size;      const int mpiSize=MPIInfo->size;
1022      const int myRank=MPIInfo->rank;      const int myRank=MPIInfo->rank;
# Line 1027  void NodeFile::createNodeMappings(int nu Line 1029  void NodeFile::createNodeMappings(int nu
1029      const int myLastNode=nodeDist[myRank+1];      const int myLastNode=nodeDist[myRank+1];
1030      const int myNumNodes=myLastNode-myFirstNode;      const int myNumNodes=myLastNode-myFirstNode;
1031    
1032      std::vector<int> maskMyReducedDOF(myNumDOF, -1);      std::vector<short> maskMyReducedDOF(myNumDOF, -1);
1033      std::vector<int> maskMyReducedNodes(myNumNodes, -1);      std::vector<short> maskMyReducedNodes(myNumNodes, -1);
1034    
1035      // mark the nodes used by the reduced mesh      // mark the nodes used by the reduced mesh
1036  #pragma omp parallel for  #pragma omp parallel for
1037      for (int i=0; i<numReducedNodes; ++i) {      for (int i=0; i<indexReducedNodes.size(); ++i) {
1038          int k=globalNodesIndex[indexReducedNodes[i]];          int k=globalNodesIndex[indexReducedNodes[i]];
1039          if (k>=myFirstNode && myLastNode>k)          if (k>=myFirstNode && myLastNode>k)
1040              maskMyReducedNodes[k-myFirstNode]=i;              maskMyReducedNodes[k-myFirstNode]=1;
1041          k=globalDegreesOfFreedom[indexReducedNodes[i]];          k=globalDegreesOfFreedom[indexReducedNodes[i]];
1042          if (k>=myFirstDOF && myLastDOF>k) {          if (k>=myFirstDOF && myLastDOF>k) {
1043              maskMyReducedDOF[k-myFirstDOF]=i;              maskMyReducedDOF[k-myFirstDOF]=1;
1044          }          }
1045      }      }
1046      std::vector<int> indexMyReducedDOF(myNumDOF);      std::vector<int> indexMyReducedDOF = util::packMask(maskMyReducedDOF);
1047      std::vector<int> indexMyReducedNodes(myNumNodes);      int myNumReducedDOF=indexMyReducedDOF.size();
1048      int myNumReducedDOF=util::packMask(myNumDOF, &maskMyReducedDOF[0], &indexMyReducedDOF[0]);      std::vector<int> indexMyReducedNodes = util::packMask(maskMyReducedNodes);
1049      int myNumReducedNodes=util::packMask(myNumNodes, &maskMyReducedNodes[0], &indexMyReducedNodes[0]);      int myNumReducedNodes=indexMyReducedNodes.size();
1050    
1051      std::vector<int> rdofDist(mpiSize+1);      std::vector<int> rdofDist(mpiSize+1);
1052      std::vector<int> rnodeDist(mpiSize+1);      std::vector<int> rnodeDist(mpiSize+1);
# Line 1070  void NodeFile::createNodeMappings(int nu Line 1072  void NodeFile::createNodeMappings(int nu
1072      rdofDist[mpiSize]=globalNumReducedDOF;      rdofDist[mpiSize]=globalNumReducedDOF;
1073    
1074      // ==== distribution of Nodes ===============================      // ==== distribution of Nodes ===============================
1075      nodesDistribution=Paso_Distribution_alloc(MPIInfo, nodeDist, 1, 0);      nodesDistribution=Paso_Distribution_alloc(MPIInfo, &nodeDist[0], 1, 0);
1076      // ==== distribution of DOFs ================================      // ==== distribution of DOFs ================================
1077      degreesOfFreedomDistribution=Paso_Distribution_alloc(MPIInfo, &dofDist[0], 1,0);      degreesOfFreedomDistribution=Paso_Distribution_alloc(MPIInfo, &dofDist[0], 1,0);
1078      // ==== distribution of reduced Nodes =======================      // ==== distribution of reduced Nodes =======================
# Line 1080  void NodeFile::createNodeMappings(int nu Line 1082  void NodeFile::createNodeMappings(int nu
1082    
1083      std::vector<int> nodeMask(numNodes);      std::vector<int> nodeMask(numNodes);
1084    
1085      if (Finley_noError()) {      if (noError()) {
1086          const int UNUSED = -1;          const int UNUSED = -1;
1087          // ==== nodes mapping which is a dummy structure ========          // ==== nodes mapping which is a dummy structure ========
1088  #pragma omp parallel for  #pragma omp parallel for
# Line 1093  void NodeFile::createNodeMappings(int nu Line 1095  void NodeFile::createNodeMappings(int nu
1095          for (int i=0; i<numNodes; ++i)          for (int i=0; i<numNodes; ++i)
1096              nodeMask[i]=UNUSED;              nodeMask[i]=UNUSED;
1097  #pragma omp parallel for  #pragma omp parallel for
1098          for (int i=0; i<numReducedNodes; ++i)          for (int i=0; i<indexReducedNodes.size(); ++i)
1099              nodeMask[indexReducedNodes[i]]=i;              nodeMask[indexReducedNodes[i]]=i;
1100          reducedNodesMapping.assign(nodeMask, UNUSED);          reducedNodesMapping.assign(nodeMask, UNUSED);
1101      }      }
1102      // ==== mapping between nodes and DOFs + DOF connector      // ==== mapping between nodes and DOFs + DOF connector
1103      if (Finley_noError())      if (noError())
1104          createDOFMappingAndCoupling(false);          createDOFMappingAndCoupling(false);
1105      // ==== mapping between nodes and reduced DOFs + reduced DOF connector      // ==== mapping between nodes and reduced DOFs + reduced DOF connector
1106      if (Finley_noError())      if (noError())
1107          createDOFMappingAndCoupling(true);          createDOFMappingAndCoupling(true);
1108    
1109      // get the Ids for DOFs and reduced nodes      // get the Ids for DOFs and reduced nodes
1110      if (Finley_noError()) {      if (noError()) {
1111  #pragma omp parallel  #pragma omp parallel
1112          {          {
1113  #pragma omp for  #pragma omp for

Legend:
Removed from v.4495  
changed lines
  Added in v.4496

  ViewVC Help
Powered by ViewVC 1.1.26