00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 
00014 
00015 
00016 
00017 
00018 
00019 
00020 
00021 
00022 
00023 
00024 
00025 
00026 
00027 
00028 
00029 
00030 
00031 #include "SundanceCFMeshPair.hpp"
00032 #include "PlayaTabs.hpp"
00033 #include "SundanceExplicitCellSet.hpp"
00034 #include "SundanceBinaryCellFilter.hpp"
00035 #include "SundanceSubsetCellFilter.hpp"
00036 #include "SundanceLabelCellPredicate.hpp"
00037 #include "SundanceNullCellFilterStub.hpp"
00038 #include "SundanceNullCellFilter.hpp"
00039 
00040 #include "Teuchos_Time.hpp"
00041 #include "Teuchos_TimeMonitor.hpp"
00042 
00043 using namespace Sundance;
00044 using namespace Sundance;
00045 using namespace Sundance;
00046 using namespace Teuchos;
00047 
00048 
00049 static Time& csPartitionTimer() 
00050 {
00051   static RCP<Time> rtn 
00052     = TimeMonitor::getNewTimer("cell set partitioning"); 
00053   return *rtn;
00054 }
00055 
00056 
00057 CFMeshPair::CFMeshPair()
00058   : filter_(),
00059     mesh_(),
00060     cellSet_(),
00061     funcs_()
00062 {}
00063 
00064 CFMeshPair::CFMeshPair(const CellFilter& cf,
00065                        const Mesh& mesh,
00066                        const Set<int>& funcs)
00067   : filter_(cf),
00068     mesh_(mesh),
00069     cellSet_(),
00070     funcs_(funcs)
00071 {
00072   if (filter_.ptr().get() != 0) cellSet_ = filter_.getCells(mesh);
00073 }
00074 
00075 bool CFMeshPair::operator<(const CFMeshPair& other) const
00076 {
00077   if (isEmpty()) return true;
00078   if (other.isEmpty()) return false;
00079 
00080   TEUCHOS_TEST_FOR_EXCEPTION(mesh_.id() != other.mesh_.id(),
00081                      std::runtime_error,
00082                      "mismatched meshes!");
00083 
00084   return cellSet_ < other.cellSet_;
00085 }
00086 
00087 bool CFMeshPair::isEmpty() const
00088 {
00089   return filter_.ptr().get() == 0 
00090     || cellSet_.ptr().get() == 0
00091     || cellSet_.begin()==cellSet_.end();
00092 }
00093 
00094 CFMeshPair CFMeshPair::setMinus(const CFMeshPair& other) const
00095 {
00096 
00097   if (isEmpty()) return CFMeshPair();
00098   if (other.isEmpty()) return *this;
00099 
00100   TEUCHOS_TEST_FOR_EXCEPTION(mesh().id() != other.mesh().id(),
00101                      std::runtime_error,
00102                      "mismatched meshes!");
00103 
00104   CellFilter diff = filter() - other.filter();
00105   return CFMeshPair(diff, mesh(), funcs_);
00106 }
00107 
00108 CFMeshPair CFMeshPair::intersection(const CFMeshPair& other) const
00109 {
00110   if (isEmpty() || other.isEmpty()) return CFMeshPair();
00111 
00112   TEUCHOS_TEST_FOR_EXCEPTION(mesh().id() != other.mesh().id(),
00113                      std::runtime_error,
00114                      "mismatched meshes!");
00115 
00116   CellFilter inter = filter().intersection(other.filter());
00117 
00118   return CFMeshPair(inter, mesh(), funcs_.setUnion(other.funcs_));
00119 }
00120 
00121 
00122 namespace Sundance
00123 {
00124   Array<CFMeshPair> resolvePair(const CFMeshPair& a, 
00125                                 const CFMeshPair& b)
00126   {
00127     Tabs tab0;
00128 
00129     CFMeshPair inter = a.intersection(b);
00130     CFMeshPair aMinusB = a.setMinus(b);
00131     CFMeshPair bMinusA = b.setMinus(a);
00132 
00133     return tuple(bMinusA, inter, aMinusB);
00134   }
00135 
00136   
00137   Array<CFMeshPair> resolveSets(const Array<CFMeshPair>& s)
00138   {
00139 
00140     if (s.size() == 1) return s;
00141 
00142     Array<CFMeshPair> T = tuple(s[0]);
00143     
00144     for (int i=0; i<s.size(); i++)
00145       {
00146         CFMeshPair A = s[i];
00147         Array<CFMeshPair> TN;
00148         for (int j=0; j<T.size(); j++)
00149           {
00150             Array<CFMeshPair> p = resolvePair(A, T[j]);
00151             if (!p[0].isEmpty())
00152               {
00153                 TN.append(p[0]);
00154               }
00155             if (!p[1].isEmpty())
00156               {
00157                 TN.append(p[1]);
00158               }
00159             A = p[2];
00160           }
00161         if (!A.isEmpty())
00162           {
00163             TN.append(A);
00164           }
00165         T = TN;
00166       }
00167     return T;
00168   }
00169 
00170 
00171   Array<CFMeshPair>
00172   findDisjointFilters(const Array<CellFilter>& filters,
00173                       const Array<Set<int> >& funcs,
00174                       const Mesh& mesh)
00175   {
00176     TimeMonitor timer(csPartitionTimer() );
00177     Array<CFMeshPair> cf(filters.size());
00178 
00179     TEUCHOS_TEST_FOR_EXCEPT(filters.size() != funcs.size());
00180     for (int i=0; i<filters.size(); i++)
00181       {
00182         cf[i] = CFMeshPair(filters[i], mesh, funcs[i]);
00183       }
00184     return resolveSets(cf);
00185   }
00186 
00187 }