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

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

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

revision 3722 by caltinay, Wed Dec 7 05:53:22 2011 UTC revision 3740 by caltinay, Tue Dec 13 01:37:21 2011 UTC
# Line 101  pair<int,int> RipleyDomain::getDataShape Line 101  pair<int,int> RipleyDomain::getDataShape
101              return pair<int,int>(1, getNumElements());              return pair<int,int>(1, getNumElements());
102          case ReducedFaceElements:          case ReducedFaceElements:
103              return pair<int,int>(1, getNumFaceElements());              return pair<int,int>(1, getNumFaceElements());
             /*  
104          case Points:          case Points:
105              return pair<int,int>(1, getNumPoints());              return pair<int,int>(1, 0); //FIXME: dirac
             */  
106          default:          default:
107              break;              break;
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 131  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;
220  }  }
221    
222    bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
223                                                  int fsType_target) const
224    {
225        if (fsType_target != Nodes &&
226                fsType_target != ReducedNodes &&
227                fsType_target != ReducedDegreesOfFreedom &&
228                fsType_target != DegreesOfFreedom &&
229                fsType_target != Elements &&
230                fsType_target != ReducedElements &&
231                fsType_target != FaceElements &&
232                fsType_target != ReducedFaceElements &&
233                fsType_target != Points) {
234            stringstream msg;
235            msg << "probeInterpolationOnDomain(): Invalid functionspace type "
236                << fsType_target << " for " << getDescription();
237            throw RipleyException(msg.str());
238        }
239    
240        switch (fsType_source) {
241            case Nodes:
242            case DegreesOfFreedom:
243                return true;
244            case ReducedNodes:
245            case ReducedDegreesOfFreedom:
246                return (fsType_target != Nodes &&
247                        fsType_target != DegreesOfFreedom);
248            case Elements:
249                return (fsType_target==Elements ||
250                        fsType_target==ReducedElements);
251            case ReducedElements:
252                return (fsType_target==ReducedElements);
253            case FaceElements:
254                return (fsType_target==FaceElements ||
255                        fsType_target==ReducedFaceElements);
256            case ReducedFaceElements:
257                return (fsType_target==ReducedFaceElements);
258            case Points:
259                return (fsType_target==Points);
260    
261            default: {
262                stringstream msg;
263                msg << "probeInterpolationOnDomain(): Invalid functionspace type "
264                    << fsType_source << " for " << getDescription();
265                throw RipleyException(msg.str());
266            }
267        }
268    }
269    
270  void RipleyDomain::interpolateOnDomain(escript::Data& target,  void RipleyDomain::interpolateOnDomain(escript::Data& target,
271                                         const escript::Data& in) const                                         const escript::Data& in) const
272  {  {
# Line 331  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 352  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 378  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 415  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 435  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 454  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 502  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 536  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 593  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;
# Line 632  const int* RipleyDomain::borrowSampleRef Line 681  const int* RipleyDomain::borrowSampleRef
681      throw RipleyException("borrowSampleReferenceIDs() not implemented");      throw RipleyException("borrowSampleReferenceIDs() not implemented");
682  }  }
683    
 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,  
                                              int fsType_target) const  
 {  
     throw RipleyException("probeInterpolationOnDomain() not implemented");  
 }  
   
684  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
685  {  {
686      throw RipleyException("interpolateACross() not implemented");      throw RipleyException("interpolateACross() not implemented");

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

  ViewVC Help
Powered by ViewVC 1.1.26