netgen/libsrc/csg/csgeom.cpp

1563 lines
34 KiB
C++
Raw Normal View History

2009-01-13 04:40:13 +05:00
#include <mystdlib.h>
#include <myadt.hpp>
#include <linalg.hpp>
#include <csg.hpp>
namespace netgen
{
2011-01-11 01:18:01 +05:00
2009-01-13 04:40:13 +05:00
int CSGeometry :: changeval = 0;
TopLevelObject ::
TopLevelObject (Solid * asolid,
Surface * asurface)
{
solid = asolid;
surface = asurface;
SetRGB (0, 0, 1);
SetTransparent (0);
SetVisible (1);
SetLayer (1);
if (!surface)
maxh = solid->GetMaxH();
else
maxh = surface->GetMaxH();
SetBCProp (-1);
bcname = "default";
}
void TopLevelObject :: GetData (ostream & ost)
{
ost << red << " " << green << " " << blue << " "
<< transp << " " << visible << " ";
}
void TopLevelObject :: SetData (istream & ist)
{
ist >> red >> green >> blue >> transp >> visible;
}
Box<3> CSGeometry::default_boundingbox (Point<3> (-1000, -1000, -1000),
Point<3> ( 1000, 1000, 1000));
CSGeometry :: CSGeometry ()
: boundingbox (default_boundingbox),
2009-04-20 03:15:26 +06:00
identicsurfaces (100), ideps(1e-9), filename("")
2009-01-13 04:40:13 +05:00
{
;
}
CSGeometry :: CSGeometry (const string & afilename)
: boundingbox (default_boundingbox),
2009-04-20 03:15:26 +06:00
identicsurfaces (100), ideps(1e-9), filename(afilename)
2009-01-13 04:40:13 +05:00
{
changeval++;
}
CSGeometry :: ~CSGeometry ()
{
Clean();
}
void CSGeometry :: Clean ()
{
2009-01-25 17:35:25 +05:00
Array< Solid* > to_delete;
2009-01-13 04:40:13 +05:00
for (int i = 0; i < solids.Size(); i++)
if(!to_delete.Contains(solids[i]->S1()))
to_delete.Append(solids[i]->S1());
for (int i = 0; i < solids.Size(); i++)
if(!to_delete.Contains(solids[i]))
to_delete.Append(solids[i]);
2013-02-27 19:49:17 +06:00
2009-01-13 04:40:13 +05:00
for(int i = 0; i < to_delete.Size(); i++)
delete to_delete[i];
solids.DeleteAll ();
for (int i = 0; i < splinecurves2d.Size(); i++)
delete splinecurves2d[i];
splinecurves2d.DeleteAll();
/*
for (int i = 0; i < surfaces.Size(); i++)
delete surfaces[i];
surfaces.DeleteAll ();
*/
for(int i = 0; i<delete_them.Size(); i++)
delete delete_them[i];
delete_them.DeleteAll();
surfaces.DeleteAll();
2013-02-27 19:49:17 +06:00
2009-01-13 04:40:13 +05:00
for (int i = 0; i < toplevelobjects.Size(); i++)
delete toplevelobjects[i];
toplevelobjects.DeleteAll ();
for (int i = 0; i < triapprox.Size(); i++)
delete triapprox[i];
triapprox.DeleteAll();
for(int i = 0; i < identifications.Size(); i++)
delete identifications[i];
identifications.DeleteAll();
for (int i = 0; i < singfaces.Size(); i++)
delete singfaces[i];
singfaces.DeleteAll();
for (int i = 0; i < singedges.Size(); i++)
delete singedges[i];
singedges.DeleteAll();
for (int i = 0; i < singpoints.Size(); i++)
delete singpoints[i];
singpoints.DeleteAll();
changeval++;
}
2011-01-11 01:18:01 +05:00
extern int CSGGenerateMesh (CSGeometry & geom,
2014-09-26 02:23:31 +06:00
shared_ptr<Mesh> & mesh, MeshingParameters & mparam,
2011-01-11 01:18:01 +05:00
int perfstepsstart, int perfstepsend);
2009-08-25 20:00:20 +06:00
2011-01-11 01:18:01 +05:00
2014-09-26 02:23:31 +06:00
int CSGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam,
2011-01-11 01:18:01 +05:00
int perfstepsstart, int perfstepsend)
2009-08-25 20:00:20 +06:00
{
2011-01-11 01:18:01 +05:00
return CSGGenerateMesh (*this, mesh, mparam, perfstepsstart, perfstepsend);
2009-08-25 20:00:20 +06:00
}
2009-01-13 04:40:13 +05:00
2009-08-25 20:00:20 +06:00
const Refinement & CSGeometry :: GetRefinement () const
{
2011-01-15 10:10:52 +05:00
// cout << "get CSGeometry - Refinement" << endl;
2009-08-25 20:00:20 +06:00
// should become class variables
RefinementSurfaces * ref = new RefinementSurfaces(*this);
ref -> Set2dOptimizer(new MeshOptimize2dSurfaces(*this));
return *ref;
}
2009-01-13 04:40:13 +05:00
class WritePrimitivesIt : public SolidIterator
{
ostream & ost;
public:
WritePrimitivesIt (ostream & aost) : ost(aost) { ; }
virtual ~WritePrimitivesIt () { ; }
virtual void Do (Solid * sol);
};
void WritePrimitivesIt :: Do (Solid * sol)
{
Primitive * prim = sol->GetPrimitive();
if (prim)
{
const char * classname;
2009-01-25 17:35:25 +05:00
Array<double> coeffs;
2009-01-13 04:40:13 +05:00
prim -> GetPrimitiveData (classname, coeffs);
if (sol->Name())
ost << "primitive "
<< sol->Name() << " "
<< classname << " " << coeffs.Size();
for (int i = 0; i < coeffs.Size(); i++)
ost << " " << coeffs[i];
ost << endl;
}
}
2011-01-11 01:18:01 +05:00
void CSGeometry :: Save (string filename) const
{
2014-09-26 02:23:31 +06:00
ofstream ost (filename.c_str());
2011-01-11 01:18:01 +05:00
Save (ost);
}
2009-01-13 04:40:13 +05:00
2011-01-11 01:18:01 +05:00
void CSGeometry :: Save (ostream & ost) const
2009-01-13 04:40:13 +05:00
{
ost << "boundingbox "
<< boundingbox.PMin()(0) << " "
<< boundingbox.PMin()(1) << " "
<< boundingbox.PMin()(2) << " "
<< boundingbox.PMax()(0) << " "
<< boundingbox.PMax()(1) << " "
<< boundingbox.PMax()(2) << endl;
WritePrimitivesIt wpi(ost);
IterateAllSolids (wpi, 1);
for (int i = 0; i < solids.Size(); i++)
{
if (!solids[i]->GetPrimitive())
{
ost << "solid " << solids.GetName(i) << " ";
solids[i] -> GetSolidData (ost);
ost << endl;
}
}
for (int i = 0; i < GetNTopLevelObjects(); i++)
{
TopLevelObject * tlo = GetTopLevelObject (i);
ost << "toplevel ";
if (tlo -> GetSurface())
ost << "surface " << tlo->GetSolid()->Name() << " "
<< tlo->GetSurface()->Name() << " ";
else
ost << "solid " << tlo->GetSolid()->Name() << " ";
tlo->GetData(ost);
ost << endl;
}
for (int i = 0; i < identifications.Size(); i++)
{
ost << "identify ";
identifications[i] -> GetData (ost);
ost << endl;
}
ost << "end" << endl;
}
void CSGeometry :: Load (istream & ist)
{
// CSGeometry * geo = new CSGeometry;
char key[100], name[100], classname[100], sname[100];
int ncoeff, i, j;
2009-01-25 17:35:25 +05:00
Array<double> coeff;
2009-01-13 04:40:13 +05:00
while (ist.good())
{
ist >> key;
if (strcmp (key, "boundingbox") == 0)
{
Point<3> pmin, pmax;
ist >> pmin(0) >> pmin(1) >> pmin(2);
ist >> pmax(0) >> pmax(1) >> pmax(2);
SetBoundingBox (Box<3> (pmin, pmax));
}
if (strcmp (key, "primitive") == 0)
{
ist >> name >> classname >> ncoeff;
coeff.SetSize (ncoeff);
for (i = 0; i < ncoeff; i++)
ist >> coeff[i];
Primitive * nprim = Primitive::CreatePrimitive (classname);
nprim -> SetPrimitiveData (coeff);
Solid * nsol = new Solid (nprim);
for (j = 0; j < nprim->GetNSurfaces(); j++)
{
sprintf (sname, "%s,%d", name, j);
AddSurface (sname, &nprim->GetSurface(j));
nprim -> SetSurfaceId (j, GetNSurf());
}
SetSolid (name, nsol);
}
else if (strcmp (key, "solid") == 0)
{
ist >> name;
Solid * nsol = Solid::CreateSolid (ist, solids);
cout << " I have found solid " << name << " = ";
nsol -> GetSolidData (cout);
cout << endl;
SetSolid (name, nsol);
}
else if (strcmp (key, "toplevel") == 0)
{
char type[20], solname[50], surfname[50];
const Solid * sol = NULL;
const Surface * surf = NULL;
int nr;
ist >> type;
if (strcmp (type, "solid") == 0)
{
ist >> solname;
sol = GetSolid (solname);
}
if (strcmp (type, "surface") == 0)
{
ist >> solname >> surfname;
sol = GetSolid (solname);
surf = GetSurface (surfname);
}
nr = SetTopLevelObject ((Solid*)sol, (Surface*)surf);
GetTopLevelObject (nr) -> SetData (ist);
}
else if (strcmp (key, "identify") == 0)
{
char type[10], surfname1[50], surfname2[50];
const Surface * surf1;
const Surface * surf2;
ist >> type >> surfname1 >> surfname2;
surf1 = GetSurface(surfname1);
surf2 = GetSurface(surfname2);
AddIdentification (new PeriodicIdentification
(GetNIdentifications(),
*this, surf1, surf2));
}
else if (strcmp (key, "end") == 0)
break;
}
changeval++;
}
2011-01-11 01:18:01 +05:00
void CSGeometry :: SaveSurfaces (ostream & out) const
2009-01-13 04:40:13 +05:00
{
if(singfaces.Size() > 0 || singedges.Size() > 0 || singpoints.Size() > 0)
{
PrintMessage(3,"Singular faces/edges/points => no csg-information in .vol file");
return;
}
2009-01-25 17:35:25 +05:00
Array<double> coeffs;
2009-01-13 04:40:13 +05:00
const char * classname;
out << "csgsurfaces " << GetNSurf() << "\n";
for(int i=0; i<GetNSurf(); i++)
{
const OneSurfacePrimitive * sp = dynamic_cast< const OneSurfacePrimitive * > (GetSurface(i));
const ExtrusionFace * ef = dynamic_cast< const ExtrusionFace * > (GetSurface(i));
const RevolutionFace * rf = dynamic_cast< const RevolutionFace * > (GetSurface(i));
2009-02-25 21:06:34 +05:00
const DummySurface * dummyf = dynamic_cast< const DummySurface * > (GetSurface(i));
2009-01-13 04:40:13 +05:00
if(sp)
{
sp->GetPrimitiveData(classname,coeffs);
out << classname << " ";
}
else if(ef)
{
out << "extrusionface ";
ef->GetRawData(coeffs);
}
else if(rf)
{
out << "revolutionface ";
rf->GetRawData(coeffs);
}
2009-02-25 21:06:34 +05:00
else if(dummyf)
{
out << "dummy ";
coeffs.SetSize(0);
}
2009-01-13 04:40:13 +05:00
else
throw NgException ("Cannot write csg surface. Please, contact developers!");
out << coeffs.Size() << "\n";
for(int j=0; j<coeffs.Size(); j++)
out << coeffs[j] << " ";
out << "\n";
}
}
void CSGeometry :: LoadSurfaces (istream & in)
{
2009-01-25 17:35:25 +05:00
Array<double> coeffs;
2009-01-13 04:40:13 +05:00
string classname;
int nsurfaces,size;
in >> classname;
if(classname == "csgsurfaces")
in >> nsurfaces;
else
nsurfaces = atoi(classname.c_str());
Point<3> dummypoint(0,0,0);
Vec<3> dummyvec(0,0,0);
double dummydouble(0.1);
for(int i=0; i<nsurfaces; i++)
{
in >> classname;
in >> size;
coeffs.SetSize(size);
for(int j=0; j<size; j++)
in >> coeffs[j];
if(classname == "plane")
{
Plane * plane = new Plane(dummypoint,dummyvec);
plane->SetPrimitiveData(coeffs);
AddSurface(plane);
delete_them.Append(plane);
}
else if(classname == "sphere")
{
Sphere * sphere = new Sphere(dummypoint,dummydouble);
sphere->SetPrimitiveData(coeffs);
AddSurface(sphere);
delete_them.Append(sphere);
}
else if(classname == "cylinder")
{
Cylinder * cylinder = new Cylinder(coeffs);
AddSurface(cylinder);
delete_them.Append(cylinder);
}
2009-10-13 03:41:54 +06:00
else if(classname == "ellipticcylinder")
{
EllipticCylinder * cylinder = new EllipticCylinder(coeffs);
AddSurface(cylinder);
delete_them.Append(cylinder);
}
2011-04-29 00:33:19 +06:00
else if(classname == "torus")
{
Torus * torus = new Torus(dummypoint,dummyvec,dummydouble, dummydouble);
torus->SetPrimitiveData(coeffs);
AddSurface(torus);
delete_them.Append(torus);
}
2009-10-13 03:41:54 +06:00
2009-01-13 04:40:13 +05:00
else if(classname == "cone")
{
Cone * cone = new Cone(dummypoint,dummypoint,dummydouble,dummydouble);
cone->SetPrimitiveData(coeffs);
AddSurface(cone);
delete_them.Append(cone);
}
else if(classname == "extrusionface")
{
ExtrusionFace * ef =
new ExtrusionFace(coeffs);
AddSurface(ef);
delete_them.Append(ef);
}
else if(classname == "revolutionface")
{
RevolutionFace * rf =
new RevolutionFace(coeffs);
AddSurface(rf);
delete_them.Append(rf);
}
2009-02-25 21:06:34 +05:00
else if(classname == "dummy")
{
Surface * surf = new DummySurface();
AddSurface(surf);
delete_them.Append(surf);
}
2009-01-13 04:40:13 +05:00
}
}
2011-01-11 01:18:01 +05:00
void CSGeometry :: SaveToMeshFile (ostream & ost) const
{
SaveSurfaces (ost);
}
2009-01-13 04:40:13 +05:00
void CSGeometry :: AddSurface (Surface * surf)
{
static int cntsurfs = 0;
cntsurfs++;
char name[15];
sprintf (name, "nnsurf%d", cntsurfs);
AddSurface (name, surf);
}
void CSGeometry :: AddSurface (char * name, Surface * surf)
{
2014-12-05 20:58:49 +05:00
(*testout) << "Adding surface " << name << ": " << (*surf) << endl;
2009-01-13 04:40:13 +05:00
surfaces.Set (name, surf);
surf->SetName (name);
changeval++;
}
void CSGeometry :: AddSurfaces (Primitive * prim)
{
for (int i = 0; i < prim->GetNSurfaces(); i++)
{
AddSurface (&prim->GetSurface(i));
prim->SetSurfaceId (i, GetNSurf()-1);
surf2prim.Append (prim);
}
}
const Surface * CSGeometry :: GetSurface (const char * name) const
{
if (surfaces.Used(name))
return surfaces.Get(name);
else
return NULL;
}
/*
const Surface * CSGeometry :: GetSurface (int i) const
{
if (i >= 0 && i < surfaces.Size())
return surfaces[i];
else
throw NgException ("CSGeometry::GetSurface out of range");
}
*/
void CSGeometry :: SetSolid (const char * name, Solid * sol)
{
Solid * oldsol = NULL;
if (solids.Used (name))
oldsol = solids.Get(name);
solids.Set (name, sol);
sol->SetName (name);
if (oldsol)
{
if (oldsol->op != Solid::ROOT ||
sol->op != Solid::ROOT)
{
cerr << "Setsolid: old or new no root" << endl;
}
oldsol -> s1 = sol -> s1;
}
changeval++;
}
const Solid * CSGeometry :: GetSolid (const char * name) const
{
if (solids.Used(name))
return solids.Get(name);
else
return NULL;
}
const Solid * CSGeometry :: GetSolid (const string & name) const
{
if (solids.Used(name.c_str()))
return solids.Get(name.c_str());
else
return NULL;
}
void CSGeometry :: SetSplineCurve (const char * name, SplineGeometry<2> * spl)
{
splinecurves2d.Set(name,spl);
}
void CSGeometry :: SetSplineCurve (const char * name, SplineGeometry<3> * spl)
{
splinecurves3d.Set(name,spl);
}
const SplineGeometry<2> * CSGeometry :: GetSplineCurve2d (const string & name) const
{
if (splinecurves2d.Used(name.c_str()))
return splinecurves2d.Get(name.c_str());
else
return NULL;
}
const SplineGeometry<3> * CSGeometry :: GetSplineCurve3d (const string & name) const
{
if (splinecurves3d.Used(name.c_str()))
return splinecurves3d.Get(name.c_str());
else
return NULL;
}
2009-04-27 23:44:28 +06:00
/*
2009-01-13 04:40:13 +05:00
class RemoveDummyIterator : public SolidIterator
{
public:
RemoveDummyIterator() { ; }
virtual ~RemoveDummyIterator() { ; }
virtual void Do(Solid * sol);
};
void RemoveDummyIterator :: Do(Solid * sol)
{
2009-04-27 23:44:28 +06:00
cerr << "remove dummy iterator is obsolete" << endl;
2009-01-13 04:40:13 +05:00
if ( (sol->op == Solid::SUB || sol->op == Solid::SECTION ||
sol->op == Solid::UNION)
&& sol->s1->op == Solid::DUMMY)
sol->s1 = sol->s1->s1;
if ( (sol->op == Solid::SECTION || sol->op == Solid::UNION)
&& sol->s2->op == Solid::DUMMY)
sol->s2 = sol->s2->s1;
}
2009-04-27 23:44:28 +06:00
*/
2009-01-13 04:40:13 +05:00
int CSGeometry :: SetTopLevelObject (Solid * sol, Surface * surf)
{
return toplevelobjects.Append (new TopLevelObject (sol, surf)) - 1;
}
TopLevelObject * CSGeometry ::
GetTopLevelObject (const Solid * sol, const Surface * surf)
{
for (int i = 0; i < toplevelobjects.Size(); i++)
{
if (toplevelobjects[i]->GetSolid() == sol &&
toplevelobjects[i]->GetSurface() == surf)
return (toplevelobjects[i]);
}
return NULL;
}
void CSGeometry :: RemoveTopLevelObject (Solid * sol, Surface * surf)
{
for (int i = 0; i < toplevelobjects.Size(); i++)
{
if (toplevelobjects[i]->GetSolid() == sol &&
toplevelobjects[i]->GetSurface() == surf)
{
delete toplevelobjects[i];
toplevelobjects.DeleteElement (i+1);
changeval++;
break;
}
}
}
void CSGeometry :: AddIdentification (Identification * ident)
{
identifications.Append (ident);
}
void CSGeometry :: SetFlags (const char * solidname, const Flags & flags)
{
Solid * solid = solids.Elem(solidname);
2009-01-25 17:35:25 +05:00
Array<int> surfind;
2009-01-13 04:40:13 +05:00
int i;
double maxh = flags.GetNumFlag ("maxh", -1);
if (maxh > 0 && solid)
{
solid->GetSurfaceIndices (surfind);
for (i = 0; i < surfind.Size(); i++)
{
if (surfaces[surfind[i]]->GetMaxH() > maxh)
surfaces[surfind[i]] -> SetMaxH (maxh);
}
solid->SetMaxH (maxh);
}
if ( flags.StringFlagDefined ("bcname") )
{
solid->GetSurfaceIndices (surfind);
string bcn = flags.GetStringFlag("bcname", "default");
for (i = 0; i < surfind.Size(); i++)
{
if(surfaces[surfind[i]]->GetBCName() == "default")
surfaces[surfind[i]]->SetBCName(bcn);
}
}
if (flags.StringListFlagDefined ("bcname"))
{
2009-01-25 17:35:25 +05:00
const Array<char*> & bcname = flags.GetStringListFlag("bcname");
2009-01-13 04:40:13 +05:00
Polyhedra * polyh;
if(solid->S1())
polyh = dynamic_cast<Polyhedra *>(solid->S1()->GetPrimitive());
else
polyh = dynamic_cast<Polyhedra *>(solid->GetPrimitive());
if(polyh)
{
2009-01-25 17:35:25 +05:00
Array < Array<int> * > polysurfs;
2009-01-13 04:40:13 +05:00
polyh->GetPolySurfs(polysurfs);
if(bcname.Size() != polysurfs.Size())
cerr << "WARNING: solid \"" << solidname << "\" has " << polysurfs.Size()
<< " surfaces and should get " << bcname.Size() << " bc-names!" << endl;
for ( i = 0; i < min2(polysurfs.Size(),bcname.Size()); i++)
{
for (int j = 0; j < polysurfs[i]->Size(); j++)
{
if(surfaces[(*polysurfs[i])[j]]->GetBCName() == "default")
surfaces[(*polysurfs[i])[j]]->SetBCName(bcname[i]);
}
delete polysurfs[i];
}
}
else
{
solid->GetSurfaceIndices (surfind);
if(bcname.Size() != surfind.Size())
cerr << "WARNING: solid \"" << solidname << "\" has " << surfind.Size()
<< " surfaces and should get " << bcname.Size() << " bc-names!" << endl;
for (i = 0; i < min2(surfind.Size(),bcname.Size()); i++)
{
if(surfaces[surfind[i]]->GetBCName() == "default")
surfaces[surfind[i]]->SetBCName(bcname[i]);
}
}
}
if (flags.NumFlagDefined ("bc"))
{
solid->GetSurfaceIndices (surfind);
int bc = int (flags.GetNumFlag("bc", -1));
for (i = 0; i < surfind.Size(); i++)
{
if (surfaces[surfind[i]]->GetBCProperty() == -1)
surfaces[surfind[i]]->SetBCProperty(bc);
}
}
if (flags.NumListFlagDefined ("bc"))
{
2009-01-25 17:35:25 +05:00
const Array<double> & bcnum = flags.GetNumListFlag("bc");
2009-01-13 04:40:13 +05:00
Polyhedra * polyh;
if(solid->S1())
polyh = dynamic_cast<Polyhedra *>(solid->S1()->GetPrimitive());
else
polyh = dynamic_cast<Polyhedra *>(solid->GetPrimitive());
if(polyh)
{
2009-01-25 17:35:25 +05:00
Array < Array<int> * > polysurfs;
2009-01-13 04:40:13 +05:00
polyh->GetPolySurfs(polysurfs);
if(bcnum.Size() != polysurfs.Size())
cerr << "WARNING: solid \"" << solidname << "\" has " << polysurfs.Size()
<< " surfaces and should get " << bcnum.Size() << " bc-numbers!" << endl;
for ( i = 0; i < min2(polysurfs.Size(),bcnum.Size()); i++)
{
for (int j = 0; j < polysurfs[i]->Size(); j++)
{
if ( surfaces[(*polysurfs[i])[j]]->GetBCProperty() == -1 )
surfaces[(*polysurfs[i])[j]]->SetBCProperty(int(bcnum[i]));
}
delete polysurfs[i];
}
}
else
{
solid->GetSurfaceIndices (surfind);
if(bcnum.Size() != surfind.Size())
cerr << "WARNING: solid \"" << solidname << "\" has " << surfind.Size()
<< " surfaces and should get " << bcnum.Size() << " bc-numbers!" << endl;
for (i = 0; i < min2(surfind.Size(),bcnum.Size()); i++)
{
if (surfaces[surfind[i]]->GetBCProperty() == -1)
surfaces[surfind[i]]->SetBCProperty(int(bcnum[i]));
}
}
}
}
void CSGeometry :: FindIdenticSurfaces (double eps)
{
int inv;
int nsurf = GetNSurf();
isidenticto.SetSize(nsurf);
for (int i = 0; i < nsurf; i++)
isidenticto[i] = i;
for (int i = 0; i < nsurf; i++)
for (int j = i+1; j < nsurf; j++)
{
if (GetSurface(j) -> IsIdentic (*GetSurface(i), inv, eps))
{
INDEX_2 i2(i, j);
identicsurfaces.Set (i2, inv);
isidenticto[j] = isidenticto[i];
}
}
(*testout) << "identicmap:" << endl;
for (int i = 0; i < isidenticto.Size(); i++)
(*testout) << i << " -> " << isidenticto[i] << endl;
}
void CSGeometry ::
GetSurfaceIndices (const Solid * sol,
const BoxSphere<3> & box,
2009-01-25 17:35:25 +05:00
Array<int> & locsurf) const
2009-01-13 04:40:13 +05:00
{
ReducePrimitiveIterator rpi(box);
UnReducePrimitiveIterator urpi;
2012-11-05 17:28:36 +06:00
const_cast<Solid*> (sol) -> IterateSolid (rpi);
2009-01-13 04:40:13 +05:00
sol -> GetSurfaceIndices (locsurf);
2012-11-05 17:28:36 +06:00
const_cast<Solid*> (sol) -> IterateSolid (urpi);
2009-01-13 04:40:13 +05:00
for (int i = locsurf.Size()-1; i >= 0; i--)
{
bool indep = 1;
for (int j = 0; j < i; j++)
if (locsurf[i] == locsurf[j])
{
indep = 0;
break;
}
if (!indep) locsurf.Delete(i);
}
}
void CSGeometry ::
GetIndependentSurfaceIndices (const Solid * sol,
const BoxSphere<3> & box,
2009-01-25 17:35:25 +05:00
Array<int> & locsurf) const
2009-01-13 04:40:13 +05:00
{
ReducePrimitiveIterator rpi(box);
UnReducePrimitiveIterator urpi;
((Solid*)sol) -> IterateSolid (rpi);
sol -> GetSurfaceIndices (locsurf);
((Solid*)sol) -> IterateSolid (urpi);
for (int i = 0; i < locsurf.Size(); i++)
locsurf[i] = isidenticto[locsurf[i]];
for (int i = locsurf.Size()-1; i >= 0; i--)
{
bool indep = 1;
for (int j = 0; j < i; j++)
if (locsurf[i] == locsurf[j])
{
indep = 0;
break;
}
if (!indep) locsurf.Delete(i);
}
/*
// delete identified
for (int i = locsurf.Size()-1; i >= 0; i--)
{
bool indep = 1;
for (int j = 0; j < i; j++)
{
if (identicsurfaces.Used (INDEX_2::Sort (locsurf[i], locsurf[j])) !=
(isidenticto[locsurf[i]] == isidenticto[locsurf[j]]))
{
cerr << "different result" << endl;
exit(1);
}
if (isidenticto[locsurf[i]] == isidenticto[locsurf[j]])
{
indep = 0;
break;
}
}
if (!indep)
locsurf.Delete(i);
}
for (int i = 0; i < locsurf.Size(); i++)
locsurf[i] = isidenticto[locsurf[i]];
*/
}
2012-11-05 17:28:36 +06:00
/*
2009-01-13 04:40:13 +05:00
void CSGeometry ::
GetIndependentSurfaceIndices (const Solid * sol,
const Point<3> & p, Vec<3> & v,
2009-01-25 17:35:25 +05:00
Array<int> & locsurf) const
2009-01-13 04:40:13 +05:00
{
cout << "very dangerous" << endl;
Point<3> p2 = p + 1e-2 * v;
BoxSphere<3> box (p2, p2);
box.Increase (1e-3);
box.CalcDiamCenter();
GetIndependentSurfaceIndices (sol, box, locsurf);
}
2012-11-05 17:28:36 +06:00
*/
2009-01-13 04:40:13 +05:00
void CSGeometry ::
2009-01-25 17:35:25 +05:00
GetIndependentSurfaceIndices (Array<int> & locsurf) const
2009-01-13 04:40:13 +05:00
{
for (int i = 0; i < locsurf.Size(); i++)
locsurf[i] = isidenticto[locsurf[i]];
for (int i = locsurf.Size()-1; i >= 0; i--)
{
bool indep = 1;
for (int j = 0; j < i; j++)
if (locsurf[i] == locsurf[j])
{
indep = 0;
break;
}
if (!indep) locsurf.Delete(i);
}
}
void CSGeometry ::
2011-01-11 01:18:01 +05:00
CalcTriangleApproximation(double detail, double facets)
2009-01-13 04:40:13 +05:00
{
PrintMessage (1, "Calc Triangle Approximation");
try
{
2009-01-13 04:40:13 +05:00
// FindIdenticSurfaces (1e-6);
int ntlo = GetNTopLevelObjects();
for (int i = 0; i < triapprox.Size(); i++)
delete triapprox[i];
triapprox.SetSize (ntlo);
2009-01-25 17:35:25 +05:00
Array<int> surfind;
2009-01-13 04:40:13 +05:00
IndexSet iset(GetNSurf());
for (int i = 0; i < ntlo; i++)
{
Solid * sol;
Surface * surf;
GetTopLevelObject (i, sol, surf);
sol -> CalcSurfaceInverse ();
TriangleApproximation * tams = new TriangleApproximation();
triapprox[i] = tams;
// sol -> GetSurfaceIndices (surfind);
for (int j = 0; j < GetNSurf(); j++)
// for (int jj = 0; jj < surfind.Size(); jj++)
{
// int j = surfind[jj];
PrintMessageCR (3, "Surface ", j, "/", GetNSurf());
// PrintMessageCR (3, "Surface ", j, "/", surfind.Size());
if (surf && surf != GetSurface(j))
continue;
TriangleApproximation tas;
2011-01-11 01:18:01 +05:00
GetSurface (j) -> GetTriangleApproximation (tas, boundingbox, facets);
2009-01-13 04:40:13 +05:00
int oldnp = tams -> GetNP();
if (!tas.GetNP())
continue;
for (int k = 0; k < tas.GetNP(); k++)
{
tams -> AddPoint (tas.GetPoint(k));
2009-04-03 20:39:52 +06:00
Vec<3> n = GetSurface(j) -> GetNormalVector (tas.GetPoint(k));
2009-01-13 04:40:13 +05:00
n.Normalize();
if (GetSurface(j)->Inverse()) n *= -1;
tams -> AddNormal (n);
}
BoxSphere<3> surfbox;
if (tas.GetNP())
surfbox.Set (tas.GetPoint(0));
for (int k = 1; k < tas.GetNP(); k++)
surfbox.Add (tas.GetPoint(k));
surfbox.Increase (1e-6);
surfbox.CalcDiamCenter();
Solid * surflocsol = sol -> GetReducedSolid (surfbox);
if (!surflocsol)
continue;
for (int k = 0; k < tas.GetNT(); k++)
{
const TATriangle & tri = tas.GetTriangle (k);
// check triangle
BoxSphere<3> box;
box.Set (tas.GetPoint (tri[0]));
box.Add (tas.GetPoint (tri[1]));
box.Add (tas.GetPoint (tri[2]));
box.Increase (1e-6);
box.CalcDiamCenter();
Solid * locsol = surflocsol -> GetReducedSolid (box);
if (locsol)
{
TATriangle tria(j,
tri[0] + oldnp,
tri[1] + oldnp,
tri[2] + oldnp);
2009-04-03 20:39:52 +06:00
// tams -> AddTriangle (tria);
2009-01-13 04:40:13 +05:00
RefineTriangleApprox (locsol, j, box, detail,
2009-04-03 20:39:52 +06:00
tria, *tams, iset, 1);
2009-01-13 04:40:13 +05:00
delete locsol;
}
}
}
tams->RemoveUnusedPoints ();
PrintMessage (2, "Object ", i, " has ", tams->GetNT(), " triangles");
}
}
catch (exception)
{
cerr << "*************************************************************" << endl
<< "**** out of memory problem in CSG visualization ****" << endl
<< "**** Restart netgen, and disable ****" << endl
<< "**** 'Draw Geometry' in Geometry -> CSG Options ****" << endl
<< "**** before loading the geometry ****" << endl
<< "**** meshing will still work ! ****" << endl
<< "*************************************************************" << endl;
exit(1);
}
2009-01-13 04:40:13 +05:00
Change();
}
void CSGeometry ::
RefineTriangleApprox (Solid * locsol,
int surfind,
const BoxSphere<3> & box,
double detail,
const TATriangle & tria,
TriangleApproximation & tams,
2009-04-03 20:39:52 +06:00
IndexSet & iset,
int level)
2009-01-13 04:40:13 +05:00
{
2009-04-05 17:39:00 +06:00
// if (level > 10) return;
2009-04-03 20:39:52 +06:00
2009-01-13 04:40:13 +05:00
//tams.AddTriangle (tria);
//(*testout) << "tria " << tams.GetPoint(tria[0]) << " - " << tams.GetPoint(tria[1]) << " - " << tams.GetPoint(tria[2])
// << " ( " << tria[0] << " " << tria[1] << " " << tria[2] << ")" <<endl;
//return;
int pinds[6];
ArrayMem<int,500> surfused(GetNSurf());
ReducePrimitiveIterator rpi(box);
UnReducePrimitiveIterator urpi;
locsol -> IterateSolid (rpi);
// locsol -> GetSurfaceIndices (lsurfi);
// IndexSet iset(GetNSurf());
locsol -> GetSurfaceIndices (iset);
2009-01-25 20:04:33 +05:00
const Array<int> & lsurfi = iset.GetArray();
2009-01-13 04:40:13 +05:00
locsol -> IterateSolid (urpi);
int surfii = -1;
for (int i = 0; i < lsurfi.Size(); i++)
if (lsurfi[i] == surfind)
{
surfii = i;
break;
}
if (surfii == -1)
return;
int cntindep = 0;
for (int i = 0; i < lsurfi.Size(); i++)
{
int linkto = isidenticto[lsurfi[i]];
surfused[linkto] = 0;
}
for (int i = 0; i < lsurfi.Size(); i++)
{
int linkto = isidenticto[lsurfi[i]];
if (!surfused[linkto])
{
surfused[linkto] = 1;
cntindep++;
}
}
int inverse = surfaces[surfind]->Inverse();
if (cntindep == 1)
{
tams.AddTriangle (tria);
//(*testout) << "pos1 " << tams.GetPoint(tria[0]) << " - " << tams.GetPoint(tria[1]) << " - " << tams.GetPoint(tria[2]) << endl;
return;
}
if (cntindep == 2)
{
// just 2 surfaces:
// if smooth, project inner points to edge and finish
int otherind = -1;
for (int i = 0; i < lsurfi.Size(); i++)
{
INDEX_2 i2 (lsurfi[i], surfind);
i2.Sort();
if (i != surfii && !identicsurfaces.Used(i2))
otherind = lsurfi[i];
}
double kappa = GetSurface(otherind)-> MaxCurvature ();
if (kappa * box.Diam() < 0.1)
{
int pnums[6];
static int between[3][3] =
{ { 1, 2, 3 },
{ 0, 2, 4 },
{ 0, 1, 5 } };
int onsurface[3];
for (int j = 0; j < 3; j++)
{
int pi = tria[j];
pnums[j] = pi;
onsurface[j] =
!locsol->IsStrictIn (tams.GetPoint (pi), 1e-6) &&
locsol->IsIn (tams.GetPoint (pi), 1e-6);
//
/*
static int nos=0;
if(!onsurface[j])
{
nos++;
cout << "NOT ON SURFACE!! "<< nos << endl;
}
*/
}
for (int j = 0; j < 3; j++)
{
int lpi1 = between[j][0];
int lpi2 = between[j][1];
int lpin = between[j][2];
if (onsurface[lpi1] == onsurface[lpi2])
pnums[lpin] = -1;
else
{
const Point<3> & p1 = tams.GetPoint (pnums[lpi1]);
const Point<3> & p2 = tams.GetPoint (pnums[lpi2]);
double f1 = GetSurface(otherind)->CalcFunctionValue (p1);
double f2 = GetSurface(otherind)->CalcFunctionValue (p2);
Point<3> pn;
double l2(100),l1(100);
if ( fabs (f1-f2) > 1e-20 )
{
l2 = -f1/(f2-f1);
l1 = f2/(f2-f1);
pn = Point<3>(l1 * p1(0) + l2 * p2(0),
l1 * p1(1) + l2 * p2(1),
l1 * p1(2) + l2 * p2(2));
}
else
pn = p1;
// if(fabs(pn(0)) > 4 || fabs(pn(1)) > 4 || fabs(pn(2)) > 4)
// {
// cout << "p1 " << p1 << " p2 " << p2
// << " f1 " << f1 << " f2 " << f2
// << " l1 " << l1 << " l2 " << l2
// << " pn " << pn << endl;
// }
//GetSurface (surfind)->Project (pn);
pnums[lpin] = tams.AddPoint (pn);
GetSurface (surfind)->Project (pn);
Vec<3> n;
n = GetSurface (surfind)->GetNormalVector (pn);
if (inverse) n *= -1;
tams.AddNormal(n);
}
}
int vcase = 0;
if (onsurface[0]) vcase++;
if (onsurface[1]) vcase+=2;
if (onsurface[2]) vcase+=4;
static int trias[8][6] =
{ { 0, 0, 0, 0, 0, 0 },
{ 1, 6, 5, 0, 0, 0 },
{ 2, 4, 6, 0, 0, 0 },
{ 1, 2, 4, 1, 4, 5 },
{ 3, 5, 4, 0, 0, 0 },
{ 1, 6, 4, 1, 4, 3 },
{ 2, 3, 6, 3, 5, 6 },
{ 1, 2, 3, 0, 0, 0 } };
static int ntrias[4] =
{ 0, 1, 2, 1 };
int nvis = 0;
for (int j = 0; j < 3; j++)
if (onsurface[j])
nvis++;
for (int j = 0; j < ntrias[nvis]; j++)
{
TATriangle ntria(tria.SurfaceIndex(),
pnums[trias[vcase][3*j]-1],
pnums[trias[vcase][3*j+1]-1],
pnums[trias[vcase][3*j+2]-1]);
//(*testout) << "pos2 " << tams.GetPoint(ntria[0]) << " - " << tams.GetPoint(ntria[1]) << " - " << tams.GetPoint(ntria[2]) << endl
// << "( " << ntria[0] << " - " << ntria[1] << " - " << ntria[2] << ")" << endl;
tams.AddTriangle (ntria);
}
/* saturn changes:
int pvis[3];
for (j = 0; j < 3; j++)
pvis[j] = !locsol->IsStrictIn (tams.GetPoint (j+1), 1e-6) &&
locsol->IsIn (tams.GetPoint (j+1), 1e-6);
int newpi[3];
for (j = 0; j < 3; j++)
{
int pi1 = j;
int pi2 = (j+1) % 3;
int pic = j;
if (pvis[pi1] != pvis[pi2])
{
Point<3> hp = Center (tams.GetPoint (tria.PNum (pi1+1)),
tams.GetPoint (tria.PNum (pi2+1)));
newpi[j] = tams.AddPoint (hp);
Vec<3> n = tams.GetNormal (pi1);
tams.AddNormal (n);
}
else
newpi[j] = 0;
}
int nvis = 0;
for (j = 0; j <= nvis; j++)
if (pvis[j]) nvis++;
int si = tria.SurfaceIndex();
switch (nvis)
{
case 0:
break;
case 1:
{
int visj;
for (j = 0; j < 3; j++)
if (pvis[j]) visj = j;
int pivis = tria.PNum (visj+1);
int pic1 = newpi[(visj+1)%3];
int pic2 = newpi[(visj+2)%3];
cout << pivis << "," << pic1 << "," << pic2 << endl;
tams.AddTriangle (TATriangle (si, pivis, pic1,pic2));
break;
}
case 2:
{
int nvisj;
for (j = 0; j < 3; j++)
if (!pvis[j]) nvisj = j;
int pivis1 = tria.PNum ((nvisj+1)%3+1);
int pivis2 = tria.PNum ((nvisj+2)%3+1);
int pic1 = newpi[nvisj];
int pic2 = newpi[(nvisj+2)%3];
tams.AddTriangle (TATriangle (si, pivis1, pic1,pic2));
tams.AddTriangle (TATriangle (si, pivis1, pic1,pivis2));
break;
}
case 3:
{
tams.AddTriangle (tria);
break;
}
}
*/
return;
}
}
// bisection
if (box.Diam() < detail)
{
//cout << "returning" << endl;
return;
}
for (int i = 0; i < 3; i++)
pinds[i] = tria[i];
static int between[3][3] =
{ { 0, 1, 5 },
{ 0, 2, 4 },
{ 1, 2, 3 } };
for (int i = 0; i < 3; i++)
{
// int pi1 = tria[between[i][0]];
Point<3> newp = Center (tams.GetPoint (tria[between[i][0]]),
tams.GetPoint (tria[between[i][1]]));
Vec<3> n;
GetSurface(surfind)->Project (newp);
2009-04-03 20:39:52 +06:00
2009-01-13 04:40:13 +05:00
n = GetSurface(surfind)->GetNormalVector (newp);
pinds[between[i][2]] = tams.AddPoint (newp);
if (inverse) n *= -1;
tams.AddNormal (n);
}
static int trias[4][4] =
{ { 0, 5, 4 },
{ 5, 1, 3 },
{ 4, 3, 2 },
{ 3, 4, 5 } };
for (int i = 0; i < 4; i++)
{
TATriangle ntri(surfind,
pinds[trias[i][0]],
pinds[trias[i][1]],
pinds[trias[i][2]]);
// check triangle
BoxSphere<3> nbox;
nbox.Set (tams.GetPoint (ntri[0]));
nbox.Add (tams.GetPoint (ntri[1]));
nbox.Add (tams.GetPoint (ntri[2]));
2009-04-03 20:39:52 +06:00
nbox.Increase (1e-8);
2009-01-13 04:40:13 +05:00
nbox.CalcDiamCenter();
Solid * nsol = locsol -> GetReducedSolid (nbox);
if (nsol)
{
RefineTriangleApprox (nsol, surfind, nbox,
2009-04-03 20:39:52 +06:00
detail, ntri, tams, iset, level+1);
2009-01-13 04:40:13 +05:00
delete nsol;
}
}
}
class ClearVisitedIt : public SolidIterator
{
public:
ClearVisitedIt () { ; }
virtual ~ClearVisitedIt () { ; }
virtual void Do (Solid * sol)
{
sol -> visited = 0;
}
};
void CSGeometry ::
2011-01-11 01:18:01 +05:00
IterateAllSolids (SolidIterator & it, bool only_once) const
2009-01-13 04:40:13 +05:00
{
if (only_once)
{
ClearVisitedIt clit;
for (int i = 0; i < solids.Size(); i++)
solids[i] -> IterateSolid (clit, 0);
}
for (int i = 0; i < solids.Size(); i++)
solids[i] -> IterateSolid (it, only_once);
}
double CSGeometry :: MaxSize () const
{
double maxs, mins;
maxs = max3 (boundingbox.PMax()(0),
boundingbox.PMax()(1),
boundingbox.PMax()(2));
mins = min3 (boundingbox.PMin()(0),
boundingbox.PMin()(1),
boundingbox.PMin()(2));
return max2 (maxs, -mins) * 1.1;
}
2012-08-20 20:10:23 +06:00
class CSGeometryRegister : public GeometryRegister
{
public:
virtual NetgenGeometry * Load (string filename) const;
virtual NetgenGeometry * LoadFromMeshFile (istream & ist) const;
// virtual VisualScene * GetVisualScene (const NetgenGeometry * geom) const;
};
2016-11-04 16:14:52 +05:00
extern CSGeometry * ParseCSG (istream & istr, CSGeometry *instance=nullptr);
2012-08-20 20:10:23 +06:00
NetgenGeometry * CSGeometryRegister :: Load (string filename) const
{
const char * cfilename = filename.c_str();
if (strcmp (&cfilename[strlen(cfilename)-3], "geo") == 0)
{
PrintMessage (1, "Load CSG geometry file ", cfilename);
ifstream infile(cfilename);
CSGeometry * hgeom = ParseCSG (infile);
if (!hgeom)
throw NgException ("geo-file should start with 'algebraic3d'");
hgeom -> FindIdenticSurfaces(1e-8 * hgeom->MaxSize());
return hgeom;
}
if (strcmp (&cfilename[strlen(cfilename)-3], "ngg") == 0)
{
PrintMessage (1, "Load new CSG geometry file ", cfilename);
ifstream infile(cfilename);
CSGeometry * hgeom = new CSGeometry("");
hgeom -> Load (infile);
return hgeom;
}
return NULL;
}
NetgenGeometry * CSGeometryRegister :: LoadFromMeshFile (istream & ist) const
{
string auxstring;
if (ist.good())
{
ist >> auxstring;
if (auxstring == "csgsurfaces")
{
CSGeometry * geometry = new CSGeometry ("");
geometry -> LoadSurfaces(ist);
return geometry;
}
// else
// ist.putback (auxstring);
}
return NULL;
}
class CSGInit
{
public:
CSGInit()
{
geometryregister.Append (new CSGeometryRegister);
}
};
CSGInit csginit;
2009-01-13 04:40:13 +05:00
}