/[escript]/branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
ViewVC logotype

Diff of /branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp

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

revision 3733 by caltinay, Fri Dec 9 04:02:56 2011 UTC revision 3740 by caltinay, Tue Dec 13 01:37:21 2011 UTC
# Line 108  pair<int,int> RipleyDomain::getDataShape Line 108  pair<int,int> RipleyDomain::getDataShape
108      }      }
109    
110      stringstream msg;      stringstream msg;
111      msg << "getDataShape(): Unsupported function space type "      msg << "getDataShape(): Unsupported function space type " << fsType
112          << functionSpaceTypeAsString(fsType) << " for " << getDescription();          << " for " << getDescription();
113      throw RipleyException(msg.str());      throw RipleyException(msg.str());
114  }  }
115    
# Line 129  bool RipleyDomain::commonFunctionSpace(c Line 129  bool RipleyDomain::commonFunctionSpace(c
129     /*     /*
130      The idea is to use equivalence classes (i.e. types which can be      The idea is to use equivalence classes (i.e. types which can be
131      interpolated back and forth):      interpolated back and forth):
132      class 1: DOF <-> Nodes      class 0: DOF <-> Nodes
133      class 2: ReducedDOF <-> ReducedNodes      class 1: ReducedDOF <-> ReducedNodes
134      class 3: Points      class 2: Points
135      class 4: Elements      class 3: Elements
136      class 5: ReducedElements      class 4: ReducedElements
137      class 6: FaceElements      class 5: FaceElements
138      class 7: ReducedFaceElements      class 6: ReducedFaceElements
     class 8: ContactElementZero <-> ContactElementOne  
     class 9: ReducedContactElementZero <-> ReducedContactElementOne  
139    
140      There is also a set of lines. Interpolation is possible down a line but not      There is also a set of lines. Interpolation is possible down a line but not
141      between lines.      between lines.
142      class 1 and 2 belong to all lines so aren't considered.      class 0 and 1 belong to all lines so aren't considered.
143      line 0: class 3      line 0: class 2
144      line 1: class 4,5      line 1: class 3,4
145      line 2: class 6,7      line 2: class 5,6
     line 3: class 8,9  
146    
147      For classes with multiple members (eg class 2) we have vars to record if      For classes with multiple members (eg class 1) we have vars to record if
148      there is at least one instance. e.g. hasnodes is true if we have at least      there is at least one instance. e.g. hasnodes is true if we have at least
149      one instance of Nodes.      one instance of Nodes.
150      */      */
151      if (fs.empty())      if (fs.empty())
152          return false;          return false;
153      vector<int> hasclass(10);      vector<bool> hasclass(7, false);
154      vector<int> hasline(4);      vector<int> hasline(3, 0);
155      bool hasnodes=false;      bool hasnodes=false;
156      bool hasrednodes=false;      bool hasrednodes=false;
157      for (int i=0;i<fs.size();++i) {      for (size_t i=0; i<fs.size(); ++i) {
158          switch (fs[i]) {          switch (fs[i]) {
159              case Nodes: hasnodes=true; // no break is deliberate              case Nodes: hasnodes=true; // fall through
160              case DegreesOfFreedom:              case DegreesOfFreedom:
161                  hasclass[1]=1;                  hasclass[0]=true;
162                  break;                  break;
163              case ReducedNodes: hasrednodes=true; // no break is deliberate              case ReducedNodes: hasrednodes=true; // fall through
164              case ReducedDegreesOfFreedom:              case ReducedDegreesOfFreedom:
165                  hasclass[2]=1;                  hasclass[1]=true;
166                  break;                  break;
167              case Points:              case Points:
168                  hasline[0]=1;                  hasline[0]=1;
169                  hasclass[3]=1;                  hasclass[2]=true;
170                  break;                  break;
171              case Elements:              case Elements:
172                  hasclass[4]=1;                  hasclass[3]=true;
173                  hasline[1]=1;                  hasline[1]=1;
174                  break;                  break;
175              case ReducedElements:              case ReducedElements:
176                  hasclass[5]=1;                  hasclass[4]=true;
177                  hasline[1]=1;                  hasline[1]=1;
178                  break;                  break;
179              case FaceElements:              case FaceElements:
180                  hasclass[6]=1;                  hasclass[5]=true;
181                  hasline[2]=1;                  hasline[2]=1;
182                  break;                  break;
183              case ReducedFaceElements:              case ReducedFaceElements:
184                  hasclass[7]=1;                  hasclass[6]=true;
185                  hasline[2]=1;                  hasline[2]=1;
186                  break;                  break;
187              default:              default:
188                  return false;                  return false;
189          }          }
190      }      }
191      int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];      int numLines=hasline[0]+hasline[1]+hasline[2];
192    
193      // fail if we have more than one leaf group      // fail if we have more than one leaf group
194      // = there are at least two branches we can't interpolate between      // = there are at least two branches we can't interpolate between
195      if (totlines>1)      if (numLines > 1)
196          return false;          return false;
197      else if (totlines==1) {      else if (numLines==1) {
198          // we have points          // we have points
199          if (hasline[0]==1)          if (hasline[0]==1)
200              resultcode=Points;              resultcode=Points;
201          else if (hasline[1]==1) {          else if (hasline[1]==1) {
202              if (hasclass[5]==1)              if (hasclass[4])
203                  resultcode=ReducedElements;                  resultcode=ReducedElements;
204              else              else
205                  resultcode=Elements;                  resultcode=Elements;
206          } else if (hasline[2]==1) {          } else { // hasline[2]==1
207              if (hasclass[7]==1)              if (hasclass[6])
208                  resultcode=ReducedFaceElements;                  resultcode=ReducedFaceElements;
209              else              else
210                  resultcode=FaceElements;                  resultcode=FaceElements;
211          } else          }
212              throw RipleyException("Internal Ripley Error!");      } else { // numLines==0
213      } else { // totlines==0          if (hasclass[1])
214          if (hasclass[2]==1)              // something from class 1
             // something from class 2  
215              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
216          else // something from class 1          else // something from class 0
217              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
218      }      }
219      return true;      return true;
# Line 237  bool RipleyDomain::probeInterpolationOnD Line 233  bool RipleyDomain::probeInterpolationOnD
233              fsType_target != Points) {              fsType_target != Points) {
234          stringstream msg;          stringstream msg;
235          msg << "probeInterpolationOnDomain(): Invalid functionspace type "          msg << "probeInterpolationOnDomain(): Invalid functionspace type "
236              << fsType_target;              << fsType_target << " for " << getDescription();
237          throw RipleyException(msg.str());          throw RipleyException(msg.str());
238      }      }
239    
# Line 265  bool RipleyDomain::probeInterpolationOnD Line 261  bool RipleyDomain::probeInterpolationOnD
261          default: {          default: {
262              stringstream msg;              stringstream msg;
263              msg << "probeInterpolationOnDomain(): Invalid functionspace type "              msg << "probeInterpolationOnDomain(): Invalid functionspace type "
264                  << fsType_source;                  << fsType_source << " for " << getDescription();
265              throw RipleyException(msg.str());              throw RipleyException(msg.str());
266          }          }
267      }      }
# Line 377  bool RipleyDomain::isCellOriented(int fs Line 373  bool RipleyDomain::isCellOriented(int fs
373              break;              break;
374      }      }
375      stringstream msg;      stringstream msg;
376      msg << "isCellOriented(): Illegal function space type " << fsType << " on " << getDescription();      msg << "isCellOriented(): Illegal function space type " << fsType
377            << " on " << getDescription();
378      throw RipleyException(msg.str());      throw RipleyException(msg.str());
379  }  }
380    
# Line 398  bool RipleyDomain::canTag(int fsType) co Line 395  bool RipleyDomain::canTag(int fsType) co
395              break;              break;
396      }      }
397      stringstream msg;      stringstream msg;
398      msg << "canTag(): Illegal function space type " << fsType << " on " << getDescription();      msg << "canTag(): Illegal function space type " << fsType << " on "
399            << getDescription();
400      throw RipleyException(msg.str());      throw RipleyException(msg.str());
401  }  }
402    
# Line 424  void RipleyDomain::setTags(const int fsT Line 422  void RipleyDomain::setTags(const int fsT
422              break;              break;
423          default: {          default: {
424              stringstream msg;              stringstream msg;
425              msg << "setTags(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "setTags(): not implemented for "
426                    << functionSpaceTypeAsString(fsType);
427              throw RipleyException(msg.str());              throw RipleyException(msg.str());
428          }          }
429      }      }
# Line 461  int RipleyDomain::getTagFromSampleNo(int Line 460  int RipleyDomain::getTagFromSampleNo(int
460              break;              break;
461          default: {          default: {
462              stringstream msg;              stringstream msg;
463              msg << "getTagFromSampleNo(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getTagFromSampleNo(): not implemented for "
464                    << functionSpaceTypeAsString(fsType);
465              throw RipleyException(msg.str());              throw RipleyException(msg.str());
466          }          }
467      }      }
# Line 481  int RipleyDomain::getNumberOfTagsInUse(i Line 481  int RipleyDomain::getNumberOfTagsInUse(i
481              return m_faceTagsInUse.size();              return m_faceTagsInUse.size();
482          default: {          default: {
483              stringstream msg;              stringstream msg;
484              msg << "getNumberOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getNumberOfTagsInUse(): not implemented for "
485                    << functionSpaceTypeAsString(fsType);
486              throw RipleyException(msg.str());              throw RipleyException(msg.str());
487          }          }
488      }      }
# Line 500  const int* RipleyDomain::borrowListOfTag Line 501  const int* RipleyDomain::borrowListOfTag
501              return &m_faceTagsInUse[0];              return &m_faceTagsInUse[0];
502          default: {          default: {
503              stringstream msg;              stringstream msg;
504              msg << "borrowListOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "borrowListOfTagsInUse(): not implemented for "
505                    << functionSpaceTypeAsString(fsType);
506              throw RipleyException(msg.str());              throw RipleyException(msg.str());
507          }          }
508      }      }
# Line 548  escript::ASM_ptr RipleyDomain::newSystem Line 550  escript::ASM_ptr RipleyDomain::newSystem
550      // is the domain right?      // is the domain right?
551      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
552      if (row_domain!=*this)      if (row_domain!=*this)
553          throw RipleyException("Domain of row function space does not match the domain of matrix generator");          throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
554      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
555      if (col_domain!=*this)      if (col_domain!=*this)
556          throw RipleyException("Domain of column function space does not match the domain of matrix generator");          throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
557      // is the function space type right?      // is the function space type right?
558      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
559          reduceRowOrder=true;          reduceRowOrder=true;
560      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
561          throw RipleyException("Illegal function space type for system matrix rows");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
562      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
563          reduceColOrder=true;          reduceColOrder=true;
564      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
565          throw RipleyException("Illegal function space type for system matrix columns");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
566    
567      // generate matrix      // generate matrix
568      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
# Line 582  void RipleyDomain::setNewX(const escript Line 584  void RipleyDomain::setNewX(const escript
584  void RipleyDomain::copyNodalData(escript::Data& out, escript::Data& in) const  void RipleyDomain::copyNodalData(escript::Data& out, escript::Data& in) const
585  {  {
586      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
587        out.requireWrite();
588  #pragma omp parallel for  #pragma omp parallel for
589      for (index_t i=0; i<in.getNumSamples(); i++) {      for (index_t i=0; i<in.getNumSamples(); i++) {
590          const double* src = in.getSampleDataRO(i);          const double* src = in.getSampleDataRO(i);
# Line 639  void RipleyDomain::updateTagsInUse(int f Line 642  void RipleyDomain::updateTagsInUse(int f
642          local_minFoundValue = minFoundValue;          local_minFoundValue = minFoundValue;
643          MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);          MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
644  #endif  #endif
         // if we found a new value add it to the tagsInUse vector  
645    
646            // if we found a new value add it to the tagsInUse vector
647          if (minFoundValue < INDEX_T_MAX) {          if (minFoundValue < INDEX_T_MAX) {
648              tagsInUse->push_back(minFoundValue);              tagsInUse->push_back(minFoundValue);
649              lastFoundValue = minFoundValue;              lastFoundValue = minFoundValue;

Legend:
Removed from v.3733  
changed lines
  Added in v.3740

  ViewVC Help
Powered by ViewVC 1.1.26