Skip to content
Snippets Groups Projects
Commit f82c25cc authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

functional "naive" implementation with a linear search

parent 4bf0ef6e
No related branches found
No related tags found
No related merge requests found
......@@ -243,6 +243,39 @@ bool PViewData::combineSpace(nameData &nd)
return false;
}
double PViewData::findClosestNode(double &xn, double &yn, double &zn, int step)
{
double dist2 = 1e200;
// if this gets used a lot we should create a kdtree to make efficient nearest
// neighbor lookups; moreover iterations on view data is currently not
// thread-safe (it uses a cache for the current element/node), which further
// motivates a better solution
#pragma omp critical
{
double x = xn, y = yn, z = zn;
if(step < 0) step = getFirstNonEmptyTimeStep();
for(int ent = 0; ent < getNumEntities(step); ent++) {
for(int ele = 0; ele < getNumElements(step, ent); ele++) {
int numNodes = getNumNodes(step, ent, ele);
for(int nod = 0; nod < numNodes; nod++) {
double xx, yy, zz;
getNode(step, ent, ele, nod, xx, yy, zz);
double d2 =
(x - xx) * (x - xx) + (y - yy) * (y - yy) + (z - zz) * (z - zz);
if(d2 < dist2) {
dist2 = d2;
xn = xx;
yn = yy;
zn = zz;
}
}
}
}
}
return sqrt(dist2);
}
bool PViewData::searchScalar(double x, double y, double z, double *values,
int step, double *size, int qn, double *qx,
double *qy, double *qz, bool grad, int dim)
......@@ -279,44 +312,13 @@ bool PViewData::searchScalarWithTol(double x, double y, double z,
bool ret = _octree->searchScalarWithTol(x, y, z, values, step, size, tol, qn, qx,
qy, qz, grad, dim);
if(!ret) {
// try a linear search in all elements and return the value at the closest
// node
double dist = 1e200;
int foundEnt = -1, foundEle = -1, foundNode = -1, foundNumComp = -1;
for(int ent = 0; ent < getNumEntities(step); ent++) {
for(int ele = 0; ele < getNumElements(step, ent); ele++) {
if(skipElement(step, ent, ele)) continue;
int numNodes = getNumNodes(step, ent, ele);
int numComp = getNumComponents(step, ent, ele);
for(int nod = 0; nod < numNodes; nod++) {
double xn, yn, zn;
getNode(step, ent, ele, nod, xn, yn, zn);
double d = std::sqrt((x - xn) * (x - xn) +
(y - yn) * (y - yn) +
(z - zn) * (z - zn));
if(d < dist) {
dist = d;
foundEnt = ent;
foundEle = ele;
foundNode = nod;
foundNumComp = numComp;
}
}
}
}
Msg::Info("Found closest node %d in ele %d", foundNode, foundEle);
if(grad) {
Msg::Warning("Gradient not coded yet for closest node probe");
return false;
}
if(foundEnt < 0 || foundEle < 0 || foundNode < 0 || foundNumComp != 1) {
Msg::Warning("Could not find closest node in probe");
return false;
}
getValue(step, foundEnt, foundEle, foundNode, 0, values[0]);
return true;
double xn = x, yn = y, zn = z;
double d = findClosestNode(xn, yn, zn, step);
ret = _octree->searchScalarWithTol(xn, yn, zn, values, step, size, tol, qn, qx,
qy, qz, grad, dim);
if(ret && d > tol)
Msg::Info("Returning value at closest node (distance = %g)", d);
}
return ret;
}
......@@ -353,8 +355,17 @@ bool PViewData::searchVectorWithTol(double x, double y, double z,
_octree = new OctreePost(this);
}
}
return _octree->searchVectorWithTol(x, y, z, values, step, size, tol, qn, qx,
bool ret = _octree->searchVectorWithTol(x, y, z, values, step, size, tol, qn, qx,
qy, qz, grad, dim);
if(!ret) {
double xn = x, yn = y, zn = z;
double d = findClosestNode(xn, yn, zn, step);
ret = _octree->searchVectorWithTol(xn, yn, zn, values, step, size, tol, qn, qx,
qy, qz, grad, dim);
if(ret && d > tol)
Msg::Info("Returning value at closest node (distance = %g)", d);
}
return ret;
}
bool PViewData::searchTensor(double x, double y, double z, double *values,
......@@ -390,6 +401,15 @@ bool PViewData::searchTensorWithTol(double x, double y, double z,
_octree = new OctreePost(this);
}
}
return _octree->searchTensorWithTol(x, y, z, values, step, size, tol, qn, qx,
bool ret = _octree->searchTensorWithTol(x, y, z, values, step, size, tol, qn, qx,
qy, qz, grad, dim);
if(!ret) {
double xn = x, yn = y, zn = z;
double d = findClosestNode(xn, yn, zn, step);
ret = _octree->searchTensorWithTol(xn, yn, zn, values, step, size, tol, qn, qx,
qy, qz, grad, dim);
if(ret && d > tol)
Msg::Info("Returning value at closest node (distance = %g)", d);
}
return ret;
}
......@@ -291,6 +291,10 @@ public:
// get MElement (if view supports it)
virtual MElement *getElement(int step, int entity, int element);
// find coordinates of closest node to point (xn, yn, zn); currently performs
// a simple linear search - we might want to use a kdtree instead
double findClosestNode(double &xn, double &yn, double &zn, int step);
// search for the value of the View at point x, y, z. Values are interpolated
// using standard first order shape functions in the post element. If several
// time steps are present, they are all interpolated unless time step is set
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment