From 63e4259d639884cda33983c8a48c75b554407047 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Tue, 10 Dec 2013 20:57:41 +0000
Subject: [PATCH] Search functions in OctreePost can now be given the node
coordinates of an element.
In that case we search for all the matches in the octree, and we try to return amongst those the element that matches the input coordinates.
We use this in MathEval so we can reliably make operations on discontinuous datasets (e.g. edge element vector fields), on identical meshes that were possibly reordered.
---
Plugin/MathEval.cpp | 9 ++-
Post/OctreePost.cpp | 171 +++++++++++++++++++++++++++++++++-----------
Post/OctreePost.h | 29 +++++---
3 files changed, 153 insertions(+), 56 deletions(-)
diff --git a/Plugin/MathEval.cpp b/Plugin/MathEval.cpp
index 7fbf6a0f43..277d1f4414 100644
--- a/Plugin/MathEval.cpp
+++ b/Plugin/MathEval.cpp
@@ -222,9 +222,12 @@ PView *GMSH_MathEvalPlugin::execute(PView *view)
values[0] = x[nod]; values[1] = y[nod]; values[2] = z[nod];
if(otherData){
if(octree){
- if(!octree->searchScalar(x[nod], y[nod], z[nod], &w[0], step2))
- if(!octree->searchVector(x[nod], y[nod], z[nod], &w[0], step2))
- octree->searchTensor(x[nod], y[nod], z[nod], &w[0], step2);
+ if(!octree->searchScalar(x[nod], y[nod], z[nod], &w[0], step2,
+ 0, numNodes, &x[0], &y[0], &z[0]))
+ if(!octree->searchVector(x[nod], y[nod], z[nod], &w[0], step2,
+ 0, numNodes, &x[0], &y[0], &z[0]))
+ octree->searchTensor(x[nod], y[nod], z[nod], &w[0], step2,
+ 0, numNodes, &x[0], &y[0], &z[0]);
}
else
for(int comp = 0; comp < otherNumComp; comp++)
diff --git a/Post/OctreePost.cpp b/Post/OctreePost.cpp
index 003182c8d8..c3039de8ab 100644
--- a/Post/OctreePost.cpp
+++ b/Post/OctreePost.cpp
@@ -364,6 +364,63 @@ void OctreePost::_create(PViewData *data)
}
}
+static void *getElement(double P[3], Octree *octree, int nbNod,
+ int qn, double *qx, double *qy, double *qz)
+{
+ if(qn && qx && qy && qz){
+ std::vector<void*> v;
+ Octree_SearchAll(P, octree, &v);
+ if(nbNod == qn){
+ // try to use the value from the same geometrical element as the one
+ // provided in qx/y/z
+ double eps = CTX::instance()->geom.tolerance;
+ for(unsigned int i = 0; i < v.size(); i++){
+ double *X = (double*)v[i], *Y = &X[qn], *Z = &X[2 * qn];
+ bool ok = true;
+ for(int j = 0; j < qn; j++){
+ ok &= (fabs(X[j] - qx[j]) < eps &&
+ fabs(Y[j] - qy[j]) < eps &&
+ fabs(Z[j] - qz[j]) < eps);
+ }
+ if(ok) return v[i];
+ }
+ }
+ if(v.size()) return v[0];
+ }
+ else{
+ return Octree_Search(P, octree);
+ }
+ return 0;
+}
+
+static MElement *getElement(double P[3], GModel *m,
+ int qn, double *qx, double *qy, double *qz)
+{
+ SPoint3 pt(P);
+ if(qn && qx && qy && qz){
+ // try to use the value from the same geometrical element as the one
+ // provided in qx/y/z
+ double eps = CTX::instance()->geom.tolerance;
+ std::vector<MElement*> elements = m->getMeshElementsByCoord(pt);
+ for(unsigned int i = 0; i < elements.size(); i++){
+ if(qn == elements[i]->getNumVertices()){
+ bool ok = true;
+ for(int j = 0; j < qn; j++){
+ MVertex *v = elements[i]->getVertex(j);
+ ok &= (fabs(v->x() - qx[j]) < eps &&
+ fabs(v->y() - qy[j]) < eps &&
+ fabs(v->z() - qz[j]) < eps);
+ }
+ }
+ }
+ if(elements.size()) return elements[0];
+ }
+ else{
+ return m->getMeshElementByCoord(pt);
+ }
+ return 0;
+}
+
bool OctreePost::_getValue(void *in, int dim, int nbNod, int nbComp,
double P[3], int step, double *values,
double *elementSize)
@@ -439,7 +496,8 @@ bool OctreePost::_getValue(void *in, int nbComp, double P[3], int timestep,
}
bool OctreePost::searchScalar(double x, double y, double z, double *values,
- int step, double *size)
+ int step, double *size,
+ int qn, double *qx, double *qy, double *qz)
{
double P[3] = {x, y, z};
@@ -455,20 +513,28 @@ bool OctreePost::searchScalar(double x, double y, double z, double *values,
values[0] = 0.0;
if(_theViewDataList){
- if(_getValue(Octree_Search(P, _SS), 3, 4, 1, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _SH), 3, 8, 1, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _SI), 3, 6, 1, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _SY), 3, 5, 1, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _ST), 2, 3, 1, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _SQ), 2, 4, 1, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _SL), 1, 2, 1, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _SPP), 0, 1, 1, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _SS, 4, qn, qx, qy, qz),
+ 3, 4, 1, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _SH, 8, qn, qx, qy, qz),
+ 3, 8, 1, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _SI, 6, qn, qx, qy, qz),
+ 3, 6, 1, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _SY, 5, qn, qx, qy, qz),
+ 3, 5, 1, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _ST, 3, qn, qx, qy, qz),
+ 2, 3, 1, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _SQ, 4, qn, qx, qy, qz),
+ 2, 4, 1, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _SL, 2, qn, qx, qy, qz),
+ 1, 2, 1, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _SPP, 1, qn, qx, qy, qz),
+ 0, 1, 1, P, step, values, size)) return true;
}
else if(_theViewDataGModel){
GModel *m = _theViewDataGModel->getModel((step < 0) ? 0 : step);
if(m){
- SPoint3 pt(P);
- if(_getValue(m->getMeshElementByCoord(pt), 1, P, step, values, size)) return true;
+ if(_getValue(getElement(P, m, qn, qx, qy, qz),
+ 1, P, step, values, size)) return true;
}
}
@@ -476,15 +542,16 @@ bool OctreePost::searchScalar(double x, double y, double z, double *values,
}
bool OctreePost::searchScalarWithTol(double x, double y, double z, double *values,
- int step, double *size, double tol)
+ int step, double *size, double tol,
+ int qn, double *qx, double *qy, double *qz)
{
- bool a = searchScalar(x, y, z, values, step, size);
+ bool a = searchScalar(x, y, z, values, step, size, qn, qx, qy, qz);
if(!a && tol != 0.){
double oldtol1 = element::getTolerance();
double oldtol2 = MElement::getTolerance();
element::setTolerance(tol);
MElement::setTolerance(tol);
- a = searchScalar(x, y, z, values, step, size);
+ a = searchScalar(x, y, z, values, step, size, qn, qx, qy, qz);
element::setTolerance(oldtol1);
MElement::setTolerance(oldtol2);
}
@@ -492,7 +559,8 @@ bool OctreePost::searchScalarWithTol(double x, double y, double z, double *value
}
bool OctreePost::searchVector(double x, double y, double z, double *values,
- int step, double *size)
+ int step, double *size,
+ int qn, double *qx, double *qy, double *qz)
{
double P[3] = {x, y, z};
@@ -508,20 +576,28 @@ bool OctreePost::searchVector(double x, double y, double z, double *values,
values[i] = 0.0;
if(_theViewDataList){
- if(_getValue(Octree_Search(P, _VS), 3, 4, 3, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _VH), 3, 8, 3, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _VI), 3, 6, 3, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _VY), 3, 5, 3, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _VT), 2, 3, 3, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _VQ), 2, 4, 3, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _VL), 1, 2, 3, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _VPP), 0, 1, 3, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _VS, 4, qn, qx, qy, qz),
+ 3, 4, 3, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _VH, 8, qn, qx, qy, qz),
+ 3, 8, 3, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _VI, 6, qn, qx, qy, qz),
+ 3, 6, 3, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _VY, 5, qn, qx, qy, qz),
+ 3, 5, 3, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _VT, 3, qn, qx, qy, qz),
+ 2, 3, 3, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _VQ, 4, qn, qx, qy, qz),
+ 2, 4, 3, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _VL, 2, qn, qx, qy, qz),
+ 1, 2, 3, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _VPP, 1, qn, qx, qy, qz),
+ 0, 1, 3, P, step, values, size)) return true;
}
else if(_theViewDataGModel){
GModel *m = _theViewDataGModel->getModel((step < 0) ? 0 : step);
if(m){
- SPoint3 pt(P);
- if(_getValue(m->getMeshElementByCoord(pt), 3, P, step, values, size)) return true;
+ if(_getValue(getElement(P, m, qn, qx, qy, qz),
+ 3, P, step, values, size)) return true;
}
}
@@ -529,15 +605,16 @@ bool OctreePost::searchVector(double x, double y, double z, double *values,
}
bool OctreePost::searchVectorWithTol(double x, double y, double z, double *values,
- int step, double *size, double tol)
+ int step, double *size, double tol,
+ int qn, double *qx, double *qy, double *qz)
{
- bool a = searchVector(x, y, z, values, step, size);
+ bool a = searchVector(x, y, z, values, step, size, qn, qx, qy, qz);
if(!a && tol != 0.){
double oldtol1 = element::getTolerance();
double oldtol2 = MElement::getTolerance();
element::setTolerance(tol);
MElement::setTolerance(tol);
- a = searchVector(x, y, z, values, step, size);
+ a = searchVector(x, y, z, values, step, size, qn, qx, qy, qz);
element::setTolerance(oldtol1);
MElement::setTolerance(oldtol2);
}
@@ -545,7 +622,8 @@ bool OctreePost::searchVectorWithTol(double x, double y, double z, double *value
}
bool OctreePost::searchTensor(double x, double y, double z, double *values,
- int step, double *size)
+ int step, double *size,
+ int qn, double *qx, double *qy, double *qz)
{
double P[3] = {x, y, z};
@@ -561,20 +639,28 @@ bool OctreePost::searchTensor(double x, double y, double z, double *values,
values[i] = 0.0;
if(_theViewDataList){
- if(_getValue(Octree_Search(P, _TS), 3, 4, 9, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _TH), 3, 8, 9, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _TI), 3, 6, 9, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _TY), 3, 5, 9, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _TT), 2, 3, 9, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _TQ), 2, 4, 9, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _TL), 1, 2, 9, P, step, values, size)) return true;
- if(_getValue(Octree_Search(P, _TPP), 0, 1, 9, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _TS, 4, qn, qx, qy, qz),
+ 3, 4, 9, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _TH, 8, qn, qx, qy, qz),
+ 3, 8, 9, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _TI, 6, qn, qx, qy, qz),
+ 3, 6, 9, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _TY, 5, qn, qx, qy, qz),
+ 3, 5, 9, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _TT, 3, qn, qx, qy, qz),
+ 2, 3, 9, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _TQ, 4, qn, qx, qy, qz),
+ 2, 4, 9, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _TL, 2, qn, qx, qy, qz),
+ 1, 2, 9, P, step, values, size)) return true;
+ if(_getValue(getElement(P, _TPP, 1, qn, qx, qy, qz),
+ 0, 1, 9, P, step, values, size)) return true;
}
else if(_theViewDataGModel){
GModel *m = _theViewDataGModel->getModel((step < 0) ? 0 : step);
if(m){
- SPoint3 pt(P);
- if(_getValue(m->getMeshElementByCoord(pt), 9, P, step, values, size)) return true;
+ if(_getValue(getElement(P, m, qn, qx, qy, qz),
+ 9, P, step, values, size)) return true;
}
}
@@ -582,15 +668,16 @@ bool OctreePost::searchTensor(double x, double y, double z, double *values,
}
bool OctreePost::searchTensorWithTol(double x, double y, double z, double *values,
- int step, double *size, double tol)
+ int step, double *size, double tol,
+ int qn, double *qx, double *qy, double *qz)
{
- bool a = searchTensor(x, y, z, values, step, size);
+ bool a = searchTensor(x, y, z, values, step, size, qn, qx, qy, qz);
if(!a && tol != 0.){
double oldtol1 = element::getTolerance();
double oldtol2 = MElement::getTolerance();
element::setTolerance(tol);
MElement::setTolerance(tol);
- a = searchTensor(x, y, z, values, step, size);
+ a = searchTensor(x, y, z, values, step, size, qn, qx, qy, qz);
element::setTolerance(oldtol1);
MElement::setTolerance(oldtol2);
}
diff --git a/Post/OctreePost.h b/Post/OctreePost.h
index 7e9993f5a9..e8bb378b5b 100644
--- a/Post/OctreePost.h
+++ b/Post/OctreePost.h
@@ -36,23 +36,30 @@ class OctreePost
OctreePost(PView *v);
OctreePost(PViewData *data);
~OctreePost();
- // 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 to a different value than
- // -1.
+ // 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
+ // to a different value than -1. If qn is given, n node coordinates stored in
+ // qx/y/z are used to select which element is used to interpolate (if the
+ // query returned more than one)
bool searchScalar(double x, double y, double z, double *values,
- int step=-1, double *size=0);
+ int step=-1, double *size=0,
+ int qn=0, double *qx=0, double *qy=0, double *qz=0);
bool searchScalarWithTol(double x, double y, double z, double *values,
- int step=-1, double *size=0, double tol=1.e-2);
+ int step=-1, double *size=0, double tol=1.e-2,
+ int qn=0, double *qx=0, double *qy=0, double *qz=0);
bool searchVector(double x, double y, double z, double *values,
- int step=-1, double *size=0);
+ int step=-1, double *size=0,
+ int qn=0, double *qx=0, double *qy=0, double *qz=0);
bool searchVectorWithTol(double x, double y, double z, double *values,
- int step=-1, double *size=0, double tol=1.e-2);
+ int step=-1, double *size=0, double tol=1.e-2,
+ int qn=0, double *qx=0, double *qy=0, double *qz=0);
bool searchTensor(double x, double y, double z, double *values,
- int step=-1, double *size=0);
+ int step=-1, double *size=0,
+ int qn=0, double *qx=0, double *qy=0, double *qz=0);
bool searchTensorWithTol(double x, double y, double z, double *values,
- int step=-1, double *size=0, double tol=1.e-2);
+ int step=-1, double *size=0, double tol=1.e-2,
+ int qn=0, double *qx=0, double *qy=0, double *qz=0);
};
#endif
--
GitLab