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 "SundanceAToCDensitySampler.hpp"
00032 #include "SundanceOut.hpp"
00033 #include "PlayaTabs.hpp"
00034 #include "SundanceGeomUtils.hpp"
00035 #include "SundanceLagrange.hpp"
00036 #include "SundanceDiscreteFuncElement.hpp"
00037 #include <queue>
00038
00039 using namespace Sundance;
00040 using namespace Sundance;
00041 using namespace Sundance;
00042 using namespace Sundance;
00043 using namespace Sundance;
00044 using namespace Sundance;
00045 using namespace Sundance;
00046 using namespace Teuchos;
00047 using namespace Playa;
00048
00049
00050 static Time& densitySamplingTimer()
00051 {
00052 static RCP<Time> rtn
00053 = TimeMonitor::getNewTimer("density sampling");
00054 return *rtn;
00055 }
00056
00057 AToCDensitySampler::AToCDensitySampler(const AToCPointLocator& locator,
00058 const VectorType<double>& vecType)
00059 : discSpace_(locator.mesh(), new Lagrange(0), locator.subdomain(), vecType),
00060 dim_(locator.mesh().spatialDim()),
00061 mesh_(locator.mesh()),
00062 elemToVecIndexMap_(),
00063 elemWeights_(new DiscreteFunction(discSpace_, 0.0)),
00064 elemWeightVec_(),
00065 locator_(locator),
00066 isAxisymmetric_(false),
00067 origin_(),
00068 axis_()
00069 {
00070 init();
00071 }
00072
00073
00074 AToCDensitySampler::AToCDensitySampler(const AToCPointLocator& locator,
00075 const std::vector<double>& origin,
00076 const std::vector<double>& rotationalAxis,
00077 const VectorType<double>& vecType)
00078 : discSpace_(locator.mesh(), new Lagrange(0), locator.subdomain(), vecType),
00079 dim_(locator.mesh().spatialDim()),
00080 mesh_(locator.mesh()),
00081 elemToVecIndexMap_(),
00082 elemWeights_(new DiscreteFunction(discSpace_, 0.0)),
00083 elemWeightVec_(),
00084 locator_(locator),
00085 isAxisymmetric_(true),
00086 origin_(vec2point(origin)),
00087 axis_(normPoint(vec2point(rotationalAxis)))
00088 {
00089 init();
00090 }
00091
00092
00093
00094 void AToCDensitySampler::init()
00095 {
00096 const CellFilter& domain = discSpace_.cellFilters(0);
00097
00098 elemWeightVec_ = DiscreteFunction::discFunc(elemWeights_)->getVector();
00099
00100 elemToVecIndexMap_ = rcp(new Array<int>(mesh_.numCells(dim_), -1));
00101
00102 Array<int>& a = *elemToVecIndexMap_;
00103
00104 CellSet cells = domain.getCells(mesh_);
00105
00106 Array<int> cellLID;
00107 cellLID.reserve(mesh_.numCells(dim_));
00108
00109 for (CellIterator i=cells.begin(); i!=cells.end(); i++)
00110 {
00111 cellLID.append(*i);
00112 }
00113
00114 const RCP<DOFMapBase>& dofMap = discSpace_.map();
00115
00116 Set<int> funcs = makeSet(0);
00117 Array<Array<int> > dofs;
00118 Array<int> nNodes;
00119 dofMap->getDOFsForCellBatch(dim_, cellLID, funcs, dofs, nNodes,0);
00120
00121 const Array<int>& dofs0 = dofs[0];
00122 for (int c=0; c<cellLID.size(); c++)
00123 {
00124 int vecIndex = dofs0[c];
00125 int lid = cellLID[c];
00126 a[lid] = vecIndex;
00127 double vol = volume(mesh_, dim_, lid);
00128 if (isAxisymmetric_)
00129 {
00130 Point xCell = mesh_.centroid(dim_, lid) - origin_;
00131 double dPerp = ::sqrt(xCell*xCell - (xCell*axis_)*(xCell*axis_));
00132 vol = vol * dPerp;
00133 }
00134 elemWeightVec_[vecIndex] = vol;
00135 }
00136 }
00137
00138 Point AToCDensitySampler::vec2point(const std::vector<double>& x) const
00139 {
00140 if (x.size()==1) return Point(x[0]);
00141 else if (x.size()==2U) return Point(x[0], x[1]);
00142 else if (x.size()==3U) return Point(x[0], x[1], x[2]);
00143 TEUCHOS_TEST_FOR_EXCEPT(x.size() < 1 || x.size() > 3U);
00144 return Point();
00145 }
00146
00147 Point AToCDensitySampler::normPoint(const Point& x) const
00148 {
00149 return x/sqrt(x*x);
00150 }
00151
00152
00153 Expr AToCDensitySampler::sample(const std::vector<double>& positions,
00154 const double& particleWeight) const
00155 {
00156 TimeMonitor timer(densitySamplingTimer());
00157
00158 TEUCHOS_TEST_FOR_EXCEPTION(positions.size() % dim_ != 0, std::runtime_error,
00159 "vector of coordinates should by an integer multiple "
00160 "of the spatial dimension");
00161
00162 Expr rtn = new DiscreteFunction(discSpace_, 0.0);
00163 Vector<double> density = DiscreteFunction::discFunc(rtn)->getVector();
00164
00165 int nPts = positions.size() / dim_;
00166
00167 for (int i=0; i<nPts; i++)
00168 {
00169 const double* x = &(positions[dim_*i]);
00170
00171 int guess = locator_.guessCell(x);
00172
00173 TEUCHOS_TEST_FOR_EXCEPTION(guess < 0, std::runtime_error, "particle #" << i << " position="
00174 << AToCPointLocator::makePoint(dim_, x)
00175 << " is out of search grid");
00176
00177 int cellLID = locator_.findEnclosingCell(guess, x);
00178
00179 int vecIndex = (*elemToVecIndexMap_)[cellLID];
00180 double vol = elemWeightVec_[vecIndex];
00181 density[vecIndex] += particleWeight/vol;
00182 }
00183
00184 return rtn;
00185 }
00186
00187
00188 Expr AToCDensitySampler::resetCounts() const
00189 {
00190 Expr rtn = new DiscreteFunction(discSpace_, 0.0);
00191
00192 return rtn;
00193 }
00194
00195 void AToCDensitySampler::addToCounts(const std::vector<double>& positions,
00196 const double& particleWeight,
00197 Expr density) const
00198 {
00199 TimeMonitor timer(densitySamplingTimer());
00200
00201 TEUCHOS_TEST_FOR_EXCEPTION(positions.size() % dim_ != 0, std::runtime_error,
00202 "vector of coordinates should by an integer multiple "
00203 "of the spatial dimension");
00204
00205 Vector<double> vec = DiscreteFunction::discFunc(density)->getVector();
00206
00207 int nPts = positions.size() / dim_;
00208
00209 for (int i=0; i<nPts; i++)
00210 {
00211 const double* x = &(positions[dim_*i]);
00212
00213 int guess = locator_.guessCell(x);
00214
00215 TEUCHOS_TEST_FOR_EXCEPTION(guess < 0, std::runtime_error, "particle #" << i << " position="
00216 << AToCPointLocator::makePoint(dim_, x)
00217 << " is out of search grid");
00218
00219 int cellLID = locator_.findEnclosingCell(guess, x);
00220
00221 int vecIndex = (*elemToVecIndexMap_)[cellLID];
00222 double vol = elemWeightVec_[vecIndex];
00223 vec[vecIndex] += particleWeight/vol;
00224 }
00225 }
00226
00227