* Added OpenCascade XDE Support to enable importing of individual surface colours from STEP Geometry

* Extended the Clipping Planes functionality to the Geometry mode for OCC Geometry
* Added the option to specify the maximum mesh size for each individual face in an OCC Geometry
This commit is contained in:
Philippose Rajan 2009-01-30 22:17:20 +00:00
parent 34bfd4a349
commit 85867fb240
5 changed files with 2005 additions and 1796 deletions

View File

@ -564,8 +564,9 @@ namespace netgen
double maxh = mparam.maxh;
// Philippose - 15/01/2009
double maxh = geom.face_maxh[k-1];
//double maxh = mparam.maxh;
mparam.checkoverlap = 0;
// int noldpoints = mesh->GetNP();
int noldsurfel = mesh.GetNSE();
@ -1085,456 +1086,465 @@ namespace netgen
int OCCGenerateMesh (OCCGeometry & geom,
Mesh *& mesh,
int perfstepsstart, int perfstepsend,
char * optstr)
{
// int i, j;
Mesh *& mesh,
int perfstepsstart, int perfstepsend,
char * optstr)
{
// int i, j;
multithread.percent = 0;
multithread.percent = 0;
if (perfstepsstart <= MESHCONST_ANALYSE)
if (perfstepsstart <= MESHCONST_ANALYSE)
{
delete mesh;
mesh = new Mesh();
mesh->geomtype = Mesh::GEOM_OCC;
mesh->SetGlobalH (mparam.maxh);
mesh->SetMinimalH (mparam.minh);
delete mesh;
mesh = new Mesh();
mesh->geomtype = Mesh::GEOM_OCC;
mesh->SetGlobalH (mparam.maxh);
mesh->SetMinimalH (mparam.minh);
Array<double> maxhdom;
maxhdom.SetSize (geom.NrSolids());
maxhdom = mparam.maxh;
Array<double> maxhdom;
maxhdom.SetSize (geom.NrSolids());
maxhdom = mparam.maxh;
mesh->SetMaxHDomain (maxhdom);
mesh->SetMaxHDomain (maxhdom);
Box<3> bb = geom.GetBoundingBox();
bb.Increase (bb.Diam()/10);
Box<3> bb = geom.GetBoundingBox();
bb.Increase (bb.Diam()/10);
mesh->SetLocalH (bb.PMin(), bb.PMax(), 0.5);
mesh->SetLocalH (bb.PMin(), bb.PMax(), 0.5);
if (mparam.uselocalh)
{
if (mparam.uselocalh)
{
const char * savetask = multithread.task;
multithread.percent = 0;
const char * savetask = multithread.task;
multithread.percent = 0;
mesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading);
mesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading);
int nedges = geom.emap.Extent();
int nedges = geom.emap.Extent();
double maxedgelen = 0;
double minedgelen = 1e99;
double maxedgelen = 0;
double minedgelen = 1e99;
multithread.task = "Setting local mesh size (elements per edge)";
// setting elements per edge
multithread.task = "Setting local mesh size (elements per edge)";
for (int i = 1; i <= nedges && !multithread.terminate; i++)
{
TopoDS_Edge e = TopoDS::Edge (geom.emap(i));
multithread.percent = 100 * (i-1)/double(nedges);
if (BRep_Tool::Degenerated(e)) continue;
// setting elements per edge
GProp_GProps system;
BRepGProp::LinearProperties(e, system);
double len = system.Mass();
for (int i = 1; i <= nedges && !multithread.terminate; i++)
{
TopoDS_Edge e = TopoDS::Edge (geom.emap(i));
multithread.percent = 100 * (i-1)/double(nedges);
if (BRep_Tool::Degenerated(e)) continue;
if (len < IGNORECURVELENGTH)
{
(*testout) << "ignored" << endl;
continue;
}
GProp_GProps system;
BRepGProp::LinearProperties(e, system);
double len = system.Mass();
double localh = len/mparam.segmentsperedge;
double s0, s1;
if (len < IGNORECURVELENGTH)
{
(*testout) << "ignored" << endl;
continue;
}
// Philippose - 23/01/2009
// Find all the parent faces of a given edge
// and limit the mesh size of the edge based on the
// mesh size limit of the face
TopTools_IndexedDataMapOfShapeListOfShape edge_face_map;
edge_face_map.Clear();
TopExp::MapShapesAndAncestors(geom.shape, TopAbs_EDGE, TopAbs_FACE, edge_face_map);
const TopTools_ListOfShape& parent_faces = edge_face_map.FindFromKey(e);
double localh = len/mparam.segmentsperedge;
double s0, s1;
Handle(Geom_Curve) c = BRep_Tool::Curve(e, s0, s1);
TopTools_ListIteratorOfListOfShape parent_face_list;
maxedgelen = max (maxedgelen, len);
minedgelen = min (minedgelen, len);
for(parent_face_list.Initialize(parent_faces); parent_face_list.More(); parent_face_list.Next())
{
TopoDS_Face parent_face = TopoDS::Face(parent_face_list.Value());
int maxj = 2 * (int) ceil (localh/len);
for (int j = 0; j <= maxj; j++)
{
gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0));
mesh->RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh);
}
}
int face_index = geom.fmap.FindIndex(parent_face);
if(face_index >= 1) localh = min(localh,geom.face_maxh[face_index - 1]);
}
Handle(Geom_Curve) c = BRep_Tool::Curve(e, s0, s1);
maxedgelen = max (maxedgelen, len);
minedgelen = min (minedgelen, len);
multithread.task = "Setting local mesh size (edge curvature)";
// Philippose - 23/01/2009
// Modified the calculation of maxj, because the
// method used so far always results in maxj = 2,
// which causes the localh to be set only at the
// starting, mid and end of the edge.
// Old Algorithm:
// int maxj = 2 * (int) ceil (localh/len);
int maxj = max((int) ceil(len/localh), 2);
for (int j = 0; j <= maxj; j++)
{
gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0));
mesh->RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh);
}
}
// setting edge curvature
multithread.task = "Setting local mesh size (edge curvature)";
int nsections = 20;
// setting edge curvature
for (int i = 1; i <= nedges && !multithread.terminate; i++)
{
double maxcur = 0;
multithread.percent = 100 * (i-1)/double(nedges);
TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
if (BRep_Tool::Degenerated(edge)) continue;
double s0, s1;
Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);
BRepAdaptor_Curve brepc(edge);
BRepLProp_CLProps prop(brepc, 2, 1e-5);
int nsections = 20;
for (int j = 1; j <= nsections; j++)
{
double s = s0 + j/(double) nsections * (s1-s0);
prop.SetParameter (s);
double curvature = prop.Curvature();
if(curvature > maxcur) maxcur = curvature;
if (curvature >= 1e99)
continue;
gp_Pnt pnt = c->Value (s);
mesh->RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()),
ComputeH (fabs(curvature)));
}
// (*testout) << "edge " << i << " max. curvature: " << maxcur << endl;
}
multithread.task = "Setting local mesh size (face curvature)";
// setting face curvature
int nfaces = geom.fmap.Extent();
for (int i = 1; i <= nfaces && !multithread.terminate; i++)
{
multithread.percent = 100 * (i-1)/double(nfaces);
TopoDS_Face face = TopoDS::Face(geom.fmap(i));
TopLoc_Location loc;
Handle(Geom_Surface) surf = BRep_Tool::Surface (face);
Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);
if (triangulation.IsNull()) continue;
BRepAdaptor_Surface sf(face, Standard_True);
BRepLProp_SLProps prop(sf, 2, 1e-5);
int ntriangles = triangulation -> NbTriangles();
for (int j = 1; j <= ntriangles; j++)
{
gp_Pnt p[3];
gp_Pnt2d par[3];
for (int k = 1; k <=3; k++)
{
int n = triangulation->Triangles()(j)(k);
p[k-1] = triangulation->Nodes()(n).Transformed(loc);
par[k-1] = triangulation->UVNodes()(n);
}
//double maxside = 0;
//maxside = max (maxside, p[0].Distance(p[1]));
//maxside = max (maxside, p[0].Distance(p[2]));
//maxside = max (maxside, p[1].Distance(p[2]));
//cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush;
RestrictHTriangle (par[0], par[1], par[2], &prop, *mesh, 0);
//cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush;
}
}
// setting close edges
if (stlparam.resthcloseedgeenable)
{
multithread.task = "Setting local mesh size (close edges)";
int sections = 100;
Array<Line> lines(sections*nedges);
Box3dTree* searchtree =
new Box3dTree (bb.PMin(), bb.PMax());
int nlines = 0;
for (int i = 1; i <= nedges && !multithread.terminate; i++)
{
TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
if (BRep_Tool::Degenerated(edge)) continue;
double s0, s1;
Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);
BRepAdaptor_Curve brepc(edge);
BRepLProp_CLProps prop(brepc, 1, 1e-5);
prop.SetParameter (s0);
gp_Vec d0 = prop.D1().Normalized();
double s_start = s0;
int count = 0;
for (int j = 1; j <= sections; j++)
{
double s = s0 + (s1-s0)*(double)j/(double)sections;
prop.SetParameter (s);
gp_Vec d1 = prop.D1().Normalized();
double cosalpha = fabs(d0*d1);
if ((j == sections) || (cosalpha < cos(10.0/180.0*M_PI)))
{
count++;
gp_Pnt p0 = c->Value (s_start);
gp_Pnt p1 = c->Value (s);
lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z());
lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z());
Box3d box;
box.SetPoint (Point3d(lines[nlines].p0));
box.AddPoint (Point3d(lines[nlines].p1));
searchtree->Insert (box.PMin(), box.PMax(), nlines+1);
nlines++;
s_start = s;
d0 = d1;
}
}
}
Array<int> linenums;
for (int i = 0; i < nlines; i++)
{
multithread.percent = (100*i)/double(nlines);
Line & line = lines[i];
Box3d box;
box.SetPoint (Point3d(line.p0));
box.AddPoint (Point3d(line.p1));
double maxhline = max (mesh->GetH(box.PMin()),
mesh->GetH(box.PMax()));
box.Increase(maxhline);
double mindist = 1e99;
linenums.SetSize(0);
searchtree->GetIntersecting(box.PMin(),box.PMax(),linenums);
for (int j = 0; j < linenums.Size(); j++)
{
int num = linenums[j]-1;
if (i == num) continue;
if ((line.p0-lines[num].p0).Length2() < 1e-15) continue;
if ((line.p0-lines[num].p1).Length2() < 1e-15) continue;
if ((line.p1-lines[num].p0).Length2() < 1e-15) continue;
if ((line.p1-lines[num].p1).Length2() < 1e-15) continue;
mindist = min (mindist, line.Dist(lines[num]));
}
mindist *= stlparam.resthcloseedgefac;
if (mindist < 1e-3)
{
(*testout) << "extremely small local h: " << mindist
<< " --> setting to 1e-3" << endl;
(*testout) << "somewhere near " << line.p0 << " - " << line.p1 << endl;
mindist = 1e-3;
}
mesh->RestrictLocalHLine(line.p0, line.p1, mindist);
}
}
multithread.task = savetask;
}
for (int i = 1; i <= nedges && !multithread.terminate; i++)
{
double maxcur = 0;
multithread.percent = 100 * (i-1)/double(nedges);
TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
if (BRep_Tool::Degenerated(edge)) continue;
double s0, s1;
Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);
BRepAdaptor_Curve brepc(edge);
BRepLProp_CLProps prop(brepc, 2, 1e-5);
for (int j = 1; j <= nsections; j++)
{
double s = s0 + j/(double) nsections * (s1-s0);
prop.SetParameter (s);
double curvature = prop.Curvature();
if(curvature> maxcur) maxcur = curvature;
if (curvature >= 1e99)
continue;
gp_Pnt pnt = c->Value (s);
mesh->RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()),
ComputeH (fabs(curvature)));
}
// (*testout) << "edge " << i << " max. curvature: " << maxcur << endl;
}
multithread.task = "Setting local mesh size (face curvature)";
// setting face curvature
int nfaces = geom.fmap.Extent();
for (int i = 1; i <= nfaces && !multithread.terminate; i++)
{
multithread.percent = 100 * (i-1)/double(nfaces);
TopoDS_Face face = TopoDS::Face(geom.fmap(i));
TopLoc_Location loc;
Handle(Geom_Surface) surf = BRep_Tool::Surface (face);
Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);
if (triangulation.IsNull()) continue;
BRepAdaptor_Surface sf(face, Standard_True);
BRepLProp_SLProps prop(sf, 2, 1e-5);
int ntriangles = triangulation -> NbTriangles();
for (int j = 1; j <= ntriangles; j++)
{
gp_Pnt p[3];
gp_Pnt2d par[3];
for (int k = 1; k <=3; k++)
{
int n = triangulation->Triangles()(j)(k);
p[k-1] = triangulation->Nodes()(n).Transformed(loc);
par[k-1] = triangulation->UVNodes()(n);
}
//double maxside = 0;
//maxside = max (maxside, p[0].Distance(p[1]));
//maxside = max (maxside, p[0].Distance(p[2]));
//maxside = max (maxside, p[1].Distance(p[2]));
//cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush;
RestrictHTriangle (par[0], par[1], par[2], &prop, *mesh, 0);
//cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush;
}
}
// setting close edges
if (stlparam.resthcloseedgeenable)
{
multithread.task = "Setting local mesh size (close edges)";
int sections = 100;
Array<Line> lines(sections*nedges);
Box3dTree* searchtree =
new Box3dTree (bb.PMin(), bb.PMax());
int nlines = 0;
for (int i = 1; i <= nedges && !multithread.terminate; i++)
{
TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
if (BRep_Tool::Degenerated(edge)) continue;
double s0, s1;
Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);
BRepAdaptor_Curve brepc(edge);
BRepLProp_CLProps prop(brepc, 1, 1e-5);
prop.SetParameter (s0);
gp_Vec d0 = prop.D1().Normalized();
double s_start = s0;
int count = 0;
for (int j = 1; j <= sections; j++)
{
double s = s0 + (s1-s0)*(double)j/(double)sections;
prop.SetParameter (s);
gp_Vec d1 = prop.D1().Normalized();
double cosalpha = fabs(d0*d1);
if ((j == sections) || (cosalpha < cos(10.0/180.0*M_PI)))
{
count++;
gp_Pnt p0 = c->Value (s_start);
gp_Pnt p1 = c->Value (s);
lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z());
lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z());
Box3d box;
box.SetPoint (Point3d(lines[nlines].p0));
box.AddPoint (Point3d(lines[nlines].p1));
searchtree->Insert (box.PMin(), box.PMax(), nlines+1);
nlines++;
s_start = s;
d0 = d1;
}
}
}
Array<int> linenums;
for (int i = 0; i < nlines; i++)
{
multithread.percent = (100*i)/double(nlines);
Line & line = lines[i];
Box3d box;
box.SetPoint (Point3d(line.p0));
box.AddPoint (Point3d(line.p1));
double maxhline = max (mesh->GetH(box.PMin()),
mesh->GetH(box.PMax()));
box.Increase(maxhline);
double mindist = 1e99;
linenums.SetSize(0);
searchtree->GetIntersecting(box.PMin(),box.PMax(),linenums);
for (int j = 0; j < linenums.Size(); j++)
{
int num = linenums[j]-1;
if (i == num) continue;
if ((line.p0-lines[num].p0).Length2() < 1e-15) continue;
if ((line.p0-lines[num].p1).Length2() < 1e-15) continue;
if ((line.p1-lines[num].p0).Length2() < 1e-15) continue;
if ((line.p1-lines[num].p1).Length2() < 1e-15) continue;
mindist = min (mindist, line.Dist(lines[num]));
}
mindist *= stlparam.resthcloseedgefac;
if (mindist < 1e-3)
{
(*testout) << "extremely small local h: " << mindist
<< " --> setting to 1e-3" << endl;
(*testout) << "somewhere near " << line.p0 << " - " << line.p1 << endl;
mindist = 1e-3;
}
mesh->RestrictLocalHLine(line.p0, line.p1, mindist);
}
}
multithread.task = savetask;
}
}
if (multithread.terminate || perfstepsend <= MESHCONST_ANALYSE)
if (multithread.terminate || perfstepsend <= MESHCONST_ANALYSE)
return TCL_OK;
if (perfstepsstart <= MESHCONST_MESHEDGES)
if (perfstepsstart <= MESHCONST_MESHEDGES)
{
FindEdges (geom, *mesh);
FindEdges (geom, *mesh);
/*
cout << "Removing redundant points" << endl;
/*
cout << "Removing redundant points" << endl;
int i, j;
int np = mesh->GetNP();
Array<int> equalto;
int i, j;
int np = mesh->GetNP();
Array<int> equalto;
equalto.SetSize (np);
equalto = 0;
equalto.SetSize (np);
equalto = 0;
for (i = 1; i <= np; i++)
{
for (j = i+1; j <= np; j++)
{
if (!equalto[j-1] && (Dist2 (mesh->Point(i), mesh->Point(j)) < 1e-12))
equalto[j-1] = i;
}
}
for (i = 1; i <= np; i++)
{
for (j = i+1; j <= np; j++)
{
if (!equalto[j-1] && (Dist2 (mesh->Point(i), mesh->Point(j)) < 1e-12))
equalto[j-1] = i;
}
}
for (i = 1; i <= np; i++)
if (equalto[i-1])
{
cout << "Point " << i << " is equal to Point " << equalto[i-1] << endl;
for (j = 1; j <= mesh->GetNSeg(); j++)
{
Segment & seg = mesh->LineSegment(j);
if (seg.p1 == i) seg.p1 = equalto[i-1];
if (seg.p2 == i) seg.p2 = equalto[i-1];
}
}
for (i = 1; i <= np; i++)
if (equalto[i-1])
{
cout << "Point " << i << " is equal to Point " << equalto[i-1] << endl;
for (j = 1; j <= mesh->GetNSeg(); j++)
{
Segment & seg = mesh->LineSegment(j);
if (seg.p1 == i) seg.p1 = equalto[i-1];
if (seg.p2 == i) seg.p2 = equalto[i-1];
}
}
cout << "Removing degenerated segments" << endl;
for (j = 1; j <= mesh->GetNSeg(); j++)
{
Segment & seg = mesh->LineSegment(j);
if (seg.p1 == seg.p2)
{
mesh->DeleteSegment(j);
cout << "Deleting Segment " << j << endl;
}
}
mesh->Compress();
*/
/*
for (int i = 1; i <= geom.fmap.Extent(); i++)
{
Handle(Geom_Surface) hf1 =
BRep_Tool::Surface(TopoDS::Face(geom.fmap(i)));
for (int j = i+1; j <= geom.fmap.Extent(); j++)
{
Handle(Geom_Surface) hf2 =
BRep_Tool::Surface(TopoDS::Face(geom.fmap(j)));
if (hf1 == hf2) cout << "face " << i << " and face " << j << " lie on same surface" << endl;
}
}
*/
cout << "Removing degenerated segments" << endl;
for (j = 1; j <= mesh->GetNSeg(); j++)
{
Segment & seg = mesh->LineSegment(j);
if (seg.p1 == seg.p2)
{
mesh->DeleteSegment(j);
cout << "Deleting Segment " << j << endl;
}
}
mesh->Compress();
*/
/*
for (int i = 1; i <= geom.fmap.Extent(); i++)
{
Handle(Geom_Surface) hf1 =
BRep_Tool::Surface(TopoDS::Face(geom.fmap(i)));
for (int j = i+1; j <= geom.fmap.Extent(); j++)
{
Handle(Geom_Surface) hf2 =
BRep_Tool::Surface(TopoDS::Face(geom.fmap(j)));
if (hf1 == hf2) cout << "face " << i << " and face " << j << " lie on same surface" << endl;
}
}
*/
#ifdef LOG_STREAM
(*logout) << "Edges meshed" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
(*logout) << "Edges meshed" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
#endif
}
if (multithread.terminate || perfstepsend <= MESHCONST_MESHEDGES)
if (multithread.terminate || perfstepsend <= MESHCONST_MESHEDGES)
return TCL_OK;
if (perfstepsstart <= MESHCONST_MESHSURFACE)
if (perfstepsstart <= MESHCONST_MESHSURFACE)
{
OCCMeshSurface (geom, *mesh, perfstepsend);
if (multithread.terminate) return TCL_OK;
OCCMeshSurface (geom, *mesh, perfstepsend);
if (multithread.terminate) return TCL_OK;
#ifdef LOG_STREAM
(*logout) << "Surfaces meshed" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
(*logout) << "Surfaces meshed" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
#endif
#ifdef STAT_STREAM
(*statout) << mesh->GetNSeg() << " & "
<< mesh->GetNSE() << " & - &"
<< GetTime() << " & " << endl;
(*statout) << mesh->GetNSeg() << " & "
<< mesh->GetNSE() << " & - &"
<< GetTime() << " & " << endl;
#endif
// MeshQuality2d (*mesh);
mesh->CalcSurfacesOfNode();
// MeshQuality2d (*mesh);
mesh->CalcSurfacesOfNode();
}
if (multithread.terminate || perfstepsend <= MESHCONST_OPTSURFACE)
if (multithread.terminate || perfstepsend <= MESHCONST_OPTSURFACE)
return TCL_OK;
if (perfstepsstart <= MESHCONST_MESHVOLUME)
if (perfstepsstart <= MESHCONST_MESHVOLUME)
{
multithread.task = "Volume meshing";
multithread.task = "Volume meshing";
MESHING3_RESULT res =
MeshVolume (mparam, *mesh);
MESHING3_RESULT res =
MeshVolume (mparam, *mesh);
ofstream problemfile("occmesh.rep",ios_base::app);
ofstream problemfile("occmesh.rep",ios_base::app);
problemfile << "VOLUMEMESHING" << endl << endl;
if(res != MESHING3_OK)
problemfile << "ERROR" << endl << endl;
else
problemfile << "OK" << endl
<< mesh->GetNE() << " elements" << endl << endl;
problemfile << "VOLUMEMESHING" << endl << endl;
if(res != MESHING3_OK)
problemfile << "ERROR" << endl << endl;
else
problemfile << "OK" << endl
<< mesh->GetNE() << " elements" << endl << endl;
problemfile.close();
problemfile.close();
if (res != MESHING3_OK) return TCL_ERROR;
if (res != MESHING3_OK) return TCL_ERROR;
if (multithread.terminate) return TCL_OK;
if (multithread.terminate) return TCL_OK;
RemoveIllegalElements (*mesh);
if (multithread.terminate) return TCL_OK;
RemoveIllegalElements (*mesh);
if (multithread.terminate) return TCL_OK;
MeshQuality3d (*mesh);
MeshQuality3d (*mesh);
#ifdef STAT_STREAM
(*statout) << GetTime() << " & ";
(*statout) << GetTime() << " & ";
#endif
#ifdef LOG_STREAM
(*logout) << "Volume meshed" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
(*logout) << "Volume meshed" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
#endif
}
if (multithread.terminate || perfstepsend <= MESHCONST_MESHVOLUME)
if (multithread.terminate || perfstepsend <= MESHCONST_MESHVOLUME)
return TCL_OK;
if (perfstepsstart <= MESHCONST_OPTVOLUME)
if (perfstepsstart <= MESHCONST_OPTVOLUME)
{
multithread.task = "Volume optimization";
multithread.task = "Volume optimization";
OptimizeVolume (mparam, *mesh);
if (multithread.terminate) return TCL_OK;
OptimizeVolume (mparam, *mesh);
if (multithread.terminate) return TCL_OK;
#ifdef STAT_STREAM
(*statout) << GetTime() << " & "
<< mesh->GetNE() << " & "
<< mesh->GetNP() << " " << '\\' << '\\' << " \\" << "hline" << endl;
(*statout) << GetTime() << " & "
<< mesh->GetNE() << " & "
<< mesh->GetNP() << " " << '\\' << '\\' << " \\" << "hline" << endl;
#endif
#ifdef LOG_STREAM
(*logout) << "Volume optimized" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
(*logout) << "Volume optimized" << endl
<< "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl;
#endif
// cout << "Optimization complete" << endl;
// cout << "Optimization complete" << endl;
}
(*testout) << "NP: " << mesh->GetNP() << endl;
for (int i = 1; i <= mesh->GetNP(); i++)
(*testout) << "NP: " << mesh->GetNP() << endl;
for (int i = 1; i <= mesh->GetNP(); i++)
(*testout) << mesh->Point(i) << endl;
(*testout) << endl << "NSegments: " << mesh->GetNSeg() << endl;
for (int i = 1; i <= mesh->GetNSeg(); i++)
(*testout) << endl << "NSegments: " << mesh->GetNSeg() << endl;
for (int i = 1; i <= mesh->GetNSeg(); i++)
(*testout) << mesh->LineSegment(i) << endl;
return TCL_OK;
}
return TCL_OK;
}
}
#endif

View File

@ -841,6 +841,14 @@ namespace netgen
facemeshstatus.SetSize (fmap.Extent());
facemeshstatus = 0;
// Philippose - 15/01/2009
face_maxh.SetSize (fmap.Extent());
face_maxh = mparam.maxh;
// Philippose - 17/01/2009
face_sel_status.SetSize (fmap.Extent());
face_sel_status = 0;
fvispar.SetSize (fmap.Extent());
evispar.SetSize (emap.Extent());
vvispar.SetSize (vmap.Extent());
@ -1127,59 +1135,135 @@ namespace netgen
return occgeo;
}
OCCGeometry * LoadOCC_STEP (const char * filename)
{
OCCGeometry * occgeo;
occgeo = new OCCGeometry;
STEPControl_Reader reader;
Standard_Integer stat = reader.ReadFile((char*)filename);
Standard_Integer nb = reader.NbRootsForTransfer();
reader.TransferRoots (); // Tranlate STEP -> OCC
// Philippose - 29/01/2009
/* Special STEP File load function including the ability
to extract individual surface colours via the extended
OpenCascade XDE and XCAF Feature set.
*/
OCCGeometry * LoadOCC_STEP (const char * filename)
{
OCCGeometry * occgeo;
occgeo = new OCCGeometry;
// Initiate a dummy XCAF Application to handle the STEP XCAF Document
static Handle_XCAFApp_Application dummy_app = XCAFApp_Application::GetApplication();
// Create an XCAF Document to contain the STEP file itself
Handle_TDocStd_Document step_doc;
// Check if a STEP File is already open under this handle, if so, close it to prevent
// Segmentation Faults when trying to create a new document
if(dummy_app->NbDocuments() > 0)
{
dummy_app->GetDocument(1,step_doc);
dummy_app->Close(step_doc);
}
dummy_app->NewDocument ("STEP-XCAF",step_doc);
STEPCAFControl_Reader reader;
Standard_Integer stat = reader.ReadFile((char*)filename);
// Enable transfer of colours
reader.SetColorMode(Standard_True);
reader.Transfer(step_doc);
// Read in the shape(s) and the colours present in the STEP File
Handle_XCAFDoc_ShapeTool step_shape_contents = XCAFDoc_DocumentTool::ShapeTool(step_doc->Main());
Handle_XCAFDoc_ColorTool step_colour_contents = XCAFDoc_DocumentTool::ColorTool(step_doc->Main());
TDF_LabelSequence step_shapes;
step_shape_contents->GetShapes(step_shapes);
/*
// List out the available colours in the STEP File as Colour Names
TDF_LabelSequence allColours;
stepColourContents->GetColors(allColours);
cout << "Number of colours in STEP = " << allColours.Length() << endl;
for(int i = 1; i <= allColours.Length(); i++)
{
Quantity_Color col;
stepColourContents->GetColor(allColours.Value(i),col);
cout << "Colour [" << i << "] = " << col.StringName(col.Name()) << endl;
}
*/
occgeo->shape = step_shape_contents->GetShape(step_shapes.Value(1));
occgeo->face_colours = step_colour_contents;
occgeo->changed = 1;
occgeo->BuildFMap();
occgeo->BuildVisualizationMesh();
PrintContents (occgeo);
return occgeo;
}
occgeo->shape = reader.OneShape();
occgeo->changed = 1;
occgeo->BuildFMap();
//
//Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
//sfs->Init(occgeo->shape);
//sfs->Perform();
//Handle(ShapeFix_Wireframe) sfwf = new ShapeFix_Wireframe(occgeo->shape);
//sfwf->FixSmallEdges();
//sfwf->FixWireGaps();
// Philippose - 29/01/2009
// The LOADOCC_STEP Function has been replaced by the one
// above, which also includes support for the OpenCascade
// XDE Features.
//
// OCCGeometry * LoadOCC_STEP (const char * filename)
// {
// OCCGeometry * occgeo;
// occgeo = new OCCGeometry;
//
// STEPControl_Reader reader;
// Standard_Integer stat = reader.ReadFile((char*)filename);
// Standard_Integer nb = reader.NbRootsForTransfer();
// reader.TransferRoots (); // Tranlate STEP -> OCC
//
//
//
// occgeo->shape = reader.OneShape();
// occgeo->changed = 1;
// occgeo->BuildFMap();
// //
// //Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
// //sfs->Init(occgeo->shape);
// //sfs->Perform();
// //Handle(ShapeFix_Wireframe) sfwf = new ShapeFix_Wireframe(occgeo->shape);
// //sfwf->FixSmallEdges();
// //sfwf->FixWireGaps();
//
//
//
// /*
// // JS
// TopoDS_Compound aRes;
// BRep_Builder aBuilder;
// aBuilder.MakeCompound (aRes);
//
// for (TopExp_Explorer exp(occgeo->shape, TopAbs_SOLID); exp.More(); exp.Next())
// {
// aBuilder.Add (aRes, exp.Current());
// cout << "solid" << endl;
// }
//
// for (TopExp_Explorer exp(aRes, TopAbs_SOLID); exp.More(); exp.Next())
// {
// cout << "compound has shapes solid" << endl;
// }
// occgeo->shape = aRes;
// occgeo->changed = 1;
// occgeo->BuildFMap();
// */
//
//
// occgeo->BuildVisualizationMesh();
// PrintContents (occgeo);
//
// return occgeo;
// }
/*
// JS
TopoDS_Compound aRes;
BRep_Builder aBuilder;
aBuilder.MakeCompound (aRes);
for (TopExp_Explorer exp(occgeo->shape, TopAbs_SOLID); exp.More(); exp.Next())
{
aBuilder.Add (aRes, exp.Current());
cout << "solid" << endl;
}
for (TopExp_Explorer exp(aRes, TopAbs_SOLID); exp.More(); exp.Next())
{
cout << "compound has shapes solid" << endl;
}
occgeo->shape = aRes;
occgeo->changed = 1;
occgeo->BuildFMap();
*/
//
occgeo->BuildVisualizationMesh();
PrintContents (occgeo);
return occgeo;
}
OCCGeometry * LoadOCC_BREP (const char * filename)
{

View File

@ -11,9 +11,9 @@
#include <meshing.hpp>
#include <BRep_Tool.hxx>
#include <Geom_Curve.hxx>
#include <Geom2d_Curve.hxx>
#include "BRep_Tool.hxx"
#include "Geom_Curve.hxx"
#include "Geom2d_Curve.hxx"
#include "Geom_Surface.hxx"
#include "GeomAPI_ProjectPointOnSurf.hxx"
#include "GeomAPI_ProjectPointOnCurve.hxx"
@ -37,6 +37,7 @@
#include "TopoDS.hxx"
#include "TopoDS_Solid.hxx"
#include "TopExp_Explorer.hxx"
#include "TopTools_ListIteratorOfListOfShape.hxx"
#include "BRep_Tool.hxx"
#include "Geom_Curve.hxx"
#include "Geom2d_Curve.hxx"
@ -62,8 +63,6 @@
#include "Poly_Triangle.hxx"
#include "GProp_GProps.hxx"
#include "BRepGProp.hxx"
#include "IGESControl_Reader.hxx"
#include "STEPControl_Reader.hxx"
#include "TopoDS_Shape.hxx"
#include "TopoDS_Face.hxx"
#include "IGESToBRep_Reader.hxx"
@ -80,8 +79,33 @@
#include "Bnd_Box.hxx"
#include "ShapeAnalysis.hxx"
#include "ShapeBuild_ReShape.hxx"
// Philippose - 29/01/2009
// OpenCascade XDE Support
// Include support for OpenCascade XDE Features
#include "TDocStd_Document.hxx"
#include "Quantity_Color.hxx"
#include "XCAFApp_Application.hxx"
#include "XCAFDoc_ShapeTool.hxx"
#include "XCAFDoc_Color.hxx"
#include "XCAFDoc_ColorTool.hxx"
#include "XCAFDoc_ColorType.hxx"
#include "XCAFDoc_LayerTool.hxx"
#include "XCAFDoc_DimTolTool.hxx"
#include "XCAFDoc_MaterialTool.hxx"
#include "XCAFDoc_DocumentTool.hxx"
#include "TDF_Label.hxx"
#include "TDF_LabelSequence.hxx"
#include "STEPCAFControl_Reader.hxx"
#include "STEPCAFControl_Writer.hxx"
#include "IGESCAFControl_Reader.hxx"
#include "IGESCAFControl_Writer.hxx"
#include "IGESControl_Reader.hxx"
#include "STEPControl_Reader.hxx"
#include "IGESControl_Writer.hxx"
#include "STEPControl_Writer.hxx"
#include "StlAPI_Writer.hxx"
#include "STEPControl_StepModelType.hxx"
@ -89,175 +113,241 @@ namespace netgen
{
#include "../visualization/vispar.hpp"
// class VisualizationParameters;
// extern VisualizationParameters vispar;
// class VisualizationParameters;
// extern VisualizationParameters vispar;
#include "occmeshsurf.hpp"
#define PROJECTION_TOLERANCE 1e-10
#define ENTITYISVISIBLE 1
#define ENTITYISHIGHLIGHTED 2
#define ENTITYISDRAWABLE 4
class EntityVisualizationCode
{
int code;
class EntityVisualizationCode
{
int code;
public:
public:
EntityVisualizationCode()
{ code = ENTITYISVISIBLE + !ENTITYISHIGHLIGHTED + ENTITYISDRAWABLE; }
EntityVisualizationCode()
{ code = ENTITYISVISIBLE + !ENTITYISHIGHLIGHTED + ENTITYISDRAWABLE;}
int IsVisible ()
{ return code & ENTITYISVISIBLE; }
int IsVisible ()
{ return code & ENTITYISVISIBLE;}
int IsHighlighted ()
{ return code & ENTITYISHIGHLIGHTED; }
int IsHighlighted ()
{ return code & ENTITYISHIGHLIGHTED;}
int IsDrawable ()
{ return code & ENTITYISDRAWABLE; }
int IsDrawable ()
{ return code & ENTITYISDRAWABLE;}
void Show ()
{ code |= ENTITYISVISIBLE; }
void Show ()
{ code |= ENTITYISVISIBLE;}
void Hide ()
{ code &= ~ENTITYISVISIBLE; }
void Hide ()
{ code &= ~ENTITYISVISIBLE;}
void Highlight ()
{ code |= ENTITYISHIGHLIGHTED; }
void Highlight ()
{ code |= ENTITYISHIGHLIGHTED;}
void Lowlight ()
{ code &= ~ENTITYISHIGHLIGHTED; }
void Lowlight ()
{ code &= ~ENTITYISHIGHLIGHTED;}
void SetDrawable ()
{ code |= ENTITYISDRAWABLE; }
void SetNotDrawable ()
{ code &= ~ENTITYISDRAWABLE; }
};
inline double Det3 (double a00, double a01, double a02,
double a10, double a11, double a12,
double a20, double a21, double a22)
{
return a00*a11*a22 + a01*a12*a20 + a10*a21*a02 - a20*a11*a02 - a10*a01*a22 - a21*a12*a00;
}
void SetDrawable ()
{ code |= ENTITYISDRAWABLE;}
void SetNotDrawable ()
{ code &= ~ENTITYISDRAWABLE;}
};
inline double Det3 (double a00, double a01, double a02,
double a10, double a11, double a12,
double a20, double a21, double a22)
{
return a00*a11*a22 + a01*a12*a20 + a10*a21*a02 - a20*a11*a02 - a10*a01*a22 - a21*a12*a00;
}
#define OCCGEOMETRYVISUALIZATIONNOCHANGE 0
#define OCCGEOMETRYVISUALIZATIONFULLCHANGE 1
// == compute transformation matrices and redraw
// == compute transformation matrices and redraw
#define OCCGEOMETRYVISUALIZATIONHALFCHANGE 2
// == redraw
// == redraw
class OCCGeometry
{
Point<3> center;
class OCCGeometry
{
Point<3> center;
public:
TopoDS_Shape shape;
TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap;
Array<bool> fsingular, esingular, vsingular;
Box<3> boundingbox;
public:
TopoDS_Shape shape;
TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap;
Array<bool> fsingular, esingular, vsingular;
Box<3> boundingbox;
int changed;
Array<int> facemeshstatus;
// Philippose - 29/01/2009
// OpenCascade XDE Support
// XCAF Handle to make the face colours available to the rest of
// the system
Handle_XCAFDoc_ColorTool face_colours;
Array<EntityVisualizationCode> fvispar, evispar, vvispar;
int changed;
Array<int> facemeshstatus;
double tolerance;
bool fixsmalledges;
bool fixspotstripfaces;
bool sewfaces;
bool makesolids;
bool splitpartitions;
// Philippose - 15/01/2009
// Maximum mesh size for a given face
// (Used to explicitly define mesh size limits on individual faces)
Array<double> face_maxh;
// Philippose - 15/01/2009
// Indicates which faces have been selected by the user in geometry mode
// (Currently handles only selection of one face at a time, but an array would
// help to extend this to multiple faces)
Array<bool> face_sel_status;
OCCGeometry()
{
somap.Clear();
shmap.Clear();
fmap.Clear();
wmap.Clear();
emap.Clear();
vmap.Clear();
}
Array<EntityVisualizationCode> fvispar, evispar, vvispar;
double tolerance;
bool fixsmalledges;
bool fixspotstripfaces;
bool sewfaces;
bool makesolids;
bool splitpartitions;
void BuildFMap();
OCCGeometry()
{
somap.Clear();
shmap.Clear();
fmap.Clear();
wmap.Clear();
emap.Clear();
vmap.Clear();
}
Box<3> GetBoundingBox()
{ return boundingbox; }
void BuildFMap();
int NrSolids()
{ return somap.Extent(); }
Box<3> GetBoundingBox()
{ return boundingbox;}
void SetCenter()
{ center = boundingbox.Center(); }
int NrSolids()
{ return somap.Extent();}
Point<3> Center()
{ return center; }
// Philippose - 17/01/2009
// Total number of faces in the geometry
int NrFaces()
{ return fmap.Extent();}
void Project (int surfi, Point<3> & p) const;
bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
void SetCenter()
{ center = boundingbox.Center();}
Point<3> Center()
{ return center;}
OCCSurface GetSurface (int surfi)
{
cout << "OCCGeometry::GetSurface using PLANESPACE" << endl;
return OCCSurface (TopoDS::Face(fmap(surfi)), PLANESPACE);
}
void Project (int surfi, Point<3> & p) const;
bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
OCCSurface GetSurface (int surfi)
{
cout << "OCCGeometry::GetSurface using PLANESPACE" << endl;
return OCCSurface (TopoDS::Face(fmap(surfi)), PLANESPACE);
}
void BuildVisualizationMesh ();
void BuildVisualizationMesh ();
void RecursiveTopologyTree (const TopoDS_Shape & sh,
stringstream & str,
TopAbs_ShapeEnum l,
bool free,
const char * lname);
void RecursiveTopologyTree (const TopoDS_Shape & sh,
stringstream & str,
TopAbs_ShapeEnum l,
bool free,
const char * lname);
void GetTopologyTree (stringstream & str);
void GetTopologyTree (stringstream & str);
void PrintNrShapes ();
void PrintNrShapes ();
void CheckIrregularEntities (stringstream & str);
void CheckIrregularEntities (stringstream & str);
void SewFaces();
void SewFaces();
void MakeSolid();
void MakeSolid();
void HealGeometry();
void HealGeometry();
void LowLightAll()
{
for (int i = 1; i <= fmap.Extent(); i++)
fvispar[i-1].Lowlight();
for (int i = 1; i <= emap.Extent(); i++)
evispar[i-1].Lowlight();
for (int i = 1; i <= vmap.Extent(); i++)
vvispar[i-1].Lowlight();
}
// Philippose - 15/01/2009
// Sets the maximum mesh size for a given face
// (Note: Local mesh size limited by the global max mesh size)
void SetFaceMaxH(int facenr, double faceh)
{
if((facenr> 0) && (facenr <= fmap.Extent()))
{
face_maxh[facenr-1] = min(mparam.maxh,faceh);
}
}
void GetUnmeshedFaceInfo (stringstream & str);
void GetNotDrawableFaces (stringstream & str);
bool ErrorInSurfaceMeshing ();
// Philippose - 15/01/2009
// Returns the local mesh size of a given face
double GetFaceMaxH(int facenr)
{
if((facenr> 0) && (facenr <= fmap.Extent()))
{
return face_maxh[facenr-1];
}
else
{
return 0.0;
}
}
void WriteOCC_STL(char * filename);
};
// Philippose - 17/01/2009
// Returns the index of the currently selected face
int SelectedFace()
{
int i;
for(i = 1; i <= fmap.Extent(); i++)
{
if(face_sel_status[i-1])
{
return i;
}
}
void PrintContents (OCCGeometry * geom);
return 0;
}
OCCGeometry * LoadOCC_IGES (const char * filename);
OCCGeometry * LoadOCC_STEP (const char * filename);
OCCGeometry * LoadOCC_BREP (const char * filename);
// Philippose - 17/01/2009
// Sets the currently selected face
void SetSelectedFace(int facenr)
{
face_sel_status = 0;
if((facenr >= 1) && (facenr <= fmap.Extent()))
{
face_sel_status[facenr-1] = 1;
}
}
void LowLightAll()
{
for (int i = 1; i <= fmap.Extent(); i++)
fvispar[i-1].Lowlight();
for (int i = 1; i <= emap.Extent(); i++)
evispar[i-1].Lowlight();
for (int i = 1; i <= vmap.Extent(); i++)
vvispar[i-1].Lowlight();
}
void GetUnmeshedFaceInfo (stringstream & str);
void GetNotDrawableFaces (stringstream & str);
bool ErrorInSurfaceMeshing ();
void WriteOCC_STL(char * filename);
};
void PrintContents (OCCGeometry * geom);
OCCGeometry * LoadOCC_IGES (const char * filename);
OCCGeometry * LoadOCC_STEP (const char * filename);
OCCGeometry * LoadOCC_BREP (const char * filename);
}

View File

@ -9,18 +9,25 @@
#include <csg.hpp>
#include <stlgeom.hpp>
#include <visual.hpp>
// #include <parallel.hpp>
#ifdef OCCGEOMETRY
// Philippose - 30/01/2009
// Required for OpenCascade XDE Support
#include <occgeom.hpp>
#endif
#include <visual.hpp>
namespace netgen
{
#ifdef OCCGEOMETRY
// Philippose - 30/01/2009
// Required for OpenCascade XDE Support
extern OCCGeometry * occgeometry;
#endif
extern AutoPtr<Mesh> mesh;
@ -1025,6 +1032,23 @@ namespace netgen
}
#endif
#ifdef OCCGEOMETRY
// Philippose - 30/01/2009
// OpenCascade XDE Support
// Update the colour of each face based on the STEP File Data
// if the advanced OpenCascade XDE Support has been enabled
if((col == 1) && (occgeometry))
{
TopoDS_Face face = TopoDS::Face(occgeometry->fmap(el.GetIndex()));
Quantity_Color face_colour;
occgeometry->face_colours->GetColor(face,XCAFDoc_ColorSurf,face_colour);
matcol[0] = face_colour.Red();
matcol[1] = face_colour.Green();
matcol[2] = face_colour.Blue();
matcol[3] = 1.0;
glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcol);
}
#endif
bool drawel = !el.IsDeleted();

File diff suppressed because it is too large Load Diff