Commit 7684cc76 authored by pugmire's avatar pugmire

Code cleanup and a bug fix in the AB integrator.


git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@11491 18c085ea-50e0-402c-830e-de6fd14e8384
parent cd085afa
......@@ -1943,9 +1943,12 @@ avtStreamlineFilter::DomainToRank(DomainType &domain)
// Dave Pugmire, Tue Mar 23 11:11:11 EDT 2010
// Moved zone-to-node centering to the streamline plot.
//
// Dave Pugmire, Wed May 26 13:48:24 EDT 2010
// New return type from avtStreamline::Advance()
//
// ****************************************************************************
avtIVPSolver::Result
void
avtStreamlineFilter::IntegrateDomain(avtStreamlineWrapper *slSeg,
vtkDataSet *ds,
double *extents,
......@@ -2004,7 +2007,7 @@ avtStreamlineFilter::IntegrateDomain(avtStreamlineWrapper *slSeg,
//slSeg->Debug();
int numSteps = slSeg->sl->size();
avtIVPSolver::Result result;
avtStreamline::Result result;
// When restarting a streamline one step is always taken. To avoid
// this unneed step check to see if the termination criteria was
......@@ -2012,7 +2015,7 @@ avtStreamlineFilter::IntegrateDomain(avtStreamlineWrapper *slSeg,
if (DebugStream::Level4())
debug4<<"IntegrateDomain: slSeg->terminated= "<<slSeg->terminated<<endl;
if( ! slSeg->terminated )
if (!slSeg->terminated)
{
if (intersectObj)
slSeg->sl->SetIntersectionObject(intersectObj);
......@@ -2026,22 +2029,25 @@ avtStreamlineFilter::IntegrateDomain(avtStreamlineWrapper *slSeg,
}
else
{
if (integrationType == STREAMLINE_INTEGRATE_M3D_C1_INTEGRATOR) {
avtIVPM3DC1Field field(velocity1);
result = slSeg->sl->Advance(&field,
slSeg->terminationType,
slSeg->termination);
} else {
avtIVPVTKField field(velocity1);
result = slSeg->sl->Advance(&field,
slSeg->terminationType,
slSeg->termination);
}
if (integrationType == STREAMLINE_INTEGRATE_M3D_C1_INTEGRATOR)
{
avtIVPM3DC1Field field(velocity1);
result = slSeg->sl->Advance(&field,
slSeg->terminationType,
slSeg->termination);
}
else
{
avtIVPVTKField field(velocity1);
result = slSeg->sl->Advance(&field,
slSeg->terminationType,
slSeg->termination);
}
}
// Termination criteria was met.
slSeg->terminated = (result == avtIVPSolver::TERMINATE);
slSeg->terminated = (result == avtStreamline::TERMINATE);
if (DebugStream::Level5())
{
debug5<<"Advance:= "<<result<<endl;
......@@ -2049,60 +2055,30 @@ avtStreamlineFilter::IntegrateDomain(avtStreamlineWrapper *slSeg,
}
}
else
result = avtIVPSolver::TERMINATE;
result = avtStreamline::TERMINATE;
//slSeg->Debug();
if (result == avtIVPSolver::OUTSIDE_DOMAIN)
//Streamline exited this domain.
if (result == avtStreamline::EXIT_DOMAIN)
{
if (DebugStream::Level5())
debug5 << numSteps << " " << slSeg->sl->size() << endl;
numSteps = slSeg->sl->size() - numSteps;
slSeg->status = avtStreamlineWrapper::OUTOFBOUNDS;
DomainType oldDomain = slSeg->domain;
//Set the new domain.
SetDomain(slSeg);
int domCnt = slSeg->seedPtDomainList.size();
// We are in the same domain.
if (slSeg->seedPtDomainList.size() > 1)
//If we land in none, or the same domain, we're done.
if (domCnt == 0 || (domCnt == 1 && slSeg->domain == oldDomain))
{
if (DebugStream::Level5())
debug5 << slSeg->domain << " " << oldDomain << " " << numSteps << endl;
// pathline terminates if timestep is out of bounds.
if (doPathlines && slSeg->domain.timeStep == -1)
{
slSeg->status = avtStreamlineWrapper::TERMINATE;
}
numSteps = slSeg->sl->size() - numSteps;
if (slSeg->domain == oldDomain && numSteps == 0)
{
slSeg->status = avtStreamlineWrapper::TERMINATE;
}
else
{
slSeg->status = avtStreamlineWrapper::OUTOFBOUNDS;
}
slSeg->status = avtStreamlineWrapper::TERMINATE;
}
// Not in any or just one domain.
else
{
slSeg->status = avtStreamlineWrapper::TERMINATE;
}
slSeg->status = avtStreamlineWrapper::OUTOFBOUNDS;
}
else
{
slSeg->status = avtStreamlineWrapper::TERMINATE;
}
velocity1->Delete();
if (DebugStream::Level4())
debug4<<"::IntegrateDomain() result= "<<result<<endl;
visitTimer->StopTimer(t0, "IntegrateDomain");
return result;
}
......@@ -2160,10 +2136,7 @@ avtStreamlineFilter::IntegrateStreamline(avtStreamlineWrapper *slSeg,
double extents[6] = { 0.,0., 0.,0., 0.,0. };
intervalTree->GetElementExtents(slSeg->domain.domain, extents);
avtIVPSolver::Result result =
IntegrateDomain(slSeg, ds, extents, maxSteps);
if (DebugStream::Level5())
debug5<<"ISL: result= "<<result<<endl;
IntegrateDomain(slSeg, ds, extents, maxSteps);
//SL exited this domain.
if (slSeg->status == avtStreamlineWrapper::OUTOFBOUNDS)
......
......@@ -329,7 +329,7 @@ class AVTFILTERS_API avtStreamlineFilter :
virtual bool CheckOnDemandViability(void);
void IntegrateStreamline(avtStreamlineWrapper *slSeg, int maxSteps=-1);
avtIVPSolver::Result IntegrateDomain(avtStreamlineWrapper *slSeg,
void IntegrateDomain(avtStreamlineWrapper *slSeg,
vtkDataSet *ds,
double *extents,
int maxSteps=-1);
......
......@@ -457,6 +457,9 @@ avtIVPAdamsBashforth::ABStep(const avtIVPField* field,
// Dave Pugmire, Tue Feb 23 09:42:25 EST 2010
// Set the velStart/velEnd direction based on integration direction.
//
// Dave Pugmire, Wed May 26 13:48:24 EDT 2010
// The velStart/velEnd direction was reversed.
//
// ****************************************************************************
avtIVPSolver::Result
......@@ -558,13 +561,13 @@ avtIVPAdamsBashforth::Step(const avtIVPField* field,
if (end < 0.0)
{
ivpstep->velStart = (*field)(t,yCur);
ivpstep->velEnd = (*field)((t+h),yNew);
ivpstep->velStart = -(*field)(t,yCur);
ivpstep->velEnd = -(*field)((t+h),yNew);
}
else
{
ivpstep->velStart = - (*field)(t,yCur);
ivpstep->velEnd = - (*field)((t+h),yNew);
ivpstep->velStart = (*field)(t,yCur);
ivpstep->velEnd = (*field)((t+h),yNew);
}
// Update the history to save the last 5 vector values.
......
......@@ -49,6 +49,9 @@
#include <DebugStream.h>
#include <vtkPlane.h>
const double avtStreamline::minH = 1e-9;
// ****************************************************************************
// Method: avtStreamline constructor
//
......@@ -198,13 +201,12 @@ avtStreamline::GetVariableIdx(const std::string &var) const
//
// ****************************************************************************
avtIVPSolver::Result
avtStreamline::Result
avtStreamline::Advance(const avtIVPField* field,
avtIVPSolver::TerminateType termType,
double end)
{
avtIVPSolver::Result res = DoAdvance(_ivpSolver, field, termType, end);
return res;
return DoAdvance(_ivpSolver, field, termType, end);
}
......@@ -273,23 +275,26 @@ avtStreamline::Advance(const avtIVPField* field,
// Don't distinguish between forward/bwd integration for steps. Always add
// add them to the back of the list.
//
// Dave Pugmire, Wed May 26 13:48:24 EDT 2010
// Use a new return code.
//
// ****************************************************************************
avtIVPSolver::Result
avtStreamline::Result
avtStreamline::DoAdvance(avtIVPSolver* ivp,
const avtIVPField* field,
avtIVPSolver::TerminateType termType,
double end)
{
avtIVPSolver::Result result;
unsigned int ousideDomainCount = 0;
// catch cases where the start position is outside the
// domain of field
if (!field->IsInside(ivp->GetCurrentT(), ivp->GetCurrentY()))
return avtIVPSolver::OUTSIDE_DOMAIN;
return avtStreamline::POINT_OUTSIDE_DOMAIN;
avtStreamline::Result rc;
bool stepTaken = false;
while (1)
{
// record state for later restore, if needed
......@@ -301,27 +306,22 @@ avtStreamline::DoAdvance(avtIVPSolver* ivp,
try
{
// if (DebugStream::Level5())
// debug5<<"Step( mode= "<<termType<<" end= "<<end<<endl;
result = ivp->Step(field, termType, end, step);
// if (DebugStream::Level5())
// debug5<<" T= "<<ivp->GetCurrentT()<<" "<<ivp->GetCurrentY()<<endl;
if (intersectionsSet)
HandleIntersections(step, termType, end, &result);
HandleIntersections(step, termType, end, &result);
debug5<<"Step to "<<ivp->GetCurrentY()<<endl;
}
catch( avtIVPField::Undefined& )
catch (avtIVPField::Undefined&)
{
if (DebugStream::Level5())
debug5<< ivp->GetCurrentY() << " not in domain "
<< ivp->GetNextStepSize() << endl;
debug5<< ivp->GetCurrentY() << " not in domain "
<< ivp->GetNextStepSize() << endl;
// integrator left the domain, retry with smaller step
// if step size is below given minimum, give up
// restore old state to before failed step
double hBeforePush = ivp->GetNextStepSize();
ivp->PutState( state );
ivp->PutState(state);
if (ivp->GetNextStepSize() == 0.)
{
// This can happen if we try to look a few points out
......@@ -335,7 +335,7 @@ avtStreamline::DoAdvance(avtIVPSolver* ivp,
double h = ivp->GetNextStepSize();
h = h/2;
if (fabs(h) < 1e-9)
if (fabs(h) < minH)
{
delete step;
if (!field->HasGhostZones())
......@@ -346,22 +346,20 @@ avtStreamline::DoAdvance(avtIVPSolver* ivp,
}
if (DebugStream::Level5())
debug5<<"avtStreamline::DoAdvance() DONE result= OUTSIDE_DOMAIN "
<< "step= " << h << " count= " << ousideDomainCount << endl;
debug5<<"avtStreamline::DoAdvance() DONE result= OUTSIDE_DOMAIN "
<< "step= " << h <<endl;
return avtIVPSolver::OUTSIDE_DOMAIN;
rc = avtStreamline::EXIT_DOMAIN;
break;
}
ivp->SetNextStepSize(h);
++ousideDomainCount;
// retry step
delete step;
continue;
}
catch( std::exception& )
catch (std::exception&)
{
}
......@@ -378,24 +376,26 @@ avtStreamline::DoAdvance(avtIVPSolver* ivp,
step->ComputeScalarVariables(scalars, field);
_steps.push_back(step);
stepTaken = true;
if (result == avtIVPSolver::TERMINATE)
{
rc = avtStreamline::TERMINATE;
break;
if (ousideDomainCount && DebugStream::Level5())
debug5<<"avtStreamline::DoAdvance() DONE result= BACKIN_DOMAIN "
<< "step= " << ivp->GetNextStepSize()
<< " count= " << ousideDomainCount << endl;
ousideDomainCount = 0;
}
}
else
{
delete step;
rc = avtStreamline::ERROR;
break;
}
}
return result;
if (!stepTaken)
rc = avtStreamline::ERROR;
return rc;
}
......
......@@ -123,11 +123,22 @@ class vtkObject;
// Dave Pugmire, Tue Feb 23 09:42:25 EST 2010
// Use a vector instead of a list for the integration steps.
//
// Dave Pugmire, Wed May 26 13:48:24 EDT 2010
// New return enum.
//
// ****************************************************************************
class IVP_API avtStreamline
{
public:
enum Result
{
TERMINATE,
POINT_OUTSIDE_DOMAIN,
EXIT_DOMAIN,
ERROR,
};
enum ScalarValueType {NONE=0, SPEED=1, VORTICITY=2, SCALAR_VARIABLE=4};
typedef std::vector<avtIVPStep*>::const_iterator iterator;
......@@ -136,9 +147,9 @@ class IVP_API avtStreamline
avtStreamline();
~avtStreamline();
avtIVPSolver::Result Advance(const avtIVPField* field,
avtIVPSolver::TerminateType termType,
double end);
avtStreamline::Result Advance(const avtIVPField* field,
avtIVPSolver::TerminateType termType,
double end);
void SetScalarValueType(ScalarValueType t) {scalarValueType = t;}
void SetIntersectionObject(vtkObject *obj);
......@@ -172,10 +183,10 @@ class IVP_API avtStreamline
avtStreamline( const avtStreamline& );
avtStreamline& operator=( const avtStreamline& );
avtIVPSolver::Result DoAdvance(avtIVPSolver* ivp,
const avtIVPField* field,
avtIVPSolver::TerminateType termType,
double end);
avtStreamline::Result DoAdvance(avtIVPSolver* ivp,
const avtIVPField* field,
avtIVPSolver::TerminateType termType,
double end);
void HandleGhostZones(bool forward, double *extents);
void HandleIntersections(avtIVPStep *step,
......@@ -202,6 +213,7 @@ class IVP_API avtStreamline
// Solver.
avtIVPSolver* _ivpSolver;
static const double minH;
public:
//Bookeeping
......@@ -210,6 +222,21 @@ class IVP_API avtStreamline
std::vector<std::string> scalars;
};
inline std::ostream& operator<<( std::ostream& out, const avtStreamline::Result &res )
{
switch (res)
{
case avtStreamline::TERMINATE: out<<"TERMINATE"; break;
case avtStreamline::POINT_OUTSIDE_DOMAIN: out<<"POINTOUTSIDE_DOMAIN"; break;
case avtStreamline::EXIT_DOMAIN: out<<"EXIT_DOMAIN"; break;
case avtStreamline::ERROR: out<<"ERROR"; break;
default:
out<<"UNKNOWN_RESULT"; break;
}
return out;
}
#endif // AVT_STREAMLINE_H
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