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 "SundanceMeshBase.hpp"
00032 #include "SundanceMesh.hpp"
00033 #include "PlayaExceptions.hpp"
00034 #include "Teuchos_Time.hpp"
00035 #include "Teuchos_TimeMonitor.hpp"
00036 #include "PlayaTabs.hpp"
00037
00038 using namespace Sundance;
00039 using namespace Sundance;
00040 using namespace Teuchos;
00041 using namespace Sundance;
00042
00043
00044
00045 MeshBase::MeshBase(int dim, const MPIComm& comm,
00046 const MeshEntityOrder& order)
00047 : dim_(dim),
00048 comm_(comm),
00049 order_(order),
00050 reorderer_(Mesh::defaultReorderer().createInstance(this)),
00051 validWeights_(true),
00052 specialWeights_(),
00053 curvePoints_Are_Valid_(true),
00054 nrCurvesForIntegral_(0),
00055 curvePoints_() ,
00056 curveDerivative_(),
00057 curveNormal_(),
00058 curveID_to_ArrayIndex_()
00059 { specialWeights_.resize(dim_);
00060 curvePoints_.resize(0);
00061 curveDerivative_.resize(0);
00062 curveNormal_.resize(0);
00063 }
00064
00065 const int* MeshBase::elemZeroFacetView(int maxCellLID) const
00066 {
00067 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00068 "MeshBase::elemZeroFacetView() is deprecated");
00069 return (const int*) 0;
00070 }
00071
00072
00073 Point MeshBase::centroid(int cellDim, int cellLID) const
00074 {
00075 if (cellDim==0) return nodePosition(cellLID);
00076 int dummy;
00077 Point x = nodePosition(facetLID(cellDim, cellLID, 0, 0, dummy));
00078 int nf = numFacets(cellDim, cellLID, 0);
00079 for (int f=1; f<nf; f++)
00080 x += nodePosition(facetLID(cellDim, cellLID, 0, f, dummy));
00081 return x / ((double) nf);
00082 }
00083
00084 void MeshBase::outwardNormals(
00085 const Array<int>& cellLIDs,
00086 Array<Point>& outwardNormals
00087 ) const
00088 {
00089 int D = spatialDim();
00090 outwardNormals.resize(cellLIDs.size());
00091 for (int c=0; c<cellLIDs.size(); c++)
00092 {
00093 int f=-1;
00094 TEUCHOS_TEST_FOR_EXCEPTION(numMaxCofacets(D-1, cellLIDs[c]) > 1,
00095 std::runtime_error,
00096 "cell #" << cellLIDs[c] << " is not a boundary cell");
00097 int maxLID = maxCofacetLID(D-1, cellLIDs[c], 0, f);
00098 Point cInterior = centroid(D, maxLID);
00099 Point cBdry = centroid(D-1, cellLIDs[c]);
00100 Point q = cBdry - cInterior;
00101 Point s;
00102 if (D==1)
00103 {
00104 s = Point(1.0);
00105 }
00106 else if (D==2)
00107 {
00108 Point A = nodePosition(facetLID(D-1, cellLIDs[c], 0, 0, f));
00109 Point B = nodePosition(facetLID(D-1, cellLIDs[c], 0, 1, f));
00110 Point t = B - A;
00111 s = Point(-t[1], t[0]);
00112 }
00113 else
00114 {
00115 Point A = nodePosition(facetLID(D-1, cellLIDs[c], 0, 0, f));
00116 Point B = nodePosition(facetLID(D-1, cellLIDs[c], 0, 1, f));
00117 Point C = nodePosition(facetLID(D-1, cellLIDs[c], 0, 2, f));
00118 s = cross(B-A, C-A);
00119 }
00120 if (q * s > 0.0)
00121 {
00122 outwardNormals[c] = s/::sqrt(s*s);
00123 }
00124 else
00125 {
00126 outwardNormals[c] = -s/::sqrt(s*s);
00127 }
00128 }
00129 }
00130
00131
00132 void MeshBase::tangentsToEdges(
00133 const Array<int>& cellLIDs,
00134 Array<Point>& tangentVectors
00135 ) const
00136 {
00137 TEUCHOS_TEST_FOR_EXCEPT(spatialDim() <= 1);
00138
00139 tangentVectors.resize(cellLIDs.size());
00140
00141 for (int c=0; c<cellLIDs.size(); c++)
00142 {
00143 int fOrient=1;
00144 Point A = nodePosition(facetLID(1, cellLIDs[c], 0, 0, fOrient));
00145 Point B = nodePosition(facetLID(1, cellLIDs[c], 0, 1, fOrient));
00146 Point t = B - A;
00147 tangentVectors[c] = t/(sqrt(t*t));
00148 }
00149 }
00150
00151
00152
00153
00154 void MeshBase::getFacetArray(int cellDim, int cellLID, int facetDim,
00155 Array<int>& facetLIDs,
00156 Array<int>& facetOrientations) const
00157 {
00158 int nf = numFacets(cellDim, cellLID, facetDim);
00159 facetLIDs.resize(nf);
00160 facetOrientations.resize(nf);
00161 for (int f=0; f<nf; f++)
00162 {
00163 facetLIDs[f] = facetLID(cellDim, cellLID, facetDim, f,
00164 facetOrientations[f]);
00165 }
00166 }
00167
00168
00169 void MeshBase::getLabels(int cellDim, const Array<int>& cellLID,
00170 Array<int>& labels) const
00171 {
00172 labels.resize(cellLID.size());
00173 for (int i=0; i<cellLID.size(); i++) labels[i] = label(cellDim, cellLID[i]);
00174 }
00175
00176
00177 void MeshBase::getMaxCofacetLIDs(
00178 const Array<int>& cellLIDs,
00179 MaximalCofacetBatch& cofacets) const
00180 {
00181 TEUCHOS_TEST_FOR_EXCEPT(cellLIDs.size()==0);
00182 int d = spatialDim() - 1;
00183 int nc = numMaxCofacets(d, cellLIDs[0]);
00184
00185 cofacets.reset(cellLIDs.size(), nc);
00186
00187 for (int i=0; i<cellLIDs.size(); i++)
00188 {
00189 int f1;
00190 int cofacet1 = maxCofacetLID(d, cellLIDs[i], 0, f1);
00191 if (nc==1) cofacets.addSingleCofacet(i, cofacet1, f1);
00192 else
00193 {
00194 int f2;
00195 int cofacet2 = maxCofacetLID(d, cellLIDs[i], 1, f2);
00196 cofacets.addTwoCofacets(i, cofacet1, f1, cofacet2, f2);
00197 }
00198 }
00199 }
00200
00201 void MeshBase::getLIDsForLabel(int cellDim, int labelToCheck,
00202 Array<int>& cellLIDs) const
00203 {
00204 cellLIDs.resize(0);
00205 int nc = numCells(cellDim);
00206
00207 for (int c=0; c<nc; c++)
00208 {
00209 int lc = label(cellDim, c);
00210 if (lc==labelToCheck) cellLIDs.append(c);
00211 }
00212 }
00213
00214 Set<int> MeshBase::getAllLabelsForDimension(int cellDim) const
00215 {
00216 int nc = numCells(cellDim);
00217 Set<int> rtn;
00218
00219 for (int c=0; c<nc; c++)
00220 {
00221 rtn.put(label(cellDim, c));
00222 }
00223 return rtn;
00224 }
00225
00226
00227
00228 bool MeshBase::hasSpecialWeight(int dim, int cellLID) const {
00229 if (specialWeights_[dim-1].containsKey(cellLID))
00230 return ((specialWeights_[dim-1].get(cellLID).size() > 0));
00231 else
00232 return false;
00233 }
00234
00235 void MeshBase::setSpecialWeight(int dim, int cellLID, Array<double>& w) const {
00236 specialWeights_[dim-1].put(cellLID,w);
00237 }
00238
00239 void MeshBase::getSpecialWeight(int dim, int cellLID, Array<double>& w) const {
00240 w = specialWeights_[dim-1].get(cellLID);
00241 }
00242
00243
00244 void MeshBase::flushSpecialWeights() const {
00245
00246 for (int d = 0 ; d < dim_ ; d++){
00247
00248 specialWeights_[d] = Sundance::Map< int , Array<double> >();
00249 }
00250 validWeights_ = true;
00251 }
00252
00253
00254
00255 void MeshBase::flushCurvePoints() const {
00256
00257
00258 curveID_to_ArrayIndex_ = Sundance::Map< int , int >();
00259 nrCurvesForIntegral_ = 0;
00260 curvePoints_Are_Valid_ = true;
00261 for (int c = 0 ; c < curvePoints_.size() ; c++ ){
00262
00263 curvePoints_[c] = Sundance::Map< int , Array<Point> >();
00264 curveDerivative_[c] = Sundance::Map< int , Array<Point> >();
00265 curveNormal_[c] = Sundance::Map< int , Array<Point> >();
00266 }
00267 }
00268
00269 bool MeshBase::hasCurvePoints(int maxCellLID , int curveID) const {
00270 if (curvePoints_[mapCurveID_to_Index(curveID)].containsKey(maxCellLID))
00271 return ( curvePoints_[mapCurveID_to_Index(curveID)].get(maxCellLID).size() > 0 );
00272 else
00273 return false;
00274 }
00275
00276 void MeshBase::setCurvePoints(int maxCellLID, int curveID ,
00277 Array<Point>& points , Array<Point>& derivs , Array<Point>& normals) const {
00278
00279 Tabs tabs;
00280 int verbo = 0;
00281 SUNDANCE_MSG3(verbo, tabs << "MeshBase::setCurvePoints , nr:" << nrCurvesForIntegral_);
00282 curvePoints_[mapCurveID_to_Index(curveID)].put( maxCellLID , points );
00283 curveDerivative_[mapCurveID_to_Index(curveID)].put( maxCellLID , derivs );
00284 curveNormal_[mapCurveID_to_Index(curveID)].put( maxCellLID , normals );
00285
00286 }
00287
00288 void MeshBase::getCurvePoints(int maxCellLID, int curveID ,
00289 Array<Point>& points , Array<Point>& derivs , Array<Point>& normals) const {
00290
00291 Tabs tabs;
00292 int verbo = 0;
00293 SUNDANCE_MSG3(verbo, tabs << "MeshBase::getCurvePoints , nr:" << nrCurvesForIntegral_);
00294 points = curvePoints_[mapCurveID_to_Index(curveID)].get( maxCellLID );
00295 derivs = curveDerivative_[mapCurveID_to_Index(curveID)].get( maxCellLID );
00296 normals = curveNormal_[mapCurveID_to_Index(curveID)].get( maxCellLID );
00297
00298 }
00299
00300 int MeshBase::mapCurveID_to_Index(int curveID) const {
00301
00302 Tabs tabs;
00303 int verbo = 0;
00304
00305 SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index curveID:" << curveID);
00306 if (curveID_to_ArrayIndex_.containsKey(curveID)){
00307 SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index value found for curveID:" << curveID << " ret:" << curveID_to_ArrayIndex_.get(curveID));
00308 return curveID_to_ArrayIndex_.get(curveID);
00309 } else {
00310 SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index create new :" << nrCurvesForIntegral_);
00311 curveID_to_ArrayIndex_.put( curveID , nrCurvesForIntegral_ );
00312 SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index , increment ");
00313 nrCurvesForIntegral_++;
00314 curvePoints_.resize(nrCurvesForIntegral_);
00315 curveDerivative_.resize(nrCurvesForIntegral_);
00316 curveNormal_.resize(nrCurvesForIntegral_);
00317 SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index create new :" << nrCurvesForIntegral_);
00318 return nrCurvesForIntegral_-1;
00319 }
00320
00321 }