Skip to content

GitLab

Nektar
Nektar
Commit 36763032 authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files
Options

Merge branch 'fix/probefld-manifold' of /opt/gitlab/repositories/nektar

parents ebc2db7b fec5dddc
Hide whitespace changes
Inline Side-by-side
Showing with 138 additions and 54 deletions
library/MultiRegions/ExpList.cpp
View file @ 36763032
......@@ -1220,60 +1220,92 @@ namespace Nektar
Array<OneD, NekDouble> &locCoords,
NekDouble tol)
{
static int start = 0;
NekDouble resid, min_resid = NekConstants::kNekMinResidInit;
int min_elmt;
Array<OneD, NekDouble> min_locCoords(locCoords.num_elements());
NekDouble resid;
// start search at previous element or 0
for (int i = start; i < (*m_exp).size(); ++i)
if (GetNumElmts() == 0)
{
if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
tol, resid))
return -1;
}
// Manifold case (point may match multiple elements)
if (GetExp(0)->GetCoordim() > GetExp(0)->GetShapeDimension())
{
std::vector<std::pair<int,NekDouble> > elmtIdDist;
SpatialDomains::PointGeomSharedPtr v;
SpatialDomains::PointGeom w;
NekDouble x, y, z;
// Scan all elements and store those which may contain the point
for (int i = 0; i < (*m_exp).size(); ++i)
{
start = i;
return i;
if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
tol, resid))
{
v = m_graph->GetVertex((*m_exp)[i]->GetGeom()->GetVid(0));
w.SetX(gloCoords[0]);
w.SetY(gloCoords[1]);
w.SetZ(gloCoords[2]);
v->GetCoords(x,y,z);
elmtIdDist.push_back(std::pair<int, NekDouble>(i, v->dist(w)));
}
}
else
// Find nearest element
if (!elmtIdDist.empty())
{
if(resid < min_resid)
NekDouble min_d = elmtIdDist[0].first;
int min_id = elmtIdDist[0].second;
for (int i = 1; i < elmtIdDist.size(); ++i)
{
min_resid = resid;
min_elmt = i;
Vmath::Vcopy(locCoords.num_elements(), locCoords, 1,
min_locCoords,1);
if (elmtIdDist[i].second < min_d) {
min_d = elmtIdDist[i].second;
min_id = elmtIdDist[i].first;
}
}
return min_id;
}
else {
return -1;
}
}
for (int i = 0; i < start; ++i)
// non-embedded mesh (point can only match one element)
else
{
if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
tol,resid))
static int start = 0;
// restart search from last found value
for (int i = start; i < (*m_exp).size(); ++i)
{
start = i;
return i;
if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
tol, resid))
{
start = i;
return i;
}
}
else
for (int i = 0; i < start; ++i)
{
if(resid < min_resid)
if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
tol, resid))
{
min_resid = resid;
min_elmt = i;
Vmath::Vcopy(locCoords.num_elements(), locCoords, 1,
min_locCoords,1);
start = i;
return i;
}
}
}
std::string msg = "Failed to find point in element to tolerance of "
std::string msg = "Failed to find point in element to tolerance of "
+ boost::lexical_cast<std::string>(resid)
+ " using nearest point found";
WARNINGL0(true,msg.c_str());
WARNINGL0(true,msg.c_str());
Vmath::Vcopy(locCoords.num_elements(),min_locCoords,1,locCoords,1);
return -1;
return min_elmt;
}
}
......
utilities/PostProcessing/ProbeFld.cpp
View file @ 36763032
......@@ -12,14 +12,25 @@ using namespace Nektar;
int main(int argc, char *argv[])
{
int i,j;
bool file = false;
if(argc != 10)
if(argc == 4)
{
file = true;
}
else if(argc != 10)
{
fprintf(stderr,
"Usage: ProbeFld meshfile fieldfile N x0 y0 z0 dx dy dz\n");
fprintf(stderr,
" Probes N points along the line from (x0,y0,z0) to "
"(x0+dx, y0+dy, z0+dz)\n");
fprintf(stderr,
"ProbeFld meshfile fieldfile points.txt\n");
fprintf(stderr,
" Probes the solution at the points in the points.txt file.\n");
fprintf(stderr,
" Points are given as space-separated x y z on each line.\n");
exit(1);
}
......@@ -204,30 +215,71 @@ int main(int argc, char *argv[])
}
//----------------------------------------------
//----------------------------------------------
// Probe data fields
NekDouble N = atoi(argv[3]);
NekDouble x0 = atof(argv[4]);
NekDouble y0 = atof(argv[5]);
NekDouble z0 = atof(argv[6]);
NekDouble dx = atof(argv[7])/(N>1 ? (N-1) : 1);
NekDouble dy = atof(argv[8])/(N>1 ? (N-1) : 1);
NekDouble dz = atof(argv[9])/(N>1 ? (N-1) : 1);
Array<OneD, NekDouble> gloCoord(3,0.0);
if (file)
{
string line;
ifstream pts(argv[3]);
while (getline(pts, line))
{
stringstream ss(line);
ss >> gloCoord[0];
ss >> gloCoord[1];
ss >> gloCoord[2];
cout << gloCoord[0] << " " << gloCoord[1] << " " << gloCoord[2];
int ExpId = Exp[0]->GetExpIndex(gloCoord,NekConstants::kGeomFactorsTol);
for (int j = 0; j < nfields; ++j)
{
if (ExpId == -1)
{
cout << " -";
}
else
{
Array<OneD, NekDouble> phys(Exp[j]->GetPhys() + Exp[j]->GetPhys_Offset(ExpId));
cout << " " << Exp[j]->GetExp(ExpId)->PhysEvaluate(gloCoord, phys);
}
}
for (int i = 0; i < N; ++i)
cout << endl;
}
}
else
{
gloCoord[0] = x0 + i*dx;
gloCoord[1] = y0 + i*dy;
gloCoord[2] = z0 + i*dz;
cout << gloCoord[0] << " " << gloCoord[1] << " " << gloCoord[2];
int ExpId = Exp[0]->GetExpIndex(gloCoord,NekConstants::kGeomFactorsTol);
for (int j = 0; j < nfields; ++j)
//----------------------------------------------
// Probe data fields
NekDouble N = atoi(argv[3]);
NekDouble x0 = atof(argv[4]);
NekDouble y0 = atof(argv[5]);
NekDouble z0 = atof(argv[6]);
NekDouble dx = atof(argv[7])/(N>1 ? (N-1) : 1);
NekDouble dy = atof(argv[8])/(N>1 ? (N-1) : 1);
NekDouble dz = atof(argv[9])/(N>1 ? (N-1) : 1);
for (int i = 0; i < N; ++i)
{
Array<OneD, NekDouble> phys(Exp[j]->GetPhys() + Exp[j]->GetPhys_Offset(j));
cout << " " << Exp[j]->GetExp(ExpId)->PhysEvaluate(gloCoord, phys);
gloCoord[0] = x0 + i*dx;
gloCoord[1] = y0 + i*dy;
gloCoord[2] = z0 + i*dz;
cout << gloCoord[0] << " " << gloCoord[1] << " " << gloCoord[2];
int ExpId = Exp[0]->GetExpIndex(gloCoord,NekConstants::kGeomFactorsTol);
for (int j = 0; j < nfields; ++j)
{
if (ExpId == -1)
{
cout << " -";
}
else
{
Array<OneD, NekDouble> phys(Exp[j]->GetPhys() + Exp[j]->GetPhys_Offset(ExpId));
cout << " " << Exp[j]->GetExp(ExpId)->PhysEvaluate(gloCoord, phys);
}
}
cout << endl;
}
cout << endl;
}
//----------------------------------------------
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment