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 "SundanceCellSet.hpp"
00032 #include "SundanceExplicitCellSet.hpp"
00033 #include "SundanceImplicitCellSet.hpp"
00034 #include "SundanceOut.hpp"
00035 #include "PlayaTabs.hpp"
00036 #include "PlayaExceptions.hpp"
00037 #include <algorithm>
00038 #include <iterator>
00039
00040 using namespace Sundance;
00041 using namespace Teuchos;
00042 using Playa::Handle;
00043 using Playa::Handleable;
00044
00045
00046 CellSet::CellSet(const Mesh& mesh, int cellDim,
00047 const CellType& cellType,
00048 const Set<int>& cellLIDs)
00049 : Handle<CellSetBase>(rcp(new ExplicitCellSet(mesh, cellDim, cellType, cellLIDs)))
00050 {}
00051
00052 CellSet CellSet::setUnion(const CellSet& other) const
00053 {
00054 if (isNull()) return other;
00055 if (other.isNull()) return *this;
00056
00057 ExplicitCellSet* rtn = new ExplicitCellSet(mesh(), dimension(), cellType());
00058
00059 checkCompatibility("union", other);
00060
00061 Set<int>& cells = rtn->cells();
00062
00063 std::set_union(this->begin(), this->end(), other.begin(), other.end(),
00064 std::insert_iterator<Set<int> >(cells, cells.begin()));
00065
00066 return rtn;
00067 }
00068
00069 CellSet CellSet::setIntersection(const CellSet& other) const
00070 {
00071 if (isNull()) return *this;
00072 if (other.isNull()) return other;
00073
00074 ExplicitCellSet* rtn = new ExplicitCellSet(mesh(), dimension(), cellType());
00075
00076 checkCompatibility("intersection", other);
00077
00078 Set<int>& cells = rtn->cells();
00079
00080 std::set_intersection(this->begin(), this->end(), other.begin(), other.end(),
00081 std::insert_iterator<Set<int> >(cells, cells.begin()));
00082
00083 return rtn;
00084 }
00085
00086 CellSet CellSet::setDifference(const CellSet& other) const
00087 {
00088 if (isNull()) return *this;
00089 if (other.isNull()) return *this;
00090
00091 ExplicitCellSet* rtn = new ExplicitCellSet(mesh(), dimension(), cellType());
00092
00093 checkCompatibility("difference", other);
00094
00095 Set<int>& cells = rtn->cells();
00096
00097 std::set_difference(this->begin(), this->end(), other.begin(), other.end(),
00098 std::insert_iterator<Set<int> >(cells, cells.begin()));
00099
00100 return rtn;
00101 }
00102
00103
00104 void CellSet::checkCompatibility(const std::string& op, const CellSet& other) const
00105 {
00106
00107 TEUCHOS_TEST_FOR_EXCEPTION(meshID() != other.meshID(), std::runtime_error,
00108 "CellSet::checkCompatibility(): "
00109 "incompatible mesh ID numbers in " << op
00110 << ". LHS=" << meshID() << " RHS=" << other.meshID());
00111
00112 TEUCHOS_TEST_FOR_EXCEPTION(dimension() != other.dimension(), std::runtime_error,
00113 "CellSet::checkCompatibility() incompatible dimensions in " << op
00114 << "LHS has "
00115 "dimension=" << dimension() << " but RHS has dimension="
00116 << other.dimension());
00117
00118 TEUCHOS_TEST_FOR_EXCEPTION(cellType() != other.cellType(), std::runtime_error,
00119 "CellSet::checkCompatibility() incompatible cell types. "
00120 " in " << op << " LHS has "
00121 "cellType=" << cellType() << " but RHS has cellType="
00122 << other.cellType());
00123
00124
00125 SUNDANCE_OUT(this->verb() > 2,
00126 "Set operation: " << op << std::endl
00127 << "LHS cells: " << *this << std::endl
00128 << "RHS cells: " << other);
00129
00130 }
00131
00132
00133 bool CellSet::areFacetsOf(const CellSet& other) const
00134 {
00135 Array<int> cofacetLIDs;
00136 int myDim = dimension();
00137 int cofacetDim = other.dimension();
00138 CellType cofacetType = other.cellType();
00139 if (myDim >= cofacetDim) return false;
00140
00141 for (CellIterator i=begin(); i!=end(); i++)
00142 {
00143 int cellLID = *i;
00144 mesh().getCofacets(myDim, cellLID, cofacetDim, cofacetLIDs);
00145 Set<int> cofacetSet;
00146 for (int c=0; c<cofacetLIDs.size(); c++)
00147 {
00148 int cf = cofacetLIDs[c];
00149 cofacetSet.put(cf);
00150 }
00151 CellSet myCofacetSet(mesh(), cofacetDim, cofacetType, cofacetSet);
00152 CellSet intersection = myCofacetSet.setIntersection(other);
00153
00154
00155 if (intersection.begin()==intersection.end()) return false;
00156 }
00157 return true;
00158 }
00159
00160 bool CellSet::operator<(const CellSet& other) const
00161 {
00162 Tabs tab;
00163 bool rtn = ptr()->lessThan(other.ptr().get());
00164 return rtn;
00165 }
00166
00167
00168 int CellSet::numCells() const
00169 {
00170 int count = 0;
00171 for (CellIterator i=begin(); i!=end(); i++)
00172 {
00173 count ++;
00174 }
00175 return count;
00176 }
00177