/[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 3711 by caltinay, Tue Dec 6 00:24:43 2011 UTC revision 3733 by caltinay, Fri Dec 9 04:02:56 2011 UTC
# Line 89  pair<int,int> RipleyDomain::getDataShape Line 89  pair<int,int> RipleyDomain::getDataShape
89      const int ptsPerSample = (m_numDim==2 ? 4 : 8);      const int ptsPerSample = (m_numDim==2 ? 4 : 8);
90      switch (fsType) {      switch (fsType) {
91          case Nodes:          case Nodes:
92            case ReducedNodes:
93          case DegreesOfFreedom:          case DegreesOfFreedom:
94            case ReducedDegreesOfFreedom:
95              return pair<int,int>(1, getNumNodes());              return pair<int,int>(1, getNumNodes());
96          case Elements:          case Elements:
97              return pair<int,int>(ptsPerSample, getNumElements());              return pair<int,int>(ptsPerSample, getNumElements());
# Line 99  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
         case ReducedNodes:  
         case ReducedDegreesOfFreedom:  
             return pair<int,int>(1, getNumReducedNodes());  
             */  
106          default:          default:
107              break;              break;
108      }      }
# Line 116  pair<int,int> RipleyDomain::getDataShape Line 113  pair<int,int> RipleyDomain::getDataShape
113      throw RipleyException(msg.str());      throw RipleyException(msg.str());
114  }  }
115    
 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const  
 {  
     throw RipleyException("getTagFromSampleNo() not implemented");  
 }  
   
116  string RipleyDomain::showTagNames() const  string RipleyDomain::showTagNames() const
117  {  {
118      stringstream ret;      stringstream ret;
# Line 231  bool RipleyDomain::commonFunctionSpace(c Line 223  bool RipleyDomain::commonFunctionSpace(c
223      return true;      return true;
224  }  }
225    
226    bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
227                                                  int fsType_target) const
228    {
229        if (fsType_target != Nodes &&
230                fsType_target != ReducedNodes &&
231                fsType_target != ReducedDegreesOfFreedom &&
232                fsType_target != DegreesOfFreedom &&
233                fsType_target != Elements &&
234                fsType_target != ReducedElements &&
235                fsType_target != FaceElements &&
236                fsType_target != ReducedFaceElements &&
237                fsType_target != Points) {
238            stringstream msg;
239            msg << "probeInterpolationOnDomain(): Invalid functionspace type "
240                << fsType_target;
241            throw RipleyException(msg.str());
242        }
243    
244        switch (fsType_source) {
245            case Nodes:
246            case DegreesOfFreedom:
247                return true;
248            case ReducedNodes:
249            case ReducedDegreesOfFreedom:
250                return (fsType_target != Nodes &&
251                        fsType_target != DegreesOfFreedom);
252            case Elements:
253                return (fsType_target==Elements ||
254                        fsType_target==ReducedElements);
255            case ReducedElements:
256                return (fsType_target==ReducedElements);
257            case FaceElements:
258                return (fsType_target==FaceElements ||
259                        fsType_target==ReducedFaceElements);
260            case ReducedFaceElements:
261                return (fsType_target==ReducedFaceElements);
262            case Points:
263                return (fsType_target==Points);
264    
265            default: {
266                stringstream msg;
267                msg << "probeInterpolationOnDomain(): Invalid functionspace type "
268                    << fsType_source;
269                throw RipleyException(msg.str());
270            }
271        }
272    }
273    
274  void RipleyDomain::interpolateOnDomain(escript::Data& target,  void RipleyDomain::interpolateOnDomain(escript::Data& target,
275                                         const escript::Data& in) const                                         const escript::Data& in) const
276  {  {
# Line 249  void RipleyDomain::interpolateOnDomain(e Line 289  void RipleyDomain::interpolateOnDomain(e
289    
290      switch (in.getFunctionSpace().getTypeCode()) {      switch (in.getFunctionSpace().getTypeCode()) {
291          case Nodes:          case Nodes:
292            case ReducedNodes:
293          case DegreesOfFreedom:          case DegreesOfFreedom:
294            case ReducedDegreesOfFreedom:
295              switch (target.getFunctionSpace().getTypeCode()) {              switch (target.getFunctionSpace().getTypeCode()) {
296                    case Nodes:
297                    case ReducedNodes:
298                  case DegreesOfFreedom:                  case DegreesOfFreedom:
299                      // FIXME!                  case ReducedDegreesOfFreedom:
300                      target=in;                      // FIXME
301                        copyNodalData(target, *const_cast<escript::Data*>(&in));
302                      break;                      break;
303    
304                  case Elements:                  case Elements:
# Line 332  bool RipleyDomain::isCellOriented(int fs Line 377  bool RipleyDomain::isCellOriented(int fs
377              break;              break;
378      }      }
379      stringstream msg;      stringstream msg;
380      msg << "Illegal function space type " << fsType << " on " << getDescription();      msg << "isCellOriented(): Illegal function space type " << fsType << " on " << getDescription();
381      throw RipleyException(msg.str());      throw RipleyException(msg.str());
382  }  }
383    
384    bool RipleyDomain::canTag(int fsType) const
385    {
386        switch(fsType) {
387            case Nodes:
388            case Elements:
389            case ReducedElements:
390            case FaceElements:
391            case ReducedFaceElements:
392                return true;
393            case DegreesOfFreedom:
394            case ReducedDegreesOfFreedom:
395            case Points:
396                return false;
397            default:
398                break;
399        }
400        stringstream msg;
401        msg << "canTag(): Illegal function space type " << fsType << " on " << getDescription();
402        throw RipleyException(msg.str());
403    }
404    
405    void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
406    {
407        IndexVector* target=NULL;
408        dim_t num=0;
409    
410        switch(fsType) {
411            case Nodes:
412                num=getNumNodes();
413                target=&m_nodeTags;
414                break;
415            case Elements:
416            case ReducedElements:
417                num=getNumElements();
418                target=&m_elementTags;
419                break;
420            case FaceElements:
421            case ReducedFaceElements:
422                num=getNumFaceElements();
423                target=&m_faceTags;
424                break;
425            default: {
426                stringstream msg;
427                msg << "setTags(): not implemented for " << functionSpaceTypeAsString(fsType);
428                throw RipleyException(msg.str());
429            }
430        }
431        if (target->size() != num) {
432            target->assign(num, -1);
433        }
434    
435        escript::Data& mask=*const_cast<escript::Data*>(&cMask);
436    #pragma omp parallel for
437        for (index_t i=0; i<num; i++) {
438            if (mask.getSampleDataRO(i)[0] > 0) {
439                (*target)[i]=newTag;
440            }
441        }
442        updateTagsInUse(fsType);
443    }
444    
445    int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
446    {
447        switch(fsType) {
448            case Nodes:
449                if (m_nodeTags.size() > sampleNo)
450                    return m_nodeTags[sampleNo];
451                break;
452            case Elements:
453            case ReducedElements:
454                if (m_elementTags.size() > sampleNo)
455                    return m_elementTags[sampleNo];
456                break;
457            case FaceElements:
458            case ReducedFaceElements:
459                if (m_faceTags.size() > sampleNo)
460                    return m_faceTags[sampleNo];
461                break;
462            default: {
463                stringstream msg;
464                msg << "getTagFromSampleNo(): not implemented for " << functionSpaceTypeAsString(fsType);
465                throw RipleyException(msg.str());
466            }
467        }
468        return -1;
469    }
470    
471    int RipleyDomain::getNumberOfTagsInUse(int fsType) const
472    {
473        switch(fsType) {
474            case Nodes:
475                return m_nodeTagsInUse.size();
476            case Elements:
477            case ReducedElements:
478                return m_elementTagsInUse.size();
479            case FaceElements:
480            case ReducedFaceElements:
481                return m_faceTagsInUse.size();
482            default: {
483                stringstream msg;
484                msg << "getNumberOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);
485                throw RipleyException(msg.str());
486            }
487        }
488    }
489    
490    const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
491    {
492        switch(fsType) {
493            case Nodes:
494                return &m_nodeTagsInUse[0];
495            case Elements:
496            case ReducedElements:
497                return &m_elementTagsInUse[0];
498            case FaceElements:
499            case ReducedFaceElements:
500                return &m_faceTagsInUse[0];
501            default: {
502                stringstream msg;
503                msg << "borrowListOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);
504                throw RipleyException(msg.str());
505            }
506        }
507    }
508    
509  void RipleyDomain::Print_Mesh_Info(const bool full) const  void RipleyDomain::Print_Mesh_Info(const bool full) const
510  {  {
511      cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "      cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
# Line 408  void RipleyDomain::setNewX(const escript Line 578  void RipleyDomain::setNewX(const escript
578      throw RipleyException("setNewX(): Operation not supported");      throw RipleyException("setNewX(): Operation not supported");
579  }  }
580    
581    //protected
582    void RipleyDomain::copyNodalData(escript::Data& out, escript::Data& in) const
583    {
584        const dim_t numComp = in.getDataPointSize();
585    #pragma omp parallel for
586        for (index_t i=0; i<in.getNumSamples(); i++) {
587            const double* src = in.getSampleDataRO(i);
588            copy(src, src+numComp, out.getSampleDataRW(i));
589        }
590    }
591    
592    //protected
593    void RipleyDomain::updateTagsInUse(int fsType) const
594    {
595        IndexVector* tagsInUse=NULL;
596        const IndexVector* tags=NULL;
597        switch(fsType) {
598            case Nodes:
599                tags=&m_nodeTags;
600                tagsInUse=&m_nodeTagsInUse;
601                break;
602            case Elements:
603            case ReducedElements:
604                tags=&m_elementTags;
605                tagsInUse=&m_elementTagsInUse;
606                break;
607            case FaceElements:
608            case ReducedFaceElements:
609                tags=&m_faceTags;
610                tagsInUse=&m_faceTagsInUse;
611                break;
612            default:
613                return;
614        }
615    
616        // gather global unique tag values from tags into tagsInUse
617        tagsInUse->clear();
618        index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
619    
620        while (true) {
621            // find smallest value bigger than lastFoundValue
622            minFoundValue = INDEX_T_MAX;
623    #pragma omp parallel private(local_minFoundValue)
624            {
625                local_minFoundValue = minFoundValue;
626    #pragma omp for schedule(static) nowait
627                for (size_t i = 0; i < tags->size(); i++) {
628                    const index_t v = (*tags)[i];
629                    if ((v > lastFoundValue) && (v < local_minFoundValue))
630                        local_minFoundValue = v;
631                }
632    #pragma omp critical
633                {
634                    if (local_minFoundValue < minFoundValue)
635                        minFoundValue = local_minFoundValue;
636                }
637            }
638    #ifdef ESYS_MPI
639            local_minFoundValue = minFoundValue;
640            MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
641    #endif
642            // if we found a new value add it to the tagsInUse vector
643    
644            if (minFoundValue < INDEX_T_MAX) {
645                tagsInUse->push_back(minFoundValue);
646                lastFoundValue = minFoundValue;
647            } else
648                break;
649        }
650    }
651    
652  //  //
653  // the methods that follow have to be implemented by the subclasses  // the methods that follow have to be implemented by the subclasses
654  //  //
# Line 437  const int* RipleyDomain::borrowSampleRef Line 678  const int* RipleyDomain::borrowSampleRef
678      throw RipleyException("borrowSampleReferenceIDs() not implemented");      throw RipleyException("borrowSampleReferenceIDs() not implemented");
679  }  }
680    
 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,  
                                              int fsType_target) const  
 {  
     throw RipleyException("probeInterpolationOnDomain() not implemented");  
 }  
   
681  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
682  {  {
683      throw RipleyException("interpolateACross() not implemented");      throw RipleyException("interpolateACross() not implemented");
# Line 479  bool RipleyDomain::ownSample(int fsType, Line 714  bool RipleyDomain::ownSample(int fsType,
714      throw RipleyException("ownSample() not implemented");      throw RipleyException("ownSample() not implemented");
715  }  }
716    
 void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& mask) const  
 {  
     throw RipleyException("setTags() not implemented");  
 }  
   
 int RipleyDomain::getNumberOfTagsInUse(int fsType) const  
 {  
     throw RipleyException("getNumberOfTagsInUse() not implemented");  
 }  
   
 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const  
 {  
     throw RipleyException("borrowListOfTagsInUse() not implemented");  
 }  
   
 bool RipleyDomain::canTag(int fsType) const  
 {  
     throw RipleyException("canTag() not implemented");  
 }  
   
   
717  void RipleyDomain::addPDEToSystem(  void RipleyDomain::addPDEToSystem(
718          escript::AbstractSystemMatrix& mat, escript::Data& rhs,          escript::AbstractSystemMatrix& mat, escript::Data& rhs,
719          const escript::Data& A, const escript::Data& B, const escript::Data& C,          const escript::Data& A, const escript::Data& B, const escript::Data& C,

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

  ViewVC Help
Powered by ViewVC 1.1.26