00001 #include "PDEOptPointData.hpp"
00002 #include "PlayaTabs.hpp"
00003
00004 namespace Sundance
00005 {
00006 using namespace Teuchos;
00007 using namespace Playa;
00008
00009 PointData::PointData(const Array<Point>& locations,
00010 const Array<double>& values,
00011 const double& tol)
00012 : sensorVals_(), sensorLocations_()
00013 {
00014 init(locations, values, tol);
00015 }
00016
00017 PointData::PointData(const XMLObject& xml,
00018 const Mesh& mesh)
00019 : sensorVals_(), sensorLocations_()
00020 {
00021 TEUCHOS_TEST_FOR_EXCEPTION(xml.getTag() != "PointData", RuntimeError,
00022 "expected tag PointData, found " << xml.getTag());
00023
00024 double tol = xml.getRequiredDouble("tol");
00025
00026 Array<Point> locations;
00027 Array<double> values;
00028
00029 for (int i=0; i<xml.numChildren(); i++)
00030 {
00031 const XMLObject& pt = xml.getChild(i);
00032 values.append( pt.getRequiredDouble("value") );
00033 const string& ptStr = pt.getRequired("pt");
00034 Array<string> tok = StrUtils::stringTokenizer(ptStr);
00035 Point p;
00036 if (tok.size()==1)
00037 {
00038 p = Point(StrUtils::atof(tok[0]));
00039 }
00040 else if (tok.size()==2)
00041 {
00042 p = Point(StrUtils::atof(tok[0]),
00043 StrUtils::atof(tok[1]));
00044 }
00045 else if (tok.size()==3)
00046 {
00047 p = Point(StrUtils::atof(tok[0]),
00048 StrUtils::atof(tok[1]),
00049 StrUtils::atof(tok[2]));
00050 }
00051 locations.append(p);
00052 }
00053
00054 locations = snapToMesh(mesh, locations);
00055 init(locations, values, tol);
00056 }
00057
00058
00059 void PointData::init(const Array<Point>& locations,
00060 const Array<double>& values,
00061 const double& tol)
00062 {
00063 TEUCHOS_TEST_FOR_EXCEPTION(locations.size() != values.size(), RuntimeError,
00064 "inconsistent measurement data: num locations = "
00065 << locations.size() << " but num readings = "
00066 << values.size());
00067
00068 TEUCHOS_TEST_FOR_EXCEPTION(locations.size() < 1, RuntimeError,
00069 "Empty data set in PointData ctor");
00070
00071
00072 int dim = locations[0].dim();
00073 for (int i=0; i<locations.size(); i++)
00074 {
00075 TEUCHOS_TEST_FOR_EXCEPTION(dim != locations[i].dim(), RuntimeError,
00076 "inconsistent point dimensions in PointData ctor. "
00077 "points are " << locations);
00078 }
00079
00080 CellFilter allPoints = new DimensionalCellFilter(0);
00081 sensorLocations_
00082 = allPoints.subset(new PointDataCellPredicateFunctor(locations, tol));
00083
00084 RefCountPtr<UserDefFunctor> op
00085 = rcp(new PointDataExprFunctor(locations, values, tol));
00086
00087
00088
00089 Expr x = new CoordExpr(0);
00090 Expr y = new CoordExpr(1);
00091 Expr z = new CoordExpr(2);
00092
00093 Expr coords;
00094 if (dim==1) coords = x;
00095 else if (dim==2) coords = List(x, y);
00096 else coords = List(x, y, z);
00097
00098 sensorVals_ = new UserDefOp(coords, op);
00099 }
00100
00101
00102 Array<Point> PointData::snapToMesh(const Mesh& mesh,
00103 const Array<Point>& locations)
00104 {
00105 Array<Point> rtn(locations.size());
00106
00107 for (int i=0; i<locations.size(); i++)
00108 {
00109 rtn[i] = nearestMeshPoint(mesh, locations[i]);
00110 }
00111 return rtn;
00112 }
00113
00114 Point PointData::nearestMeshPoint(const Mesh& mesh, const Point& x)
00115 {
00116 double d2Min = 1.0e100;
00117 int iPt = 0;
00118 for (int i=0; i<mesh.numCells(0); i++)
00119 {
00120 const Point& p = mesh.nodePosition(i);
00121 Point dx = x - p;
00122 double d2 = dx*dx;
00123 if (d2 < d2Min)
00124 {
00125 d2Min = d2;
00126 iPt = i;
00127 }
00128 }
00129 return mesh.nodePosition(iPt);
00130 }
00131
00132
00133
00134 PointDataExprFunctor::PointDataExprFunctor(const Array<Point>& locations,
00135 const Array<double>& values,
00136 const double& tol)
00137 : PointwiseUserDefFunctor0("pointData", locations[0].dim(), 1), pointToValueMap_(tol),
00138 dim_(locations[0].dim())
00139 {
00140 Tabs tabs;
00141 for (int i=0; i<locations.size(); i++)
00142 {
00143 TEUCHOS_TEST_FOR_EXCEPTION(pointToValueMap_.find(locations[i]) != pointToValueMap_.end(),
00144 RuntimeError,
00145 "Data set contains duplicate point " << locations[i]
00146 << " to within tolerance "
00147 << tol << ". Points are "
00148 << locations);
00149 pointToValueMap_[locations[i]] = values[i];
00150
00151 TEUCHOS_TEST_FOR_EXCEPTION(pointToValueMap_.find(locations[i]) == pointToValueMap_.end(),
00152 InternalError,
00153 "Bad map integrity in PointDataExprFunctor ctor: point "
00154 << locations[i] << " could not be recovered from map");
00155 }
00156 }
00157
00158 void PointDataExprFunctor::eval0(const double* in, double* out) const
00159 {
00160 Point p;
00161 if (dim_==1) p = Point(in[0]);
00162 else if (dim_==2) p = Point(in[0], in[1]);
00163 else p = Point(in[0], in[1], in[2]);
00164
00165 typedef std::map<Point, double, SloppyPointComparitor>::const_iterator iter;
00166 iter pos = pointToValueMap_.find(p);
00167 TEUCHOS_TEST_FOR_EXCEPTION(pos == pointToValueMap_.end(), RuntimeError,
00168 "Point " << p << " not found in list of data values");
00169
00170
00171 Tabs tabs;
00172
00173 *out = pos->second;
00174 }
00175
00176
00177 PointDataCellPredicateFunctor
00178 ::PointDataCellPredicateFunctor(const Array<Point>& locations,
00179 const double& tol)
00180 : pointSet_(SloppyPointComparitor(tol))
00181 {
00182 for (int i=0; i<locations.size(); i++)
00183 {
00184 pointSet_.insert(locations[i]);
00185 }
00186 }
00187
00188 bool PointDataCellPredicateFunctor::operator()(const Point& x) const
00189 {
00190 Tabs tabs;
00191 bool rtn = pointSet_.find(x) != pointSet_.end();
00192 return rtn;
00193 }
00194
00195
00196 bool SloppyPointComparitor::operator()(const Point& p1, const Point& p2) const
00197 {
00198
00199 if (p1.dim() < p2.dim()) return true;
00200 if (p2.dim() < p1.dim()) return false;
00201
00202
00203
00204 for (int i=0; i<p1.dim(); i++)
00205 {
00206 if (p1[i] < p2[i] - tol_) return true;
00207 if (p1[i] > p2[i] + tol_) return false;
00208 }
00209 return false;
00210 }
00211
00212 }