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 "SundanceCToAInterpolator.hpp"
00032 #include "SundanceOut.hpp"
00033 #include "PlayaTabs.hpp"
00034 #include "SundanceLagrange.hpp"
00035 #include "SundanceDiscreteFuncElement.hpp"
00036
00037 using namespace Sundance;
00038 using namespace Sundance;
00039 using namespace Sundance;
00040 using namespace Sundance;
00041 using namespace Sundance;
00042 using namespace Sundance;
00043 using namespace Sundance;
00044 using namespace Teuchos;
00045 using namespace Playa;
00046
00047
00048 static Time& particleInterpolationTimer()
00049 {
00050 static RCP<Time> rtn
00051 = TimeMonitor::getNewTimer("particle interpolation");
00052 return *rtn;
00053 }
00054
00055 CToAInterpolator::CToAInterpolator(const AToCPointLocator& locator,
00056 const Expr& field)
00057 : dim_(locator.mesh().spatialDim()),
00058 nFacets_(dim_+1),
00059 rangeDim_(-1),
00060 elemToVecValuesMap_(),
00061 locator_(locator)
00062 {
00063 updateField(field);
00064 }
00065
00066 void CToAInterpolator::updateField(const Expr& field)
00067 {
00068 int newRangeDim = field.size();
00069 if (newRangeDim != rangeDim_)
00070 {
00071 rangeDim_ = newRangeDim;
00072 elemToVecValuesMap_
00073 = rcp(new Array<double>(rangeDim_ * locator_.mesh().numCells(dim_) * nFacets_));
00074 }
00075
00076 int nCells = locator_.mesh().numCells(dim_);
00077
00078 const DiscreteFunction* df = DiscreteFunction::discFunc(field);
00079 TEUCHOS_TEST_FOR_EXCEPTION(df == 0, std::runtime_error,
00080 "discrete function expected in "
00081 "CToAInterpolator::updateField()");
00082
00083 Vector<double> vec = df->getVector();
00084
00085 Array<int> cellLID(nCells);
00086 for (int i=0; i<nCells; i++)
00087 {
00088 cellLID[i] = i;
00089 }
00090 Array<Array<int> > dofs;
00091 Array<int> nNodes;
00092 Set<int> funcs;
00093 for (int i=0; i<rangeDim_; i++) funcs.put(i);
00094
00095
00096 df->map()->getDOFsForCellBatch(dim_, cellLID, funcs, dofs, nNodes,0);
00097
00098 const Array<int>& dofs0 = dofs[0];
00099
00100 for (int c=0; c<cellLID.size(); c++)
00101 {
00102 for (int n=0; n<nFacets_; n++)
00103 {
00104 for (int f=0; f<rangeDim_; f++)
00105 {
00106 int dof = dofs0[(c*rangeDim_ + f)*nFacets_ + n];
00107 (*elemToVecValuesMap_)[(cellLID[c]*nFacets_ + n)*rangeDim_ + f]
00108 = vec[dof];
00109 }
00110 }
00111 }
00112 }
00113
00114
00115 void CToAInterpolator::interpolate(const Teuchos::Array<double>& positions,
00116 Teuchos::Array<double>& results) const
00117 {
00118 TimeMonitor timer(particleInterpolationTimer());
00119
00120 TEUCHOS_TEST_FOR_EXCEPTION(positions.size() % dim_ != 0, std::runtime_error,
00121 "vector of coordinates should by an integer multiple "
00122 "of the spatial dimension");
00123
00124 int nPts = positions.size() / dim_;
00125
00126 results.resize(rangeDim_ * nPts);
00127
00128 Array<double> xLocal(dim_);
00129
00130 for (int i=0; i<nPts; i++)
00131 {
00132 const double* x = &(positions[dim_*i]);
00133
00134 int guess = locator_.guessCell(x);
00135
00136 TEUCHOS_TEST_FOR_EXCEPTION(guess < 0, std::runtime_error, "particle position "
00137 << x << " out of search grid");
00138
00139
00140 int cellLID = locator_.findEnclosingCell(guess, x, &(xLocal[0]));
00141
00142 if (dim_==2)
00143 {
00144 double s = xLocal[0];
00145 double t = xLocal[1];
00146 Array<double> phi(nFacets_);
00147 phi[0] = 1.0-s-t;
00148 phi[1] = s;
00149 phi[2] = t;
00150 for (int f=0; f<rangeDim_; f++) results[rangeDim_*i+f] = 0.0;
00151 for (int n=0; n<nFacets_; n++)
00152 {
00153 for (int f=0; f<rangeDim_; f++)
00154 {
00155 results[rangeDim_*i+f] += phi[n]*(*elemToVecValuesMap_)[(cellLID*nFacets_ + n)*rangeDim_ + f];
00156 }
00157 }
00158 }
00159 }
00160
00161 }
00162