mirror of
https://github.com/NGSolve/netgen.git
synced 2025-01-11 21:50:34 +05:00
make fieldlines code available in non-gui builds
This commit is contained in:
Matthias Hochsteger
parent
1f793e66ea
commit
9ceb2baeaa
@ -2,9 +2,9 @@ add_definitions(-DNGINTERFACE_EXPORTS)
|
||||
install(FILES soldata.hpp DESTINATION ${NG_INSTALL_DIR_INCLUDE} COMPONENT netgen_devel )
|
||||
|
||||
if(USE_GUI)
|
||||
set( LIB_VISUAL_SOURCES meshdoc.cpp mvdraw.cpp vsfieldlines.cpp vsmesh.cpp vssolution.cpp importsolution.cpp )
|
||||
set( LIB_VISUAL_SOURCES meshdoc.cpp mvdraw.cpp vsfieldlines.cpp fieldlines.cpp vsmesh.cpp vssolution.cpp importsolution.cpp )
|
||||
else(USE_GUI)
|
||||
set( LIB_VISUAL_SOURCES visual_dummy.cpp )
|
||||
set( LIB_VISUAL_SOURCES visual_dummy.cpp fieldlines.cpp)
|
||||
endif(USE_GUI)
|
||||
|
||||
add_library(visual ${NG_LIB_TYPE} ${LIB_VISUAL_SOURCES})
|
||||
@ -14,6 +14,6 @@ install( TARGETS visual ${NG_INSTALL_DIR})
|
||||
|
||||
install(FILES
|
||||
meshdoc.hpp mvdraw.hpp
|
||||
vispar.hpp visual.hpp vssolution.hpp vsfieldlines.hpp
|
||||
vispar.hpp visual.hpp vssolution.hpp fieldlines.hpp
|
||||
DESTINATION ${NG_INSTALL_DIR_INCLUDE}/visualization COMPONENT netgen_devel
|
||||
)
|
||||
|
||||
384
libsrc/visualization/fieldlines.cpp
Normal file
384
libsrc/visualization/fieldlines.cpp
Normal file
@ -0,0 +1,384 @@
|
||||
#include <mystdlib.h>
|
||||
|
||||
#include <myadt.hpp>
|
||||
#include <meshing.hpp>
|
||||
#include <csg.hpp>
|
||||
#include <stlgeom.hpp>
|
||||
|
||||
#include "fieldlines.hpp"
|
||||
|
||||
namespace netgen
|
||||
{
|
||||
RKStepper :: ~RKStepper()
|
||||
{
|
||||
delete a;
|
||||
}
|
||||
|
||||
RKStepper :: RKStepper(int type) : a(NULL), tolerance(1e100)
|
||||
{
|
||||
notrestarted = 0;
|
||||
|
||||
if (type == 0) // explicit Euler
|
||||
{
|
||||
c.SetSize(1); c[0] = 0;
|
||||
b.SetSize(1); b[0] = 1;
|
||||
steps = order = 1;
|
||||
}
|
||||
else if (type == 1) // Euler-Cauchy
|
||||
{
|
||||
c.SetSize(2); c[0] = 0; c[1] = 0.5;
|
||||
b.SetSize(2); b[0] = 0; b[1] = 1;
|
||||
NgArray<int> size(2);
|
||||
size[0] = 0; size[1] = 1;
|
||||
a = new TABLE<double>(size);
|
||||
a->Set(2,1,0.5); // Set, Get: 1-based!
|
||||
steps = order = 2;
|
||||
}
|
||||
else if (type == 2) // Simpson
|
||||
{
|
||||
c.SetSize(3); c[0] = 0; c[1] = 1; c[2] = 0.5;
|
||||
b.SetSize(3); b[0] = b[1] = 1./6.; b[2] = 2./3.;
|
||||
NgArray<int> size(3);
|
||||
size[0] = 0; size[1] = 1; size[2] = 2;
|
||||
a = new TABLE<double>(size);
|
||||
a->Set(2,1,1);
|
||||
a->Set(3,1,0.25); a->Set(3,2,0.25);
|
||||
steps = order = 3;
|
||||
}
|
||||
else if (type == 3) // classical Runge-Kutta
|
||||
{
|
||||
c.SetSize(4); c[0] = 0; c[1] = c[2] = 0.5; c[3] = 1;
|
||||
b.SetSize(4); b[0] = b[3] = 1./6.; b[1] = b[2] = 1./3.;
|
||||
NgArray<int> size(4);
|
||||
size[0] = 0; size[1] = 1; size[2] = 2; size[3] = 3;
|
||||
a = new TABLE<double>(size);
|
||||
a->Set(2,1,0.5);
|
||||
a->Set(3,1,0); a->Set(3,2,0.5);
|
||||
a->Set(4,1,0); a->Set(4,2,0); a->Set(4,3,1);
|
||||
steps = order = 4;
|
||||
}
|
||||
|
||||
K.SetSize(steps);
|
||||
}
|
||||
|
||||
void RKStepper :: StartNextValCalc(const Point<3> & astartval, const double astartt, const double ah, const bool aadaptive)
|
||||
{
|
||||
//cout << "Starting RK-Step with h=" << ah << endl;
|
||||
|
||||
stepcount = 0;
|
||||
h = ah;
|
||||
startt = astartt;
|
||||
startval = astartval;
|
||||
adaptive = aadaptive;
|
||||
adrun = 0;
|
||||
}
|
||||
|
||||
bool RKStepper :: GetNextData(Point<3> & val, double & t, double & ah)
|
||||
{
|
||||
bool finished = false;
|
||||
|
||||
if(stepcount <= steps)
|
||||
{
|
||||
t = startt + c[stepcount-1]*h;
|
||||
val = startval;
|
||||
for(int i=0; i<stepcount-1; i++)
|
||||
val += h * a->Get(stepcount,i+1) * K[i];
|
||||
}
|
||||
|
||||
|
||||
if(stepcount == steps)
|
||||
{
|
||||
val = startval;
|
||||
for(int i=0; i<steps; i++)
|
||||
val += h * b[i] * K[i];
|
||||
|
||||
if(adaptive)
|
||||
{
|
||||
if(adrun == 0)
|
||||
{
|
||||
stepcount = 0;
|
||||
h *= 0.5;
|
||||
adrun = 1;
|
||||
valh = val;
|
||||
}
|
||||
else if (adrun == 1)
|
||||
{
|
||||
stepcount = 0;
|
||||
startval_bak = startval;
|
||||
startval = val;
|
||||
startt_bak = startt;
|
||||
startt += h;//0.5*h;
|
||||
adrun = 2;
|
||||
}
|
||||
else if (adrun == 2)
|
||||
{
|
||||
Point<3> valh2 = val;
|
||||
val = valh2 + 1./(pow(2.,order)-1.) * (valh2 - valh);
|
||||
auto errvec = val - valh;
|
||||
|
||||
double err = errvec.Length();
|
||||
|
||||
double fac = 0.7 * pow(tolerance/err,1./(order+1.));
|
||||
if(fac > 1.3) fac = 1.3;
|
||||
|
||||
if(fac < 1 || notrestarted >= 2)
|
||||
ah = 2.*h * fac;
|
||||
|
||||
if(err < tolerance)
|
||||
{
|
||||
finished = true;
|
||||
notrestarted++;
|
||||
//(*testout) << "finished RK-Step, new h=" << ah << " tolerance " << tolerance << " err " << err << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
//ah *= 0.9;
|
||||
notrestarted = 0;
|
||||
//(*testout) << "restarting h " << 2.*h << " ah " << ah << " tolerance " << tolerance << " err " << err << endl;
|
||||
StartNextValCalc(startval_bak,startt_bak, ah, adaptive);
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
t = startt + h;
|
||||
finished = true;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
if(stepcount == 0)
|
||||
{
|
||||
t = startt + c[stepcount]*h;
|
||||
val = startval;
|
||||
for(int i=0; i<stepcount; i++)
|
||||
val += h * a->Get(stepcount,i) * K[i];
|
||||
}
|
||||
|
||||
return finished;
|
||||
}
|
||||
|
||||
|
||||
bool RKStepper :: FeedNextF(const Vec<3> & f)
|
||||
{
|
||||
K[stepcount] = f;
|
||||
stepcount++;
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
|
||||
void FieldLineCalc :: GenerateFieldLines(Array<Point<3>> & potential_startpoints, const int numlines)
|
||||
{
|
||||
|
||||
|
||||
Array<Point<3>> line_points;
|
||||
Array<double> line_values;
|
||||
Array<bool> drawelems;
|
||||
Array<int> dirstart;
|
||||
pstart.SetSize0();
|
||||
pend.SetSize0();
|
||||
values.SetSize0();
|
||||
|
||||
double crit = 1.0;
|
||||
|
||||
if(randomized)
|
||||
{
|
||||
double sum = 0;
|
||||
double lami[3];
|
||||
Vec<3> v;
|
||||
|
||||
for(int i=0; i<potential_startpoints.Size(); i++)
|
||||
{
|
||||
int elnr = mesh.GetElementOfPoint(potential_startpoints[i],lami,true) - 1;
|
||||
if(elnr == -1)
|
||||
continue;
|
||||
|
||||
mesh.SetPointSearchStartElement(elnr);
|
||||
|
||||
func(elnr, lami, v);
|
||||
sum += v.Length();
|
||||
}
|
||||
|
||||
crit = sum/double(numlines);
|
||||
}
|
||||
|
||||
|
||||
int calculated = 0;
|
||||
|
||||
cout << endl;
|
||||
|
||||
|
||||
for(int i=0; i<potential_startpoints.Size(); i++)
|
||||
{
|
||||
cout << "\rFieldline Calculation " << int(100.*i/potential_startpoints.Size()) << "%"; cout.flush();
|
||||
|
||||
if(randomized)
|
||||
SetCriticalValue((double(rand())/RAND_MAX)*crit);
|
||||
|
||||
if(calculated >= numlines) break;
|
||||
|
||||
Calc(potential_startpoints[i],line_points,line_values,drawelems,dirstart);
|
||||
|
||||
bool usable = false;
|
||||
|
||||
for(int j=1; j<dirstart.Size(); j++)
|
||||
for(int k=dirstart[j-1]; k<dirstart[j]-1; k++)
|
||||
{
|
||||
if(!drawelems[k] || !drawelems[k+1]) continue;
|
||||
|
||||
usable = true;
|
||||
pstart.Append(line_points[k]);
|
||||
pend.Append(line_points[k+1]);
|
||||
values.Append( 0.5*(line_values[k]+line_values[k+1]) );
|
||||
}
|
||||
|
||||
if(usable) calculated++;
|
||||
}
|
||||
cout << "\rFieldline Calculation " << 100 << "%" << endl;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
FieldLineCalc :: FieldLineCalc(const Mesh & amesh, const VectorFunction & afunc,
|
||||
const double rel_length, const int amaxpoints,
|
||||
const double rel_thickness, const double rel_tolerance, const int rk_type, const int adirection) :
|
||||
mesh(amesh), func(afunc), stepper(rk_type)
|
||||
{
|
||||
mesh.GetBox (pmin, pmax);
|
||||
rad = 0.5 * Dist (pmin, pmax);
|
||||
|
||||
|
||||
maxlength = (rel_length > 0) ? rel_length : 0.5;
|
||||
maxlength *= 2.*rad;
|
||||
|
||||
thickness = (rel_thickness > 0) ? rel_thickness : 0.0015;
|
||||
thickness *= 2.*rad;
|
||||
|
||||
double auxtolerance = (rel_tolerance > 0) ? rel_tolerance : 1.5e-3;
|
||||
auxtolerance *= 2.*rad;
|
||||
|
||||
stepper.SetTolerance(auxtolerance);
|
||||
|
||||
direction = adirection;
|
||||
|
||||
|
||||
maxpoints = amaxpoints;
|
||||
|
||||
if(direction == 0)
|
||||
{
|
||||
maxlength *= 0.5;
|
||||
maxpoints /= 2;
|
||||
}
|
||||
|
||||
|
||||
critical_value = -1;
|
||||
|
||||
randomized = false;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
void FieldLineCalc :: Calc(const Point<3> & startpoint, Array<Point<3>> & points, Array<double> & vals, Array<bool> & drawelems, Array<int> & dirstart)
|
||||
{
|
||||
Vec<3> v = 0.0;
|
||||
double startlami[3] = {0.0, 0.0, 0.0};
|
||||
|
||||
points.SetSize(0);
|
||||
vals.SetSize(0);
|
||||
drawelems.SetSize(0);
|
||||
|
||||
dirstart.SetSize(0);
|
||||
dirstart.Append(0);
|
||||
|
||||
|
||||
int startelnr = mesh.GetElementOfPoint(startpoint,startlami,true) - 1;
|
||||
(*testout) << "p = " << startpoint << "; elnr = " << startelnr << endl;
|
||||
if (startelnr == -1)
|
||||
return;
|
||||
|
||||
mesh.SetPointSearchStartElement(startelnr);
|
||||
|
||||
Vec<3> startv;
|
||||
bool startdraw = func(startelnr, startlami, startv);
|
||||
|
||||
double startval = startv.Length();
|
||||
|
||||
if(critical_value > 0 && fabs(startval) < critical_value)
|
||||
return;
|
||||
|
||||
//cout << "p = " << startpoint << "; elnr = " << startelnr << endl;
|
||||
|
||||
|
||||
|
||||
for(int dir = 1; dir >= -1; dir -= 2)
|
||||
{
|
||||
if(dir*direction < 0) continue;
|
||||
|
||||
points.Append(startpoint);
|
||||
vals.Append(startval);
|
||||
drawelems.Append(startdraw);
|
||||
|
||||
double h = 0.001*rad/startval; // otherwise no nice lines; should be made accessible from outside
|
||||
|
||||
v = startv;
|
||||
if(dir == -1) v *= -1.;
|
||||
|
||||
int elnr = startelnr;
|
||||
double lami[3] = { startlami[0], startlami[1], startlami[2]};
|
||||
|
||||
|
||||
for(double length = 0; length < maxlength; length += h*vals.Last())
|
||||
{
|
||||
if(v.Length() < 1e-12*rad)
|
||||
{
|
||||
(*testout) << "Current fieldlinecalculation came to a stillstand at " << points.Last() << endl;
|
||||
break;
|
||||
}
|
||||
|
||||
double dummyt;
|
||||
stepper.StartNextValCalc(points.Last(),dummyt,h,true);
|
||||
stepper.FeedNextF(v);
|
||||
bool drawelem = false;
|
||||
|
||||
Point<3> newp;
|
||||
while(!stepper.GetNextData(newp,dummyt,h) && elnr != -1)
|
||||
{
|
||||
elnr = mesh.GetElementOfPoint(newp,lami,true) - 1;
|
||||
if(elnr != -1)
|
||||
{
|
||||
mesh.SetPointSearchStartElement(elnr);
|
||||
drawelem = func(elnr, lami, v);
|
||||
if(dir == -1) v *= -1.;
|
||||
stepper.FeedNextF(v);
|
||||
}
|
||||
}
|
||||
|
||||
if (elnr == -1)
|
||||
{
|
||||
//cout << "direction " <<dir << " reached the wall." << endl;
|
||||
break;
|
||||
}
|
||||
|
||||
points.Append(newp);
|
||||
vals.Append(v.Length());
|
||||
drawelems.Append(drawelem);
|
||||
|
||||
if(points.Size() % 40 == 0 && points.Size() > 1)
|
||||
(*testout) << "Points in current fieldline: " << points.Size() << ", current position: " << newp << endl;
|
||||
|
||||
if(maxpoints > 0 && points.Size() >= maxpoints)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
//cout << "length " << length << " h " << h << " vals.Last() " << vals.Last() << " maxlength " << maxlength << endl;
|
||||
}
|
||||
dirstart.Append(points.Size());
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
@ -17,384 +17,6 @@ namespace netgen
|
||||
|
||||
// extern shared_ptr<Mesh> mesh;
|
||||
|
||||
|
||||
|
||||
RKStepper :: ~RKStepper()
|
||||
{
|
||||
delete a;
|
||||
}
|
||||
|
||||
RKStepper :: RKStepper(int type) : a(NULL), tolerance(1e100)
|
||||
{
|
||||
notrestarted = 0;
|
||||
|
||||
if (type == 0) // explicit Euler
|
||||
{
|
||||
c.SetSize(1); c[0] = 0;
|
||||
b.SetSize(1); b[0] = 1;
|
||||
steps = order = 1;
|
||||
}
|
||||
else if (type == 1) // Euler-Cauchy
|
||||
{
|
||||
c.SetSize(2); c[0] = 0; c[1] = 0.5;
|
||||
b.SetSize(2); b[0] = 0; b[1] = 1;
|
||||
NgArray<int> size(2);
|
||||
size[0] = 0; size[1] = 1;
|
||||
a = new TABLE<double>(size);
|
||||
a->Set(2,1,0.5); // Set, Get: 1-based!
|
||||
steps = order = 2;
|
||||
}
|
||||
else if (type == 2) // Simpson
|
||||
{
|
||||
c.SetSize(3); c[0] = 0; c[1] = 1; c[2] = 0.5;
|
||||
b.SetSize(3); b[0] = b[1] = 1./6.; b[2] = 2./3.;
|
||||
NgArray<int> size(3);
|
||||
size[0] = 0; size[1] = 1; size[2] = 2;
|
||||
a = new TABLE<double>(size);
|
||||
a->Set(2,1,1);
|
||||
a->Set(3,1,0.25); a->Set(3,2,0.25);
|
||||
steps = order = 3;
|
||||
}
|
||||
else if (type == 3) // classical Runge-Kutta
|
||||
{
|
||||
c.SetSize(4); c[0] = 0; c[1] = c[2] = 0.5; c[3] = 1;
|
||||
b.SetSize(4); b[0] = b[3] = 1./6.; b[1] = b[2] = 1./3.;
|
||||
NgArray<int> size(4);
|
||||
size[0] = 0; size[1] = 1; size[2] = 2; size[3] = 3;
|
||||
a = new TABLE<double>(size);
|
||||
a->Set(2,1,0.5);
|
||||
a->Set(3,1,0); a->Set(3,2,0.5);
|
||||
a->Set(4,1,0); a->Set(4,2,0); a->Set(4,3,1);
|
||||
steps = order = 4;
|
||||
}
|
||||
|
||||
K.SetSize(steps);
|
||||
}
|
||||
|
||||
void RKStepper :: StartNextValCalc(const Point<3> & astartval, const double astartt, const double ah, const bool aadaptive)
|
||||
{
|
||||
//cout << "Starting RK-Step with h=" << ah << endl;
|
||||
|
||||
stepcount = 0;
|
||||
h = ah;
|
||||
startt = astartt;
|
||||
startval = astartval;
|
||||
adaptive = aadaptive;
|
||||
adrun = 0;
|
||||
}
|
||||
|
||||
bool RKStepper :: GetNextData(Point<3> & val, double & t, double & ah)
|
||||
{
|
||||
bool finished = false;
|
||||
|
||||
if(stepcount <= steps)
|
||||
{
|
||||
t = startt + c[stepcount-1]*h;
|
||||
val = startval;
|
||||
for(int i=0; i<stepcount-1; i++)
|
||||
val += h * a->Get(stepcount,i+1) * K[i];
|
||||
}
|
||||
|
||||
|
||||
if(stepcount == steps)
|
||||
{
|
||||
val = startval;
|
||||
for(int i=0; i<steps; i++)
|
||||
val += h * b[i] * K[i];
|
||||
|
||||
if(adaptive)
|
||||
{
|
||||
if(adrun == 0)
|
||||
{
|
||||
stepcount = 0;
|
||||
h *= 0.5;
|
||||
adrun = 1;
|
||||
valh = val;
|
||||
}
|
||||
else if (adrun == 1)
|
||||
{
|
||||
stepcount = 0;
|
||||
startval_bak = startval;
|
||||
startval = val;
|
||||
startt_bak = startt;
|
||||
startt += h;//0.5*h;
|
||||
adrun = 2;
|
||||
}
|
||||
else if (adrun == 2)
|
||||
{
|
||||
Point<3> valh2 = val;
|
||||
val = valh2 + 1./(pow(2.,order)-1.) * (valh2 - valh);
|
||||
auto errvec = val - valh;
|
||||
|
||||
double err = errvec.Length();
|
||||
|
||||
double fac = 0.7 * pow(tolerance/err,1./(order+1.));
|
||||
if(fac > 1.3) fac = 1.3;
|
||||
|
||||
if(fac < 1 || notrestarted >= 2)
|
||||
ah = 2.*h * fac;
|
||||
|
||||
if(err < tolerance)
|
||||
{
|
||||
finished = true;
|
||||
notrestarted++;
|
||||
//(*testout) << "finished RK-Step, new h=" << ah << " tolerance " << tolerance << " err " << err << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
//ah *= 0.9;
|
||||
notrestarted = 0;
|
||||
//(*testout) << "restarting h " << 2.*h << " ah " << ah << " tolerance " << tolerance << " err " << err << endl;
|
||||
StartNextValCalc(startval_bak,startt_bak, ah, adaptive);
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
t = startt + h;
|
||||
finished = true;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
if(stepcount == 0)
|
||||
{
|
||||
t = startt + c[stepcount]*h;
|
||||
val = startval;
|
||||
for(int i=0; i<stepcount; i++)
|
||||
val += h * a->Get(stepcount,i) * K[i];
|
||||
}
|
||||
|
||||
return finished;
|
||||
}
|
||||
|
||||
|
||||
bool RKStepper :: FeedNextF(const Vec<3> & f)
|
||||
{
|
||||
K[stepcount] = f;
|
||||
stepcount++;
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
|
||||
void FieldLineCalc :: GenerateFieldLines(Array<Point<3>> & potential_startpoints, const int numlines)
|
||||
{
|
||||
|
||||
|
||||
Array<Point<3>> line_points;
|
||||
Array<double> line_values;
|
||||
Array<bool> drawelems;
|
||||
Array<int> dirstart;
|
||||
pstart.SetSize0();
|
||||
pend.SetSize0();
|
||||
values.SetSize0();
|
||||
|
||||
double crit = 1.0;
|
||||
|
||||
if(randomized)
|
||||
{
|
||||
double sum = 0;
|
||||
double lami[3];
|
||||
Vec<3> v;
|
||||
|
||||
for(int i=0; i<potential_startpoints.Size(); i++)
|
||||
{
|
||||
int elnr = mesh.GetElementOfPoint(potential_startpoints[i],lami,true) - 1;
|
||||
if(elnr == -1)
|
||||
continue;
|
||||
|
||||
mesh.SetPointSearchStartElement(elnr);
|
||||
|
||||
func(elnr, lami, v);
|
||||
sum += v.Length();
|
||||
}
|
||||
|
||||
crit = sum/double(numlines);
|
||||
}
|
||||
|
||||
|
||||
int calculated = 0;
|
||||
|
||||
cout << endl;
|
||||
|
||||
|
||||
for(int i=0; i<potential_startpoints.Size(); i++)
|
||||
{
|
||||
cout << "\rFieldline Calculation " << int(100.*i/potential_startpoints.Size()) << "%"; cout.flush();
|
||||
|
||||
if(randomized)
|
||||
SetCriticalValue((double(rand())/RAND_MAX)*crit);
|
||||
|
||||
if(calculated >= numlines) break;
|
||||
|
||||
Calc(potential_startpoints[i],line_points,line_values,drawelems,dirstart);
|
||||
|
||||
bool usable = false;
|
||||
|
||||
for(int j=1; j<dirstart.Size(); j++)
|
||||
for(int k=dirstart[j-1]; k<dirstart[j]-1; k++)
|
||||
{
|
||||
if(!drawelems[k] || !drawelems[k+1]) continue;
|
||||
|
||||
usable = true;
|
||||
pstart.Append(line_points[k]);
|
||||
pend.Append(line_points[k+1]);
|
||||
values.Append( 0.5*(line_values[k]+line_values[k+1]) );
|
||||
}
|
||||
|
||||
if(usable) calculated++;
|
||||
}
|
||||
cout << "\rFieldline Calculation " << 100 << "%" << endl;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
FieldLineCalc :: FieldLineCalc(const Mesh & amesh, const VectorFunction & afunc,
|
||||
const double rel_length, const int amaxpoints,
|
||||
const double rel_thickness, const double rel_tolerance, const int rk_type, const int adirection) :
|
||||
mesh(amesh), func(afunc), stepper(rk_type)
|
||||
{
|
||||
mesh.GetBox (pmin, pmax);
|
||||
rad = 0.5 * Dist (pmin, pmax);
|
||||
|
||||
|
||||
maxlength = (rel_length > 0) ? rel_length : 0.5;
|
||||
maxlength *= 2.*rad;
|
||||
|
||||
thickness = (rel_thickness > 0) ? rel_thickness : 0.0015;
|
||||
thickness *= 2.*rad;
|
||||
|
||||
double auxtolerance = (rel_tolerance > 0) ? rel_tolerance : 1.5e-3;
|
||||
auxtolerance *= 2.*rad;
|
||||
|
||||
stepper.SetTolerance(auxtolerance);
|
||||
|
||||
direction = adirection;
|
||||
|
||||
|
||||
maxpoints = amaxpoints;
|
||||
|
||||
if(direction == 0)
|
||||
{
|
||||
maxlength *= 0.5;
|
||||
maxpoints /= 2;
|
||||
}
|
||||
|
||||
|
||||
critical_value = -1;
|
||||
|
||||
randomized = false;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
void FieldLineCalc :: Calc(const Point<3> & startpoint, Array<Point<3>> & points, Array<double> & vals, Array<bool> & drawelems, Array<int> & dirstart)
|
||||
{
|
||||
Vec<3> v = 0.0;
|
||||
double startlami[3] = {0.0, 0.0, 0.0};
|
||||
|
||||
points.SetSize(0);
|
||||
vals.SetSize(0);
|
||||
drawelems.SetSize(0);
|
||||
|
||||
dirstart.SetSize(0);
|
||||
dirstart.Append(0);
|
||||
|
||||
|
||||
int startelnr = mesh.GetElementOfPoint(startpoint,startlami,true) - 1;
|
||||
(*testout) << "p = " << startpoint << "; elnr = " << startelnr << endl;
|
||||
if (startelnr == -1)
|
||||
return;
|
||||
|
||||
mesh.SetPointSearchStartElement(startelnr);
|
||||
|
||||
Vec<3> startv;
|
||||
bool startdraw = func(startelnr, startlami, startv);
|
||||
|
||||
double startval = startv.Length();
|
||||
|
||||
if(critical_value > 0 && fabs(startval) < critical_value)
|
||||
return;
|
||||
|
||||
//cout << "p = " << startpoint << "; elnr = " << startelnr << endl;
|
||||
|
||||
|
||||
|
||||
for(int dir = 1; dir >= -1; dir -= 2)
|
||||
{
|
||||
if(dir*direction < 0) continue;
|
||||
|
||||
points.Append(startpoint);
|
||||
vals.Append(startval);
|
||||
drawelems.Append(startdraw);
|
||||
|
||||
double h = 0.001*rad/startval; // otherwise no nice lines; should be made accessible from outside
|
||||
|
||||
v = startv;
|
||||
if(dir == -1) v *= -1.;
|
||||
|
||||
int elnr = startelnr;
|
||||
double lami[3] = { startlami[0], startlami[1], startlami[2]};
|
||||
|
||||
|
||||
for(double length = 0; length < maxlength; length += h*vals.Last())
|
||||
{
|
||||
if(v.Length() < 1e-12*rad)
|
||||
{
|
||||
(*testout) << "Current fieldlinecalculation came to a stillstand at " << points.Last() << endl;
|
||||
break;
|
||||
}
|
||||
|
||||
double dummyt;
|
||||
stepper.StartNextValCalc(points.Last(),dummyt,h,true);
|
||||
stepper.FeedNextF(v);
|
||||
bool drawelem = false;
|
||||
|
||||
Point<3> newp;
|
||||
while(!stepper.GetNextData(newp,dummyt,h) && elnr != -1)
|
||||
{
|
||||
elnr = mesh.GetElementOfPoint(newp,lami,true) - 1;
|
||||
if(elnr != -1)
|
||||
{
|
||||
mesh.SetPointSearchStartElement(elnr);
|
||||
drawelem = func(elnr, lami, v);
|
||||
if(dir == -1) v *= -1.;
|
||||
stepper.FeedNextF(v);
|
||||
}
|
||||
}
|
||||
|
||||
if (elnr == -1)
|
||||
{
|
||||
//cout << "direction " <<dir << " reached the wall." << endl;
|
||||
break;
|
||||
}
|
||||
|
||||
points.Append(newp);
|
||||
vals.Append(v.Length());
|
||||
drawelems.Append(drawelem);
|
||||
|
||||
if(points.Size() % 40 == 0 && points.Size() > 1)
|
||||
(*testout) << "Points in current fieldline: " << points.Size() << ", current position: " << newp << endl;
|
||||
|
||||
if(maxpoints > 0 && points.Size() >= maxpoints)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
//cout << "length " << length << " h " << h << " vals.Last() " << vals.Last() << " maxlength " << maxlength << endl;
|
||||
}
|
||||
dirstart.Append(points.Size());
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void VisualSceneSolution :: BuildFieldLinesFromBox(Array<Point<3>> & startpoints)
|
||||
{
|
||||
shared_ptr<Mesh> mesh = GetMesh();
|
||||
|
||||
Loading…
Reference in New Issue
Block a user