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
00032
00033
00034
00035
00036
00037
00038 #include "SundanceCellCurvePredicate.hpp"
00039 #include "SundanceStdMathOps.hpp"
00040
00041 using namespace Sundance;
00042
00043 #define SIGN(X) ((X>0.0)?1:-1)
00044
00045 bool CellCurvePredicate::lessThan(const CellPredicateBase* other) const
00046 {
00047 const CellCurvePredicate* S = dynamic_cast<const CellCurvePredicate*>(other);
00048
00049 TEUCHOS_TEST_FOR_EXCEPTION( S== 0,
00050 std::logic_error,
00051 "argument " << other->toXML()
00052 << " to CellCurvePredicate::lessThan() should be "
00053 "a CellCurvePredicate pointer.");
00054
00055 return OrderedPair<ParametrizedCurve, int>(curve_, (int)filterMode_)
00056 < OrderedPair<ParametrizedCurve, int>(S->curve_, (int)S->filterMode_);
00057 }
00058
00059 void CellCurvePredicate::testBatch(const Array<int>& cellLID,
00060 Array<int>& results) const
00061 {
00062 results.resize(cellLID.size());
00063 double eps = 1e-14;
00064
00065 if (cellDim()==0)
00066 {
00067 switch (filterMode_){
00068 case Outside_Curve:
00069 for (int i=0; i<cellLID.size(); i++)
00070 if ( curve_.curveEquation( mesh().nodePosition(cellLID[i])) >= 0.0 )
00071 results[i] = true;
00072 else
00073 results[i] = false;
00074 break;
00075 case Inside_Curve:
00076 for (int i=0; i<cellLID.size(); i++)
00077 if ( curve_.curveEquation( mesh().nodePosition(cellLID[i])) <= 0.0 )
00078 results[i] = true;
00079 else
00080 results[i] = false;
00081 break;
00082 case On_Curve:
00083 for (int i=0; i<cellLID.size(); i++)
00084 if ( fabs(curve_.curveEquation( mesh().nodePosition(cellLID[i]))) < eps )
00085 results[i] = true;
00086 else
00087 results[i] = false;
00088 break;
00089 }
00090 }
00091 else
00092 {
00093
00094 switch (mesh().spatialDim()) {
00095 case 2:{
00096
00097
00098
00099 Array<int> facetLIDs;
00100 Array<int> facetSigns;
00101 Array<double> equSum( cellLID.size() , 0.0 );
00102 int nf = mesh().numFacets(cellDim(), cellLID[0], 0);
00103 int spaceDim = mesh().spatialDim() ;
00104 mesh().getFacetLIDs(cellDim(), cellLID, 0, facetLIDs, facetSigns);
00105 for (int c=0; c<cellLID.size(); c++)
00106 {
00107 results[c] = true;
00108 if (filterMode_ == On_Curve) results[c] = false;
00109 int curve_sign = 0;
00110 for (int f=0; f<nf; f++)
00111 {
00112 int fLID = facetLIDs[c*nf + f];
00113 double curveEqu = curve_.curveEquation( mesh().nodePosition(fLID) );
00114
00115 switch (filterMode_){
00116 case Outside_Curve:{
00117 if ( curveEqu <= 0.0 ) {
00118 results[c] = false;
00119 }
00120 equSum[c] = equSum[c] + curveEqu;
00121 break;}
00122 case Inside_Curve: {
00123 if ( curveEqu >= 0.0 ) {
00124 results[c] = false;
00125 }
00126 equSum[c] = equSum[c] + curveEqu;
00127 break; }
00128 case On_Curve:
00129 if (f == 0){
00130 curve_sign = SIGN(curveEqu);
00131 equSum[c] = equSum[c] + curveEqu;
00132 } else {
00133 if ( curve_sign != SIGN(curveEqu) ){
00134 results[c] = true;
00135 }
00136 equSum[c] = equSum[c] + curveEqu;
00137 }
00138 break;
00139 }
00140 }
00141 }
00142
00143 if ( (spaceDim == 2) && (cellDim() == 2) )
00144 {
00145 int nrEdge = mesh().numFacets(cellDim(), cellLID[0], 1 ) , nrPoints = 0 ;
00146 Array<int> edgeLIDs;
00147 Array<int> edgeLID(1);
00148 Array<int> cLID(1);
00149 Array<int> edgeSigns;
00150 Array<int> pointsOfEdgeLIDs;
00151 Array<bool> hasIntersectionPoints( cellLID.size() , false );
00152 Array<Point> intPoints;
00153
00154 for (int c=0; c<cellLID.size(); c++)
00155 {
00156 Point p0 , start, end;
00157 int nrIntPoint = 0;
00158 cLID[0] = cellLID[c];
00159 mesh().getFacetLIDs(cellDim(), cLID, 1 , edgeLIDs, edgeSigns);
00160
00161 for (int edgeI = 0 ; edgeI < nrEdge ; edgeI++ ){
00162 edgeLID[0] = edgeLIDs[edgeI];
00163 mesh().getFacetLIDs(1, edgeLID, 0 , pointsOfEdgeLIDs, edgeSigns);
00164 start = mesh().nodePosition(pointsOfEdgeLIDs[0]);
00165 end = mesh().nodePosition(pointsOfEdgeLIDs[1]);
00166
00167 curve_.returnIntersectPoints( start , end , nrPoints , intPoints);
00168 for (int p = 0 ;p < nrPoints ; p++){
00169
00170 if (nrIntPoint == 0){
00171 p0 = intPoints[p]; nrIntPoint++;
00172 } else {
00173 double dist = ::sqrt( (p0 - intPoints[p])*(p0 - intPoints[p]) );
00174 if (dist > eps){
00175
00176 hasIntersectionPoints[c] = true;
00177 continue;
00178 }
00179 }
00180 }
00181
00182 if (hasIntersectionPoints[c]) continue;
00183 }
00184 }
00185
00186
00187 int count = 0;
00188 switch (filterMode_){
00189 case Outside_Curve:{
00190 for (int c=0; c<cellLID.size(); c++){
00191
00192 if (hasIntersectionPoints[c]) results[c] = false;
00193
00194 else results[c] = (equSum[c] >= 0);
00195
00196
00197 } }
00198 break;
00199 case Inside_Curve:{
00200 for (int c=0; c<cellLID.size(); c++){
00201
00202 if (hasIntersectionPoints[c]) results[c] = false;
00203
00204 else results[c] = (equSum[c] <= 0);
00205
00206
00207 } }
00208 break;
00209 case On_Curve:{
00210 for (int c=0; c<cellLID.size(); c++){
00211
00212 if (hasIntersectionPoints[c]) results[c] = true;
00213 else results[c] = false;
00214 if (results[c]) {count++;}
00215
00216
00217
00218 } }
00219 break;
00220 }
00221 }
00222 } break;
00223
00224
00225 default:{
00226 Array<int> facetLIDs;
00227 Array<int> facetSigns;
00228 int nf = mesh().numFacets(cellDim(), cellLID[0], 0);
00229 mesh().getFacetLIDs(cellDim(), cellLID, 0, facetLIDs, facetSigns);
00230 for (int c=0; c<cellLID.size(); c++)
00231 {
00232 results[c] = true;
00233 if (filterMode_ == On_Curve) results[c] = false;
00234 int curve_sign = 0;
00235 for (int f=0; f<nf; f++)
00236 {
00237 int fLID = facetLIDs[c*nf + f];
00238 switch (filterMode_){
00239 case Outside_Curve:
00240 if ( curve_.curveEquation( mesh().nodePosition(fLID) ) <= 0.0 ) {
00241 results[c] = false;
00242 continue;
00243 }
00244 break;
00245 case Inside_Curve:
00246 if ( curve_.curveEquation( mesh().nodePosition(fLID) ) >= 0.0 ) {
00247 results[c] = false;
00248 continue;
00249 }
00250 break;
00251 case On_Curve:
00252 if (f == 0){
00253 curve_sign = SIGN(curve_.curveEquation( mesh().nodePosition(fLID)));
00254 } else {
00255 if ( curve_sign != SIGN(curve_.curveEquation( mesh().nodePosition(fLID))) ){
00256 results[c] = true;
00257 continue;
00258 }
00259 }
00260 break;
00261 }
00262 }
00263
00264 }
00265 } break;
00266 }
00267
00268
00269 }
00270 }
00271
00272 XMLObject CellCurvePredicate::toXML() const
00273 {
00274 XMLObject rtn("SundanceCellCurvePredicate");
00275 return rtn;
00276 }