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 "SundanceGeomUtils.hpp"
00032 #include "SundanceMesh.hpp"
00033 #include "SundancePoint.hpp"
00034 #include "SundanceSet.hpp"
00035 #include <queue>
00036
00037 using namespace Sundance;
00038 using namespace Teuchos;
00039 using namespace Sundance;
00040
00041
00042 namespace Sundance
00043 {
00044 bool cellContainsPoint(const Mesh& mesh, int cellDim, int cellLID,
00045 const double* x, Array<int>& facetLID)
00046 {
00047 if (cellDim==1)
00048 {
00049 static Array<int> facetLID(2);
00050 static Array<int> tmp(2);
00051 mesh.getFacetArray(cellDim, cellLID, 0, facetLID, tmp);
00052 Point A = mesh.nodePosition(facetLID[0]);
00053 Point B = mesh.nodePosition(facetLID[1]);
00054 if (A[0] < B[0]) return (x[0] >= A[0] && x[0] <= B[0]);
00055 else return (x[0] <= A[0] && x[0] >= B[0]);
00056 }
00057 else if (cellDim==2)
00058 {
00059 static Array<int> tmp(3);
00060 mesh.getFacetArray(cellDim, cellLID, 0, facetLID, tmp);
00061 const double* A = mesh.nodePositionView(facetLID[0]);
00062 const double* B = mesh.nodePositionView(facetLID[1]);
00063 const double* C = mesh.nodePositionView(facetLID[2]);
00064
00065
00066 double sign = orient2D(A, B, C);
00067 if (sign > 0.0)
00068 {
00069 double s1 = orient2D(A, B, x);
00070 double s2 = orient2D(B, C, x);
00071 double s3 = orient2D(C, A, x);
00072 if (s1 >= 0.0 && s2 >= 0.0 && s3 >= 0.0) return true;
00073 return false;
00074 }
00075 else
00076 {
00077 double s1 = orient2D(A, C, x);
00078 double s2 = orient2D(C, B, x);
00079 double s3 = orient2D(B, A, x);
00080 if (s1 >= 0.0 && s2 >= 0.0 && s3 >= 0.0) return true;
00081 return false;
00082 }
00083 }
00084 else if (cellDim==3)
00085 {
00086 TEUCHOS_TEST_FOR_EXCEPT(true);
00087 return false;
00088 }
00089 else
00090 {
00091 TEUCHOS_TEST_FOR_EXCEPTION(cellDim<=0 || cellDim>3, std::runtime_error,
00092 "invalid point dimension " << cellDim);
00093 return false;
00094 }
00095 }
00096
00097 double volume(const Mesh& mesh, int cellDim, int cellLID)
00098 {
00099 if (cellDim==1)
00100 {
00101 static Array<int> facetLID(2);
00102 static Array<int> tmp(2);
00103 mesh.getFacetArray(cellDim, cellLID, 0, facetLID, tmp);
00104 Point A = mesh.nodePosition(facetLID[0]);
00105 Point B = mesh.nodePosition(facetLID[1]);
00106 return fabs(A[0] - B[0]);
00107 }
00108 else if (cellDim==2)
00109 {
00110 static Array<int> facetLID(3);
00111 static Array<int> tmp(3);
00112 mesh.getFacetArray(cellDim, cellLID, 0, facetLID, tmp);
00113 const double* A = mesh.nodePositionView(facetLID[0]);
00114 const double* B = mesh.nodePositionView(facetLID[1]);
00115 const double* C = mesh.nodePositionView(facetLID[2]);
00116 return fabs(0.5*orient2D(A, B, C));
00117 }
00118 TEUCHOS_TEST_FOR_EXCEPT(true);
00119 }
00120
00121
00122 double orient2D(const double* A, const double* B, const double* x)
00123 {
00124 double acx, bcx, acy, bcy;
00125
00126 acx = A[0] - x[0];
00127 bcx = B[0] - x[0];
00128 acy = A[1] - x[1];
00129 bcy = B[1] - x[1];
00130 return acx * bcy - acy * bcx;
00131 }
00132
00133 int findEnclosingCell(const Mesh& mesh, int cellDim,
00134 int initialGuessLID,
00135 const double* x)
00136 {
00137 std::queue<int> Q;
00138 Set<int> repeats;
00139 static Array<int> facets;
00140
00141 Q.push(initialGuessLID);
00142
00143
00144
00145 while (!Q.empty())
00146 {
00147 int next = Q.front();
00148 Q.pop();
00149 if (repeats.contains(next)) continue;
00150
00151 if (cellContainsPoint(mesh, cellDim, next, x, facets)) return next;
00152 repeats.put(next);
00153
00154 std::list<int> neighbors;
00155 maximalNeighbors(mesh, cellDim, next, facets, neighbors);
00156 for (std::list<int>::const_iterator
00157 i=neighbors.begin(); i!=neighbors.end(); i++)
00158 {
00159 Q.push(*i);
00160 }
00161 }
00162 return -1;
00163 }
00164
00165 Point pullback(const Mesh& mesh, int cellDim, int cellLID, const double* x)
00166 {
00167 int dim = cellDim;
00168 if (dim==2)
00169 {
00170 static Array<int> facetLID(3);
00171 static Array<int> tmp(3);
00172 mesh.getFacetArray(cellDim, cellLID, 0, facetLID, tmp);
00173 const double* A = mesh.nodePositionView(facetLID[0]);
00174 const double* B = mesh.nodePositionView(facetLID[1]);
00175 const double* C = mesh.nodePositionView(facetLID[2]);
00176
00177 double bax = B[0] - A[0];
00178 double bay = B[1] - A[1];
00179 double cax = C[0] - A[0];
00180 double cay = C[1] - A[1];
00181 double delta = bax*cay - bay*cax;
00182
00183 double xax = x[0] - A[0];
00184 double xay = x[1] - A[1];
00185
00186 return Point( (cay*xax - cax*xay)/delta,
00187 (-bay*xax + bax*xay)/delta );
00188
00189 }
00190 else
00191 {
00192 TEUCHOS_TEST_FOR_EXCEPTION(cellDim != 2, std::runtime_error,
00193 "invalid point dimension " << cellDim);
00194 return Point();
00195 }
00196 }
00197
00198 void printCell(const Mesh& mesh, int cellLID)
00199 {
00200 if (mesh.spatialDim()==2)
00201 {
00202 static Array<int> facetLID(3);
00203 static Array<int> tmp(3);
00204 mesh.getFacetArray(mesh.spatialDim(), cellLID, 0, facetLID, tmp);
00205 Point A = mesh.nodePosition(facetLID[0]);
00206 Point B = mesh.nodePosition(facetLID[1]);
00207 Point C = mesh.nodePosition(facetLID[2]);
00208 std::cout << "{" << A << ", " << B << ", " << C << "}" << std::endl;
00209 }
00210 }
00211
00212 void maximalNeighbors(const Mesh& mesh, int cellDim, int cellLID,
00213 const Array<int>& facetLID,
00214 std::list<int>& rtn)
00215 {
00216 for (int f=0; f<facetLID.size(); f++)
00217 {
00218 Array<int> cofacets;
00219 mesh.getCofacets(0, facetLID[f], cellDim, cofacets);
00220 for (int c=0; c<cofacets.size(); c++)
00221 {
00222 if (cofacets[c] != cellLID) rtn.push_back(cofacets[c]);
00223 }
00224 }
00225 }
00226
00227
00228 }