ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-24 04:50:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

extrusion fixes, reduce warnings

Browse Source
This commit is contained in:
Joachim Schoeberl 2009-04-19 21:15:26 +00:00
parent ff3eaf1119
commit 2584903baa
26 changed files with 4192 additions and 4166 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
libsrc/csg/csgeom.cpp
View File

@ -53,14 +53,14 @@ namespace netgen
CSGeometry :: CSGeometry ()
: boundingbox (default_boundingbox),
identicsurfaces (100), filename(""), ideps(1e-9)
identicsurfaces (100), ideps(1e-9), filename("")
{
;
}
CSGeometry :: CSGeometry (const string & afilename)
: boundingbox (default_boundingbox),
identicsurfaces (100), filename(afilename), ideps(1e-9)
identicsurfaces (100), ideps(1e-9), filename(afilename)
{
changeval++;
}

2
libsrc/csg/csgeom.hpp
View File

@ -134,11 +134,9 @@ private:
double ideps;
/// filename of inputfile
string filename;
public:
CSGeometry ();
CSGeometry (const string & afilename);

113
libsrc/csg/extrusion.cpp
View File

@ -356,25 +356,27 @@ namespace netgen
int seg;
CalcProj (point, p2d, seg, t_path);
Point<3> phi;
Vec<3> phip, phipp, phi_minus_point;
path->GetSpline(seg).GetDerivatives(t_path, phi, phip, phipp);
phi_minus_point = phi-point;
Vec<3> grad_t = (1.0/(phipp*phi_minus_point + phip*phip)) * phip;
Vec<3> grad_xbar, grad_ybar;
/*
Vec<3> dphi_dX[3];
Vec<3> dy_dir_dX[3];
Vec<3> dx_dir_dX[3];
for(int i = 0; i < 3; i++)
dphi_dX[i] = grad_t(i)*phip;
Vec<3> dx_dir_dX[3];
Vec<3> dy_dir_dX[3];
Vec<3> dz_dir_dX[3];
double lphip = phip.Length();
Vec<3> dy_dir_dt = (1./lphip) * phipp - ((phip*phipp)/pow(lphip,3)) * phip;
Vec<3> dx_dir_dt = Cross(dy_dir_dt,z_dir[seg]);
for(int i = 0; i < 3; i++)
@ -383,22 +385,26 @@ namespace netgen
dx_dir_dX[i] = grad_t(i) * dx_dir_dt;
Vec<3> grad_xbar;
for(int i=0; i<3; i++)
grad_xbar(i) = -1.*(phi_minus_point * dx_dir_dX[i]) + x_dir[seg](i) - x_dir[seg] * dphi_dX[i];
double zy = z_dir[seg]*y_dir[seg];
Vec<3> grad_ybar;
Vec<3> aux = z_dir[seg] - zy*y_dir[seg];
for(int i=0; i<3; i++)
grad_ybar(i) = ( (z_dir[seg]*dy_dir_dX[i])*y_dir[seg] + zy*dy_dir_dX[i] ) * phi_minus_point +
aux[i] -
aux * dphi_dX[i];
*/
// new version by JS
Vec<3> hex, hey, hez, dex, dey, dez;
CalcLocalCoordinatesDeriv (seg, t_path, hex, hey, hez, dex, dey, dez);
grad_xbar = hex - (phi_minus_point*dex + hex*phip) * grad_t;
grad_ybar = hez - (phi_minus_point*dez + hez*phip) * grad_t;
double dFdxbar = 2.*profile_spline_coeff(0)*p2d(0) +
profile_spline_coeff(2)*p2d(1) + profile_spline_coeff(3);
@ -409,31 +415,8 @@ namespace netgen
grad = dFdxbar * grad_xbar + dFdybar * grad_ybar;
/*
{
cout << "grad_xbar = " << grad_xbar << endl;
cout << "grad_ybar = " << grad_ybar << endl;
Vec<3> numgradx, numgrady;
for (int i = 0; i < 3; i++)
{
Point<3> hpl = point, hpr = point;
hpl(i) -= 1e-6; hpr(i) += 1e-6;
Point<2> p2dl, p2dr;
CalcProj (hpl, p2dl, seg, t_path);
CalcProj (hpr, p2dr, seg, t_path);
numgradx(i) = (p2dr(0)-p2dl(0)) / (2e-6);
numgrady(i) = (p2dr(1)-p2dl(1)) / (2e-6);
}
cout << "num grad_xbar" << numgradx << endl;
cout << "num grad_ybar" << numgrady << endl;
}
{
cout << "grad = " << grad << " =?= ";
Vec<3> numgrad;
@ -449,24 +432,6 @@ namespace netgen
}
cout << " numgrad = " << numgrad << endl;
for (int i = 0; i < 3; i++)
{
Point<3> hpl = point;
Point<3> hpr = point;
hpl(i) -= 1e-4;
hpr(i) += 1e-4;
double vall = CalcFunctionValue (hpl);
double valr = CalcFunctionValue (hpr);
numgrad(i) = (valr - vall) / (2e-4);
}
cout << " numgrad2 = " << numgrad << endl;
}
static int cnt = 0;
cnt++;
if (cnt == 10000)
exit (0);
*/
}
@ -778,6 +743,52 @@ namespace netgen
}
void ExtrusionFace ::
CalcLocalCoordinates (int seg, double t,
Vec<3> & ex, Vec<3> & ey, Vec<3> & ez) const
{
ey = path->GetSpline(seg).GetTangent(t);
ey /= ey.Length();
ex = Cross (ey, glob_z_direction);
ex /= ex.Length();
ez = Cross (ex, ey);
}
void ExtrusionFace ::
CalcLocalCoordinatesDeriv (int seg, double t,
Vec<3> & ex, Vec<3> & ey, Vec<3> & ez,
Vec<3> & dex, Vec<3> & dey, Vec<3> & dez) const
{
Point<3> point;
Vec<3> first, second;
path->GetSpline(seg).GetDerivatives (t, point, first, second);
ey = first;
ex = Cross (ey, glob_z_direction);
ez = Cross (ex, ey);
dey = second;
dex = Cross (dey, glob_z_direction);
dez = Cross (dex, ey) + Cross (ex, dey);
double lenx = ex.Length();
double leny = ey.Length();
double lenz = ez.Length();
ex /= lenx;
ey /= leny;
ez /= lenz;
dex /= lenx;
dex -= (dex * ex) * ex;
dey /= leny;
dey -= (dey * ey) * ey;
dez /= lenz;
dez -= (dez * ez) * ez;
}
Extrusion :: Extrusion(const SplineGeometry<3> & path_in,
const SplineGeometry<2> & profile_in,

8
libsrc/csg/extrusion.hpp
View File

@ -92,6 +92,14 @@ public:
double GetProfilePar(void) const {return profile_par;}
void GetRawData(Array<double> & data) const;
void CalcLocalCoordinates (int seg, double t,
Vec<3> & ex, Vec<3> & ey, Vec<3> & ez) const;
void CalcLocalCoordinatesDeriv (int seg, double t,
Vec<3> & ex, Vec<3> & ey, Vec<3> & ez,
Vec<3> & dex, Vec<3> & dey, Vec<3> & dez) const;
};

34
libsrc/csg/identify.cpp
View File

@ -256,7 +256,7 @@ Identifyable (const Point<3> & p1, const Point<3> & p2) const
int PeriodicIdentification ::
GetIdentifiedPoint (class Mesh & mesh, int pi)
{
const Surface * sold, *snew;
const Surface *snew;
const Point<3> & p = mesh.Point (pi);
if (s1->PointOnSurface (p))
@ -280,9 +280,8 @@ GetIdentifiedPoint (class Mesh & mesh, int pi)
Point<3> hp = p;
snew->Project (hp);
int i;
int newpi = 0;
for (i = 1; i <= mesh.GetNP(); i++)
for (int i = 1; i <= mesh.GetNP(); i++)
if (Dist2 (mesh.Point(i), hp) < 1e-12)
{
newpi = i;
@ -442,7 +441,7 @@ BuildSurfaceElements (Array<Segment> & segs,
Mesh & mesh, const Surface * surf)
{
int found = 0;
int fother;
int fother = -1;
int facei = segs.Get(1).si;
int surfnr = mesh.GetFaceDescriptor(facei).SurfNr();
@ -876,7 +875,7 @@ ShortEdge (const SpecialPoint & sp1, const SpecialPoint & sp2) const
int CloseSurfaceIdentification ::
GetIdentifiedPoint (class Mesh & mesh, int pi)
{
const Surface * sold, *snew;
const Surface *snew;
const Point<3> & p = mesh.Point (pi);
Array<int,PointIndex::BASE> identmap(mesh.GetNP());
@ -1062,7 +1061,7 @@ void CloseSurfaceIdentification :: IdentifyPoints (Mesh & mesh)
void CloseSurfaceIdentification :: IdentifyFaces (class Mesh & mesh)
{
int fi1, fi2, side;
int s1rep, s2rep;
int s1rep = -1, s2rep = -1;
for (int i = 0; i < geom.GetNSurf(); i++)
{
@ -1305,9 +1304,9 @@ BuildSurfaceElements2 (Array<Segment> & segs,
if (!segs.Size()) return;
bool found = 0;
int fother = -1;
int fother;
int facei = segs.Get(1).si;
int facei = segs[0].si;
int surfnr = mesh.GetFaceDescriptor(facei).SurfNr();
@ -1356,11 +1355,7 @@ BuildSurfaceElements2 (Array<Segment> & segs,
Element2d newel(sel.GetType());
newel.SetIndex (facei);
for (int k = 1; k <= sel.GetNP(); k++)
{
newel.PNum(k) =
GetIdentifiedPoint (mesh, sel.PNum(k));
// cout << "id-point = " << sel.PNum(k) << ", np = " << newel.PNum(k) << endl;
}
newel.PNum(k) = GetIdentifiedPoint (mesh, sel.PNum(k));
Vec<3> nt = Cross (Point<3> (mesh.Point (newel.PNum(2)))-
Point<3> (mesh.Point (newel.PNum(1))),
@ -1572,12 +1567,9 @@ Identifyable (const SpecialPoint & sp1, const SpecialPoint & sp2,
void CloseEdgesIdentification :: IdentifyPoints (Mesh & mesh)
{
int i, j;
int i1, i2;
int np = mesh.GetNP();
for (i1 = 1; i1 <= np; i1++)
for (i2 = 1; i2 <= np; i2++)
for (int i1 = 1; i1 <= np; i1++)
for (int i2 = 1; i2 <= np; i2++)
{
if (i2 == i1)
continue;
@ -1617,15 +1609,13 @@ void CloseEdgesIdentification ::
BuildSurfaceElements (Array<Segment> & segs,
Mesh & mesh, const Surface * surf)
{
int i1, i2;
int found = 0;
int i, j, k;
if (surf != facet)
return;
for (i1 = 1; i1 <= segs.Size(); i1++)
for (i2 = 1; i2 < i1; i2++)
for (int i1 = 1; i1 <= segs.Size(); i1++)
for (int i2 = 1; i2 < i1; i2++)
{
const Segment & s1 = segs.Get(i1);
const Segment & s2 = segs.Get(i2);

2
libsrc/csg/revolution.cpp
View File

@ -45,7 +45,7 @@ namespace netgen
bool first,
bool last,
const int id_in) :
spline(&spline_in), p0(p), v_axis(vec), isfirst(first), islast(last), id(id_in)
isfirst(first), islast(last), spline(&spline_in), p0(p), v_axis(vec), id(id_in)
{
deletable = false;
Init();

14
libsrc/csg/singularref.cpp
View File

@ -9,11 +9,11 @@ namespace netgen
{
SingularEdge :: SingularEdge (double abeta, int adomnr,
const CSGeometry & ageom,
const Solid * asol1,
const Solid * asol2, double sf,
const double maxh_at_initialization)
: geom(ageom), domnr(adomnr)
const CSGeometry & ageom,
const Solid * asol1,
const Solid * asol2, double sf,
const double maxh_at_initialization)
: domnr(adomnr), geom(ageom)
{
beta = abeta;
maxhinit = maxh_at_initialization;
@ -86,8 +86,8 @@ void SingularEdge :: FindPointsOnEdge (class Mesh & mesh)
int surfi1 = geom.GetSurfaceClassRepresentant(mesh[si].surfnr1);
int surfi2 = geom.GetSurfaceClassRepresentant(mesh[si].surfnr2);
if (si1.Contains(surfi1) && si2.Contains(surfi2) ||
si1.Contains(surfi2) && si2.Contains(surfi1))
if ( (si1.Contains(surfi1) && si2.Contains(surfi2)) ||
(si1.Contains(surfi2) && si2.Contains(surfi1)) )
// if (onedge)
{

2
libsrc/csg/specpoin.cpp
View File

@ -1368,7 +1368,7 @@ namespace netgen
Array<int> surfind, rep_surfind, surfind2, rep_surfind2, surfind3;
Array<Vec<3> > normalvecs;
Vec<3> nsurf;
Vec<3> nsurf = 0.0;
Array<int> specpoint2point;
specpoints.SetSize (0);

9
libsrc/geom2d/genmesh2d.cpp
View File

@ -80,14 +80,11 @@ namespace netgen
mesh->SetNBCNames(maxsegmentindex);
for ( int sindex = 0; sindex < maxsegmentindex; sindex++ )
{
mesh->SetBCName ( sindex, geometry.GetBCName( sindex+1 ) );
}
mesh->SetBCName ( sindex, geometry.GetBCName( sindex+1 ) );
for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
{
(*mesh)[si].SetBCName ( (*mesh).GetBCNamePtr( (*mesh)[si].si-1 ) );
}
(*mesh)[si].SetBCName ( (*mesh).GetBCNamePtr( (*mesh)[si].si-1 ) );
Point3d pmin(bbox.PMin()(0), bbox.PMin()(1), -bbox.Diam());
Point3d pmax(bbox.PMax()(0), bbox.PMax()(1), bbox.Diam());

8
libsrc/geom2d/splinegeometry.cpp
View File

@ -1226,12 +1226,10 @@ void SplineGeometry<D> :: AppendDiscretePointsSegment (const Array< Point<D> > &
template<int D>
string SplineGeometry<D> :: GetBCName( const int bcnr ) const
{
string bcname;
if ( bcnames.Size() >= bcnr)
if ( bcnames[bcnr-1] )
bcname = *bcnames[bcnr-1];
else bcname = "default";
return bcname;
if (bcnames[bcnr-1] )
return *bcnames[bcnr-1];
return "default";
}
template<int D>

2
libsrc/gprim/geomobjects.hpp
View File

@ -22,7 +22,7 @@ protected:
public:
Point () { ; }
Point (double ax) { x[0] = ax; }
Point (double ax) { for (int i = 0; i < D; i++) x[i] = ax; }
Point (double ax, double ay) { x[0] = ax; x[1] = ay; }
Point (double ax, double ay, double az)
{ x[0] = ax; x[1] = ay; x[2] = az; }

2
libsrc/gprim/geomtest3d.cpp
View File

@ -429,7 +429,7 @@ int IntersectTetTriangleRef (const Point<3> ** tri, const int * tripi)
{
int i, j;
double eps = 1e-8;
double eps2 = eps * eps;
// double eps2 = eps * eps;
static Point<3> rtetp1(0, 0, 0);
static Point<3> rtetp2(1, 0, 0);

2
libsrc/interface/readtetmesh.cpp
View File

@ -50,7 +50,7 @@ namespace netgen
int numObj3D,numObj2D,numObj1D,numObj0D;
bool nullstarted;
Array<int> eldom;
int minId3D,minId2D;
int minId3D = -1, minId2D = -1;
int maxId3D(-1), maxId2D(-1), maxId1D(-1), maxId0D(-1);
Array<Array<int> *> segmentdata;
Array<Element2d* > tris;

33
libsrc/interface/readuser.cpp
View File

@ -22,8 +22,6 @@ namespace netgen
const char * filename = hfilename.c_str();
int i, j;
char reco[100];
int np, nbe;
@ -41,7 +39,7 @@ namespace netgen
in >> reco;
in >> np;
for (i = 1; i <= np; i++)
for (int i = 1; i <= np; i++)
{
Point3d p;
in >> p.X() >> p.Y() >> p.Z();
@ -53,15 +51,12 @@ namespace netgen
in >> nbe;
// int invert = globflags.GetDefineFlag ("invertsurfacemesh");
for (i = 1; i <= nbe; i++)
for (int i = 1; i <= nbe; i++)
{
Element2d el;
int hi;
el.SetIndex(1);
// in >> hi;
for (j = 1; j <= 3; j++)
for (int j = 1; j <= 3; j++)
{
in >> el.PNum(j);
// el.PNum(j)++;
@ -93,11 +88,11 @@ namespace netgen
if ( (strlen (filename) > 4) &&
strcmp (&filename[strlen (filename)-4], ".unv") == 0 )
{
int i, j, k;
// int i, j, k;
double h;
// double h;
char reco[100];
int np, nbe;
// int np, nbe;
int invert;
@ -116,13 +111,13 @@ namespace netgen
if (strcmp (reco, "NODES") == 0)
{
cout << "nodes found" << endl;
for (j = 1; j <= 4; j++)
for (int j = 1; j <= 4; j++)
in >> reco; // read dummy
while (1)
{
int pi, hi;
double x, y, z;
// double x, y, z;
Point3d p;
in >> pi;
@ -144,7 +139,7 @@ namespace netgen
if (strcmp (reco, "ELEMENTS") == 0)
{
cout << "elements found" << endl;
for (j = 1; j <= 4; j++)
for (int j = 1; j <= 4; j++)
in >> reco; // read dummy
while (1)
@ -152,7 +147,7 @@ namespace netgen
int hi;
in >> hi;
if (hi == -1) break;
for (j = 1; j <= 7; j++)
for (int j = 1; j <= 7; j++)
in >> hi;
Element2d el;
@ -163,7 +158,7 @@ namespace netgen
swap (el.PNum(2), el.PNum(3));
mesh.AddSurfaceElement (el);
for (j = 1; j <= 5; j++)
for (int j = 1; j <= 5; j++)
in >> hi;
}
}
@ -260,7 +255,7 @@ namespace netgen
in >> nse;
for (i = 1; i <= nse; i++)
{
int mat, nelp;
int mat; // , nelp;
in >> mat;
Element2d el (TRIG);
el.SetIndex (mat);
@ -301,7 +296,7 @@ namespace netgen
cout << "pktfile = " << pktfile << endl;
int np, nse, i;
int num, bcprop;
int bcprop;
ifstream inpkt (pktfile.c_str());
inpkt >> np;
Array<double> values(np);

12
libsrc/interface/writedolfin.cpp
View File

@ -25,14 +25,14 @@ namespace netgen
int np = mesh.GetNP();
int ne = mesh.GetNE();
int nse = mesh.GetNSE();
// int nse = mesh.GetNSE();
int nsd = mesh.GetDimension();
int invertsurf = mparam.inverttrigs;
int i, j;
// int invertsurf = mparam.inverttrigs;
// int i, j;
ofstream outfile (filename.c_str());
char str[100];
// char str[100];
outfile.precision(8);
outfile.setf (ios::fixed, ios::floatfield);
outfile.setf (ios::showpoint);
@ -45,7 +45,7 @@ namespace netgen
outfile << "<dolfin xmlns:dolfin=\"http://www.phi.chalmers.se/dolfin/\">"<<endl;
outfile << " <mesh celltype=\"tetrahedron\" dim=\"3\">" <<endl;
outfile << " <vertices size=\""<<np<<"\">"<<endl;
for (i = 1; i <= np; i++) {
for (int i = 1; i <= np; i++) {
const Point3d & p = mesh.Point(i);
outfile << " <vertex index=\""<<i-1<<"\" x=\""<<p.X()<<"\" y=\""<<p.Y()<<"\" z=\""<<p.Z()<<"\"/>"<<endl;
}
@ -54,7 +54,7 @@ namespace netgen
outfile << " <cells size=\""<<ne<<"\">"<<endl;
for (i = 1; i <= ne; i++) {
for (int i = 1; i <= ne; i++) {
const Element & el = mesh.VolumeElement(i);
outfile << " <tetrahedron index=\""<<i-1<<"\" v0=\""<<el.PNum(1)-1<<"\" v1=\""<<el.PNum(2)-1<<"\" v2=\""<<el.PNum(3)-1<<"\" v3=\""<<el.PNum(4)-1<<"\"/>"<<endl;

3
libsrc/interface/writeelmer.cpp
View File

@ -37,7 +37,8 @@ void WriteElmerFormat (const Mesh &mesh,
sprintf( a, "mkdir %s", filename.c_str() );
system( a );
#else
int rc = mkdir(filename.c_str(), S_IRWXU|S_IRWXG);
// int rc =
mkdir(filename.c_str(), S_IRWXU|S_IRWXG);
#endif
sprintf( str, "%s/mesh.header", filename.c_str() );

494
libsrc/interface/wuchemnitz.cpp
View File

@ -12,302 +12,306 @@
namespace netgen
{
class POINT3D
class POINT3D
{
public:
POINT3D () { };
double x, y, z;
POINT3D () { };
double x, y, z;
};
class VOLELEMENT
class VOLELEMENT
{
public:
VOLELEMENT () {};
int domnr, p1, p2, p3, p4;
int faces[4];
int domnr, p1, p2, p3, p4;
int faces[4];
VOLELEMENT ()
{ for (int i = 0; i < 4; i++) faces[i] = 0; }
};
class SURFELEMENT
class SURFELEMENT
{
public:
SURFELEMENT () { };
int snr, p1, p2, p3;
SURFELEMENT () { };
int snr, p1, p2, p3;
};
class FACE
class FACE
{
public:
FACE () { };
int p1, p2, p3;
int edges[3];
int p1, p2, p3;
int edges[3];
FACE ()
{ for (int i = 0; i < 3; i++) edges[i] = 0; }
};
class EDGE
class EDGE
{
public:
EDGE () { };
int p1, p2;
EDGE () { };
int p1, p2;
};
static Array<POINT3D> points;
static Array<VOLELEMENT> volelements;
static Array<SURFELEMENT> surfelements;
static Array<POINT3D> points;
static Array<VOLELEMENT> volelements;
static Array<SURFELEMENT> surfelements;
static Array<FACE> faces;
static Array<EDGE> edges;
static Array<FACE> faces;
static Array<EDGE> edges;
void ReadFile (char * filename)
void ReadFile (char * filename)
{
int i, n;
ifstream infile(filename);
char reco[100];
int i, n;
ifstream infile(filename);
char reco[100];
infile >> reco; // file format recognition
infile >> reco; // file format recognition
infile >> n; // number of surface elements
cout << n << " Surface elements" << endl;
infile >> n; // number of surface elements
cout << n << " Surface elements" << endl;
for (i = 1; i <= n; i++)
{
SURFELEMENT sel;
infile >> sel.snr >> sel.p1 >> sel.p2 >> sel.p3;
surfelements.Append (sel);
}
infile >> n; // number of volume elements
cout << n << " Volume elements" << endl;
for (i = 1; i <= n; i++)
{
VOLELEMENT el;
infile >> el.p1 >> el.p2 >> el.p3 >> el.p4;
volelements.Append (el);
}
infile >> n; // number of points
cout << n << " Points" << endl;
for (i = 1; i <= n; i++)
{
POINT3D p;
infile >> p.x >> p.y >> p.z;
points.Append (p);
}
}
void ReadFileMesh (const Mesh & mesh)
{
int i, n;
n = mesh.GetNSE(); // number of surface elements
cout << n << " Surface elements" << endl;
for (i = 1; i <= n; i++)
{
SURFELEMENT sel;
const Element2d & el = mesh.SurfaceElement(i);
sel.snr = el.GetIndex();
sel.p1 = el.PNum(1);
sel.p2 = el.PNum(2);
sel.p3 = el.PNum(3);
surfelements.Append (sel);
}
n = mesh.GetNE(); // number of volume elements
cout << n << " Volume elements" << endl;
for (i = 1; i <= n; i++)
{
VOLELEMENT el;
const Element & nel = mesh.VolumeElement(i);
el.p1 = nel.PNum(1);
el.p2 = nel.PNum(2);
el.p3 = nel.PNum(3);
el.p4 = nel.PNum(4);
// infile >> el.p1 >> el.p2 >> el.p3 >> el.p4;
volelements.Append (el);
}
n = mesh.GetNP(); // number of points
cout << n << " Points" << endl;
for (i = 1; i <= n; i++)
{
POINT3D p;
Point3d mp = mesh.Point(i);
p.x = mp.X();
p.y = mp.Y();
p.z = mp.Z();
// infile >> p.x >> p.y >> p.z;
points.Append (p);
}
}
void Convert ()
{
int i, j, facei, edgei;
INDEX_3 i3;
INDEX_2 i2;
INDEX_3_HASHTABLE<int> faceindex(volelements.Size()/5 + 1);
INDEX_2_HASHTABLE<int> edgeindex(volelements.Size()/5 + 1);
for (i = 1; i <= volelements.Size(); i++)
{
for (j = 1; j <= 4; j++)
for (i = 1; i <= n; i++)
{
switch (j)
{
case 1:
i3.I1() = volelements.Get(i).p2;
i3.I2() = volelements.Get(i).p3;
i3.I3() = volelements.Get(i).p4;
break;
case 2:
i3.I1() = volelements.Get(i).p1;
i3.I2() = volelements.Get(i).p3;
i3.I3() = volelements.Get(i).p4;
break;
case 3:
i3.I1() = volelements.Get(i).p1;
i3.I2() = volelements.Get(i).p2;
i3.I3() = volelements.Get(i).p4;
break;
case 4:
i3.I1() = volelements.Get(i).p1;
i3.I2() = volelements.Get(i).p2;
i3.I3() = volelements.Get(i).p3;
break;
default:
i3.I1()=i3.I2()=i3.I3()=0;
}
i3.Sort();
if (faceindex.Used (i3))
facei = faceindex.Get(i3);
else
{
FACE fa;
fa.p1 = i3.I1();
fa.p2 = i3.I2();
fa.p3 = i3.I3();
facei = faces.Append (fa);
faceindex.Set (i3, facei);
}
volelements.Elem(i).faces[j-1] = facei;
}
SURFELEMENT sel;
infile >> sel.snr >> sel.p1 >> sel.p2 >> sel.p3;
surfelements.Append (sel);
}
}
infile >> n; // number of volume elements
cout << n << " Volume elements" << endl;
for (i = 1; i <= n; i++)
{
VOLELEMENT el;
infile >> el.p1 >> el.p2 >> el.p3 >> el.p4;
volelements.Append (el);
}
infile >> n; // number of points
cout << n << " Points" << endl;
for (i = 1; i <= n; i++)
{
POINT3D p;
infile >> p.x >> p.y >> p.z;
points.Append (p);
}
}
void ReadFileMesh (const Mesh & mesh)
{
int i, n;
n = mesh.GetNSE(); // number of surface elements
cout << n << " Surface elements" << endl;
for (i = 1; i <= n; i++)
{
SURFELEMENT sel;
const Element2d & el = mesh.SurfaceElement(i);
sel.snr = el.GetIndex();
sel.p1 = el.PNum(1);
sel.p2 = el.PNum(2);
sel.p3 = el.PNum(3);
surfelements.Append (sel);
}
n = mesh.GetNE(); // number of volume elements
cout << n << " Volume elements" << endl;
for (i = 1; i <= n; i++)
{
VOLELEMENT el;
const Element & nel = mesh.VolumeElement(i);
el.p1 = nel.PNum(1);
el.p2 = nel.PNum(2);
el.p3 = nel.PNum(3);
el.p4 = nel.PNum(4);
// infile >> el.p1 >> el.p2 >> el.p3 >> el.p4;
volelements.Append (el);
}
n = mesh.GetNP(); // number of points
cout << n << " Points" << endl;
for (i = 1; i <= n; i++)
{
POINT3D p;
Point3d mp = mesh.Point(i);
p.x = mp.X();
p.y = mp.Y();
p.z = mp.Z();
// infile >> p.x >> p.y >> p.z;
points.Append (p);
}
}
void Convert ()
{
int i, j, facei, edgei;
INDEX_3 i3;
INDEX_2 i2;
INDEX_3_HASHTABLE<int> faceindex(volelements.Size()/5 + 1);
INDEX_2_HASHTABLE<int> edgeindex(volelements.Size()/5 + 1);
for (i = 1; i <= volelements.Size(); i++)
{
for (j = 1; j <= 4; j++)
{
switch (j)
{
case 1:
i3.I1() = volelements.Get(i).p2;
i3.I2() = volelements.Get(i).p3;
i3.I3() = volelements.Get(i).p4;
break;
case 2:
i3.I1() = volelements.Get(i).p1;
i3.I2() = volelements.Get(i).p3;
i3.I3() = volelements.Get(i).p4;
break;
case 3:
i3.I1() = volelements.Get(i).p1;
i3.I2() = volelements.Get(i).p2;
i3.I3() = volelements.Get(i).p4;
break;
case 4:
i3.I1() = volelements.Get(i).p1;
i3.I2() = volelements.Get(i).p2;
i3.I3() = volelements.Get(i).p3;
break;
default:
i3.I1()=i3.I2()=i3.I3()=0;
}
i3.Sort();
if (faceindex.Used (i3))
facei = faceindex.Get(i3);
else
{
FACE fa;
fa.p1 = i3.I1();
fa.p2 = i3.I2();
fa.p3 = i3.I3();
facei = faces.Append (fa);
faceindex.Set (i3, facei);
}
volelements.Elem(i).faces[j-1] = facei;
}
}
for (i = 1; i <= faces.Size(); i++)
{
for (j = 1; j <= 3; j++)
for (i = 1; i <= faces.Size(); i++)
{
switch (j)
{
case 1:
i2.I1() = faces.Get(i).p2;
i2.I2() = faces.Get(i).p3;
break;
case 2:
i2.I1() = faces.Get(i).p1;
i2.I2() = faces.Get(i).p3;
break;
case 3:
i2.I1() = faces.Get(i).p1;
i2.I2() = faces.Get(i).p2;
break;
default:
i2.I1()=i2.I2()=0;
}
if (i2.I1() > i2.I2()) swap (i2.I1(), i2.I2());
if (edgeindex.Used (i2))
edgei = edgeindex.Get(i2);
else
{
EDGE ed;
ed.p1 = i2.I1();
ed.p2 = i2.I2();
edgei = edges.Append (ed);
edgeindex.Set (i2, edgei);
}
for (j = 1; j <= 3; j++)
{
switch (j)
{
case 1:
i2.I1() = faces.Get(i).p2;
i2.I2() = faces.Get(i).p3;
break;
case 2:
i2.I1() = faces.Get(i).p1;
i2.I2() = faces.Get(i).p3;
break;
case 3:
i2.I1() = faces.Get(i).p1;
i2.I2() = faces.Get(i).p2;
break;
default:
i2.I1()=i2.I2()=0;
}
if (i2.I1() > i2.I2()) swap (i2.I1(), i2.I2());
if (edgeindex.Used (i2))
edgei = edgeindex.Get(i2);
else
{
EDGE ed;
ed.p1 = i2.I1();
ed.p2 = i2.I2();
edgei = edges.Append (ed);
edgeindex.Set (i2, edgei);
}
faces.Elem(i).edges[j-1] = edgei;
}
faces.Elem(i).edges[j-1] = edgei;
}
}
}
}
void WriteFile (ostream & outfile)
void WriteFile (ostream & outfile)
{
int i;
int i;
outfile
<< "#VERSION: 1.0" << endl
<< "#PROGRAM: NETGEN" << endl
<< "#EQN_TYPE: POISSON" << endl
<< "#DIMENSION: 3D" << endl
<< "#DEG_OF_FREE: 1" << endl
<< "#DESCRIPTION: I don't know" << endl
<< "##RENUM: not done" << endl
<< "#USER: Kleinzen" << endl
<< "DATE: 10.06.1996" << endl;
outfile
<< "#VERSION: 1.0" << endl
<< "#PROGRAM: NETGEN" << endl
<< "#EQN_TYPE: POISSON" << endl
<< "#DIMENSION: 3D" << endl
<< "#DEG_OF_FREE: 1" << endl
<< "#DESCRIPTION: I don't know" << endl
<< "##RENUM: not done" << endl
<< "#USER: Kleinzen" << endl
<< "DATE: 10.06.1996" << endl;
outfile << "#HEADER: 8" << endl
<< points.Size() << " " << edges.Size() << " "
<< faces.Size() << " " << volelements.Size() << " 0 0 0 0" << endl;
outfile << "#HEADER: 8" << endl
<< points.Size() << " " << edges.Size() << " "
<< faces.Size() << " " << volelements.Size() << " 0 0 0 0" << endl;
outfile << "#VERTEX: " << points.Size() << endl;
for (i = 1; i <= points.Size(); i++)
outfile << " " << i << " " << points.Get(i).x << " " << points.Get(i).y
<< " " << points.Get(i).z << endl;
outfile << "#VERTEX: " << points.Size() << endl;
for (i = 1; i <= points.Size(); i++)
outfile << " " << i << " " << points.Get(i).x << " " << points.Get(i).y
<< " " << points.Get(i).z << endl;
outfile << "#EDGE: " << edges.Size() << endl;
for (i = 1; i <= edges.Size(); i++)
outfile << " " << i << " 1 "
<< edges.Get(i).p1 << " "
<< edges.Get(i).p2
<< " 0" << endl;
outfile << "#EDGE: " << edges.Size() << endl;
for (i = 1; i <= edges.Size(); i++)
outfile << " " << i << " 1 "
<< edges.Get(i).p1 << " "
<< edges.Get(i).p2
<< " 0" << endl;
outfile << "#FACE: " << faces.Size() << endl;
for (i = 1; i <= faces.Size(); i++)
outfile << " " << i << " 1 3 "
<< faces.Get(i).edges[0] << " "
<< faces.Get(i).edges[1] << " "
<< faces.Get(i).edges[2] << endl;
outfile << "#FACE: " << faces.Size() << endl;
for (i = 1; i <= faces.Size(); i++)
outfile << " " << i << " 1 3 "
<< faces.Get(i).edges[0] << " "
<< faces.Get(i).edges[1] << " "
<< faces.Get(i).edges[2] << endl;
outfile << "#SOLID: " << volelements.Size() << endl;
for (i = 1; i <= volelements.Size(); i++)
outfile << " " << i << " 1 4 "
<< volelements.Get(i).faces[0] << " "
<< volelements.Get(i).faces[1] << " "
<< volelements.Get(i).faces[2] << " "
<< volelements.Get(i).faces[3] << endl;
outfile << "#SOLID: " << volelements.Size() << endl;
for (i = 1; i <= volelements.Size(); i++)
outfile << " " << i << " 1 4 "
<< volelements.Get(i).faces[0] << " "
<< volelements.Get(i).faces[1] << " "
<< volelements.Get(i).faces[2] << " "
<< volelements.Get(i).faces[3] << endl;
outfile << "#END_OF_DATA" << endl;
outfile << "#END_OF_DATA" << endl;
}
void WriteUserChemnitz (const Mesh & mesh,
const string & filename)
{
ofstream outfile (filename.c_str());
void WriteUserChemnitz (const Mesh & mesh,
const string & filename)
{
ofstream outfile (filename.c_str());
ReadFileMesh (mesh);
Convert ();
ReadFileMesh (mesh);
Convert ();
WriteFile (outfile);
cout << "Wrote Chemnitz standard file" << endl;
}
WriteFile (outfile);
cout << "Wrote Chemnitz standard file" << endl;
}
}

16
libsrc/meshing/bisect.cpp
View File

@ -362,6 +362,8 @@ namespace netgen
ned = 6;
break;
}
default:
throw NgException("Bisect, element type not handled in switch");
}
for (j = 0; j < ned; j++)
@ -462,7 +464,7 @@ namespace netgen
{ { 1, 2, 4, 3 },
{ 1, 5, 4, 5 },
{ 2, 5, 3, 5 } };
int (*pairs)[4] = NULL;
switch (el.GetType())
{
@ -477,6 +479,8 @@ namespace netgen
pairs = pyramidpairs;
break;
}
default:
throw NgException("Bisect, element type not handled in switch, 2");
}
for (j = 0; j < 3; j++)
@ -744,6 +748,8 @@ namespace netgen
ned = 6;
break;
}
default:
throw NgException("Bisect, element type not handled in switch, 3");
}
for (j = 0; j < ned; j++)
@ -802,6 +808,8 @@ namespace netgen
pairs = pyramidpairs;
break;
}
default:
throw NgException("Bisect, element type not handled in switch, 3a");
}
for (j = 0; j < 3; j++)
@ -2100,6 +2108,8 @@ namespace netgen
mprisms.Append (mp);
break;
}
default:
throw NgException("Bisect, element type not handled in switch, 4");
}
}
@ -2523,6 +2533,8 @@ namespace netgen
<< endl;
break;
}
default:
throw NgException("Bisect, element type not handled in switch, 6");
}
}
@ -2563,6 +2575,8 @@ namespace netgen
mquads.Append (mt);
break;
}
default:
throw NgException("Bisect, element type not handled in switch, 5");
}

4
libsrc/meshing/classifyhpel.hpp
View File

@ -792,6 +792,7 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edge
if(cornerpoint.Test(el.PNum(p[k])))
point_sing[p[k]-1] = 3;
*testout << "point_sing = " << point_sing[0] << point_sing[1] << point_sing[2] << endl;
if(edge_sing[0] + edge_sing[1] + edge_sing[2] == 0)
{
@ -851,8 +852,11 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edge
// cout << " run for " << j << " gives type " << type << endl;
//*testout << " run for " << j << " gives type " << type << endl;
if(type!=HP_NONE) break;
}
*testout << "type = " << type << endl;
for(int k=0;k<3;k++) el[k] = pnums[k];
/*if(type != HP_NONE)

15
libsrc/meshing/curvedelems.cpp
View File

@ -367,7 +367,6 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta,
AutoDiff<3> adt(t, 2);
AutoDiff<3> res[2000];
CalcScaledTrigShape (n, adx, ady, adt, &res[0]);
double dshape1[6000];
int ndof = (n-1)*(n-2)/2;
for (int i = 0; i < ndof; i++)
{
@ -377,6 +376,7 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta,
}
/*
double dshape1[6000];
if (n < 3) return;
double hvl[1000], hvr[1000];
int nd = (n-1)*(n-2)/2;
@ -1486,6 +1486,9 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta,
break;
}
default:
throw NgException("CurvedElements::CalcShape 2d, element type not handled");
};
}
@ -1746,6 +1749,9 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta,
*/
break;
}
default:
throw NgException("CurvedElements::CalcDShape 2d, element type not handled");
};
}
@ -2189,6 +2195,10 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta,
shapes[7] = (1-x)* y *(z);
break;
}
default:
throw NgException("CurvedElements::CalcShape 3d, element type not handled");
};
}
@ -2616,6 +2626,9 @@ inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta,
break;
}
default:
throw NgException("CurvedElements::CalcDShape 3d, element type not handled");
}
/*

2
libsrc/meshing/findip.hpp
View File

@ -78,7 +78,7 @@ inline int FindInnerPoint (POINTArray & points,
Array<Vec3d> a;
Array<double> c;
Mat<3> m, inv;
Vec<3> rs, x, center;
Vec<3> rs, x = 0.0, center;
double f;
int nf = faces.Size();

7357
libsrc/meshing/meshclass.cpp
View File

File diff suppressed because it is too large Load Diff

13
libsrc/meshing/meshing2.cpp
View File

@ -13,7 +13,7 @@ namespace netgen
static Array<Point2d> plainpoints;
static Array<int> plainzones;
static Array<INDEX_2> loclines;
static int geomtrig;
// static int geomtrig;
//static const char * rname;
static int cntelem, trials, nfaces;
static int oldnl;
@ -400,13 +400,10 @@ namespace netgen
debugflag =
debugparam.haltsegment &&
( (debugparam.haltsegmentp1 == gpi1) &&
(debugparam.haltsegmentp2 == gpi2) ||
(debugparam.haltsegmentp1 == gpi2) &&
(debugparam.haltsegmentp2 == gpi1)) ||
( ((debugparam.haltsegmentp1 == gpi1) && (debugparam.haltsegmentp2 == gpi2)) ||
((debugparam.haltsegmentp1 == gpi2) && (debugparam.haltsegmentp2 == gpi1))) ||
debugparam.haltnode &&
( (debugparam.haltsegmentp1 == gpi1) ||
(debugparam.haltsegmentp2 == gpi1));
( (debugparam.haltsegmentp1 == gpi1) || (debugparam.haltsegmentp2 == gpi1));
if (debugparam.haltface && debugparam.haltfacenr == facenr)
@ -566,7 +563,7 @@ namespace netgen
}
}
else if (z1 > 0 && z2 > 0 && (z1 != z2) || (z1 < 0) && (z2 < 0) )
else if ( (z1 > 0 && z2 > 0 && (z1 != z2)) || ((z1 < 0) && (z2 < 0)) )
{
loclines.DeleteElement(i);
lindex.DeleteElement(i);

2
libsrc/stlgeom/stlgeommesh.cpp
View File

@ -678,7 +678,7 @@ int STLGeometry :: ProjectOnWholeSurface(Point<3> & p3d) const
int STLGeometry :: ProjectNearest(Point<3> & p3d) const
{
Point<3> p, pf;
Point<3> p, pf = 0.0;
//set new chart
const STLChart& chart = GetChart(meshchart);

6
libsrc/stlgeom/stltool.cpp
View File

@ -569,12 +569,10 @@ double STLTriangle :: GetNearestPoint(const Array<Point<3> >& ap,
if (PointInside(ap, p)) {p3d = p; return dist;}
else
{
Point<3> pf;
Point<3> pf = 0.0;
double nearest = 1E50;
//int fi = 0;
int j;
for (j = 1; j <= 3; j++)
for (int j = 1; j <= 3; j++)
{
p = p3d;
dist = GetDistFromLine(ap.Get(PNum(j)), ap.Get(PNumMod(j+1)), p);

199
libsrc/visualization/vssolution.cpp
View File

@ -2166,6 +2166,7 @@ namespace netgen
}
else if (el.GetType() == QUAD)
{
/*
Array < Point<3> > lp(3);
lp[0] = mesh->Point(el[0]);
@ -2179,106 +2180,106 @@ namespace netgen
lp[2] = mesh->Point(el[3]);
DrawTrigSurfaceVectors(lp,pmin,pmax,sei,vsol);
/*
Point<3> lp[4];
Point<2> p2d[4];
for (k = 0; k < 4; k++)
lp[k] = mesh->Point (el[k]);
Vec<3> n = Cross (lp[1]-lp[0], lp[2]-lp[0]);
Vec<3> na (fabs (n(0)), fabs(n(1)), fabs(n(2)));
if (na(0) > na(1) && na(0) > na(2))
dir = 1;
else if (na(1) > na(2))
dir = 2;
else
dir = 3;
dir1 = (dir % 3) + 1;
dir2 = (dir1 % 3) + 1;
for (k = 0; k < 4; k++)
{
p2d[k] = Point<2> ((lp[k](dir1-1) - pmin(dir1-1)) / (2*rad),
(lp[k](dir2-1) - pmin(dir2-1)) / (2*rad));
}
double minx2d, maxx2d, miny2d, maxy2d;
minx2d = maxx2d = p2d[0](0);
miny2d = maxy2d = p2d[0](1);
for (k = 1; k < 4; k++)
{
minx2d = min2 (minx2d, p2d[k](0));
maxx2d = max2 (maxx2d, p2d[k](0));
miny2d = min2 (miny2d, p2d[k](1));
maxy2d = max2 (maxy2d, p2d[k](1));
}
for (s = xoffset/gridsize; s <= 1+xoffset/gridsize; s += 1.0 / gridsize)
if (s >= minx2d && s <= maxx2d)
for (t = yoffset/gridsize; t <= 1+yoffset/gridsize; t += 1.0 / gridsize)
if (t >= miny2d && t <= maxy2d)
{
double lami[3];
Point3d p3d(2*rad*s+pmin(0), 2*rad*t+pmin(1),0);
if (mesh->PointContainedIn2DElement (p3d, lami, sei+1))
{
Point<3> cp = p3d;
double lam1 = lami[0];
double lam2 = lami[1];
//for (k = 0; k < 3; k++)
//cp(k) = lp[0](k) +
//lam1 * (lp[1](k)-lp[0](k)) +
//lam2 * (lp[2](k)-lp[0](k));
Vec<3> v;
double values[6];
drawelem = GetSurfValues (vsol, sei, lam1, lam2, values);
(*testout) << "sei " << sei << " lam1 " << lam1 << " lam2 " << lam2 << " drawelem " << drawelem << endl;
if (!vsol->iscomplex)
for (k = 0; k < 3; k++)
v(k) = values[k];
else
{
if (!imag_part)
for (k = 0; k < 3; k++)
v(k) = values[2*k];
else
for (k = 0; k < 3; k++)
v(k) = values[2*k+1];
}
if (mesh->GetDimension() == 2)
if ( (!vsol->iscomplex && vsol->components != 3) ||
(vsol->iscomplex && vsol->components != 6) )
v(2) = 0;
double val = v.Length();
SetOpenGlColor (val, minval, maxval, logscale);
(*testout) << "v " << v << endl;
if (val > 1e-10 * maxval)
v *= (rad / val / gridsize * 0.5);
(*testout) << "v " << v << endl;
if ( drawelem )
{
DrawCone (cp, cp+4*v, 0.8*rad / gridsize);
(*testout) << "cp " << cp << " rad " << rad << " gridsize " << gridsize << endl;
}
}
}
*/
Point<3> lp[4];
Point<2> p2d[4];
for (int k = 0; k < 4; k++)
lp[k] = mesh->Point (el[k]);
Vec<3> n = Cross (lp[1]-lp[0], lp[2]-lp[0]);
Vec<3> na (fabs (n(0)), fabs(n(1)), fabs(n(2)));
int dir, dir1, dir2;
if (na(0) > na(1) && na(0) > na(2))
dir = 1;
else if (na(1) > na(2))
dir = 2;
else
dir = 3;
dir1 = (dir % 3) + 1;
dir2 = (dir1 % 3) + 1;
for (int k = 0; k < 4; k++)
{
p2d[k] = Point<2> ((lp[k](dir1-1) - pmin(dir1-1)) / (2*rad),
(lp[k](dir2-1) - pmin(dir2-1)) / (2*rad));
}
double minx2d, maxx2d, miny2d, maxy2d;
minx2d = maxx2d = p2d[0](0);
miny2d = maxy2d = p2d[0](1);
for (int k = 1; k < 4; k++)
{
minx2d = min2 (minx2d, p2d[k](0));
maxx2d = max2 (maxx2d, p2d[k](0));
miny2d = min2 (miny2d, p2d[k](1));
maxy2d = max2 (maxy2d, p2d[k](1));
}
for (double s = xoffset/gridsize; s <= 1+xoffset/gridsize; s += 1.0 / gridsize)
if (s >= minx2d && s <= maxx2d)
for (double t = yoffset/gridsize; t <= 1+yoffset/gridsize; t += 1.0 / gridsize)
if (t >= miny2d && t <= maxy2d)
{
double lami[3];
Point3d p3d(2*rad*s+pmin(0), 2*rad*t+pmin(1),0);
if (mesh->PointContainedIn2DElement (p3d, lami, sei+1))
{
Point<3> cp = p3d;
double lam1 = lami[0];
double lam2 = lami[1];
//for (k = 0; k < 3; k++)
//cp(k) = lp[0](k) +
//lam1 * (lp[1](k)-lp[0](k)) +
//lam2 * (lp[2](k)-lp[0](k));
Vec<3> v;
double values[6];
bool drawelem = GetSurfValues (vsol, sei, lam1, lam2, values);
(*testout) << "sei " << sei << " lam1 " << lam1 << " lam2 " << lam2 << " drawelem " << drawelem << endl;
if (!vsol->iscomplex)
for (int k = 0; k < 3; k++)
v(k) = values[k];
else
{
if (!imag_part)
for (int k = 0; k < 3; k++)
v(k) = values[2*k];
else
for (int k = 0; k < 3; k++)
v(k) = values[2*k+1];
}
if (mesh->GetDimension() == 2)
if ( (!vsol->iscomplex && vsol->components != 3) ||
(vsol->iscomplex && vsol->components != 6) )
v(2) = 0;
double val = v.Length();
SetOpenGlColor (val, minval, maxval, logscale);
(*testout) << "v " << v << endl;
if (val > 1e-10 * maxval)
v *= (rad / val / gridsize * 0.5);
(*testout) << "v " << v << endl;
if ( drawelem )
{
DrawCone (cp, cp+4*v, 0.8*rad / gridsize);
(*testout) << "cp " << cp << " rad " << rad << " gridsize " << gridsize << endl;
}
}
}
}
}
}