From 4be31562cb6ecfacfea99a09553ab0c6bd79fdf4 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 12 May 2010 09:42:51 +0000
Subject: [PATCH] the solid modeler is now beginning to be usable : OCC :
Pipes, NURBS, Constraint surfaces, planar surfaces, ruled surfaces are there
The gluer is less bugged (yet not really usable) ACIS : interface is a
little more complete brep IO's FIXED for OSX
---
Common/Gmsh.cpp | 2 +-
Geo/GFace.cpp | 82 ++++++++++++++++++
Geo/GModel.cpp | 48 ++++++++++-
Geo/GModel.h | 4 +
Geo/GModelFactory.cpp | 149 ++++++++++++++++++++++++++++++++-
Geo/GModelFactory.h | 15 ++++
Geo/GModelIO_OCC.cpp | 22 ++++-
Geo/GRegion.cpp | 5 ++
Geo/OCCRegion.cpp | 3 +-
Solver/solverAlgorithms.h | 26 ++++++
benchmarks/boolean/square1.lua | 34 +++++++-
11 files changed, 375 insertions(+), 15 deletions(-)
diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index ed28f1559a..39be0d9578 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -167,7 +167,7 @@ int GmshBatch()
GModel::current()->checkMeshCoherence(CTX::instance()->geom.tolerance);
}
else if(CTX::instance()->batch == -1){
- CreateOutputFile(CTX::instance()->outputFileName, FORMAT_GEO);
+ CreateOutputFile(CTX::instance()->outputFileName, FORMAT_AUTO);
}
else if(CTX::instance()->batch > 0){
#if defined(HAVE_MESH)
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 6e2174bcce..f1572a01fa 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -484,6 +484,88 @@ end:
}
}
+void computeMeanPlane(const std::vector<SPoint3> &points, mean_plane &meanPlane)
+{
+ double xm = 0., ym = 0., zm = 0.;
+ int ndata = points.size();
+ int na = 3;
+ for(int i = 0; i < ndata; i++) {
+ xm += points[i].x();
+ ym += points[i].y();
+ zm += points[i].z();
+ }
+ xm /= (double)ndata;
+ ym /= (double)ndata;
+ zm /= (double)ndata;
+
+ fullMatrix<double> U(ndata, na), V(na, na);
+ fullVector<double> sigma(na);
+ for(int i = 0; i < ndata; i++) {
+ U(i, 0) = points[i].x() - xm;
+ U(i, 1) = points[i].y() - ym;
+ U(i, 2) = points[i].z() - zm;
+ }
+ U.svd(V, sigma);
+ double res[4], svd[3];
+ svd[0] = sigma(0);
+ svd[1] = sigma(1);
+ svd[2] = sigma(2);
+ int min;
+ if(fabs(svd[0]) < fabs(svd[1]) && fabs(svd[0]) < fabs(svd[2]))
+ min = 0;
+ else if(fabs(svd[1]) < fabs(svd[0]) && fabs(svd[1]) < fabs(svd[2]))
+ min = 1;
+ else
+ min = 2;
+ res[0] = V(0, min);
+ res[1] = V(1, min);
+ res[2] = V(2, min);
+ norme(res);
+
+ double ex[3], t1[3], t2[3];
+
+ ex[0] = ex[1] = ex[2] = 0.0;
+ if(res[0] == 0.)
+ ex[0] = 1.0;
+ else if(res[1] == 0.)
+ ex[1] = 1.0;
+ else
+ ex[2] = 1.0;
+
+ prodve(res, ex, t1);
+ norme(t1);
+ prodve(t1, res, t2);
+ norme(t2);
+
+ res[3] = (xm * res[0] + ym * res[1] + zm * res[2]);
+
+ for(int i = 0; i < 3; i++)
+ meanPlane.plan[0][i] = t1[i];
+ for(int i = 0; i < 3; i++)
+ meanPlane.plan[1][i] = t2[i];
+ for(int i = 0; i < 3; i++)
+ meanPlane.plan[2][i] = res[i];
+
+ meanPlane.a = res[0];
+ meanPlane.b = res[1];
+ meanPlane.c = res[2];
+ meanPlane.d = res[3];
+
+ meanPlane.x = meanPlane.y = meanPlane.z = 0.;
+ if(fabs(meanPlane.a) >= fabs(meanPlane.b) &&
+ fabs(meanPlane.a) >= fabs(meanPlane.c) ){
+ meanPlane.x = meanPlane.d / meanPlane.a;
+ }
+ else if(fabs(meanPlane.b) >= fabs(meanPlane.a) &&
+ fabs(meanPlane.b) >= fabs(meanPlane.c)){
+ meanPlane.y = meanPlane.d / meanPlane.b;
+ }
+ else{
+ meanPlane.z = meanPlane.d / meanPlane.c;
+ }
+}
+
+
void GFace::getMeanPlaneData(double VX[3], double VY[3],
double &x, double &y, double &z) const
{
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index fc92ff85db..2b873a0a60 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1470,6 +1470,23 @@ GEdge *GModel::addNURBS(GVertex *start, GVertex *end,
return 0;
}
+void GModel::addRuledFaces (std::vector<std::vector<GEdge *> > edges){
+ if(_factory)
+ _factory->addRuledFaces(this, edges);
+}
+
+GFace* GModel::addFace (std::vector<GEdge *> edges, std::vector< std::vector<double > > points){
+ if(_factory)
+ return _factory->addFace(this, edges, points);
+ return 0;
+}
+
+GFace* GModel::addPlanarFace (std::vector<std::vector<GEdge *> > edges){
+ if(_factory)
+ return _factory->addPlanarFace(this, edges);
+ return 0;
+}
+
GEntity *GModel::revolve(GEntity *e, std::vector<double> p1, std::vector<double> p2,
double angle)
{
@@ -1484,6 +1501,12 @@ GEntity *GModel::extrude(GEntity *e, std::vector<double> p1, std::vector<double>
return _factory->extrude(this, e, p1, p2);
return 0;
}
+GEntity *GModel::addPipe(GEntity *e, std::vector<GEdge *> edges){
+ if(_factory)
+ return _factory->addPipe(this,e,edges);
+ return 0;
+}
+
GEntity *GModel::addSphere(double cx, double cy, double cz, double radius)
{
@@ -1735,7 +1758,12 @@ static void glueFacesInRegions(GModel *model,
bool aDifferenceExists = false;
std::list<GFace*> old = r->faces(), fnew;
for (std::list<GFace*>::iterator fit = old.begin(); fit != old.end(); fit++){
- GFace *temp = Duplicates2Unique[*fit];
+ std::map<GFace*, GFace*>::iterator itR = Duplicates2Unique.find(*fit);
+ if (itR == Duplicates2Unique.end())
+ {
+ Msg::Fatal("error in the gluing process");
+ }
+ GFace *temp = itR->second;;
fnew.push_back(temp);
if (temp != *fit){
aDifferenceExists = true;
@@ -1834,6 +1862,18 @@ void GModel::registerBindings(binding *b)
cm = cb->addMethod("addNURBS", &GModel::addNURBS);
cm->setDescription("creates a NURBS curve from v1 to v2 with control Points, knots, weights and multiplicities");
cm->setArgNames("v1", "v2", "{{poles}}","{knots}","{weights}","{mult}",NULL);
+ cm = cb->addMethod("addFace", &GModel::addFace);
+ cm->setDescription("creates a face that is constraint by edges and points");
+ cm->setArgNames("{list of edges}","{{x,y,z},{x,y,z},...}",NULL);
+ cm = cb->addMethod("addRuledFaces", &GModel::addRuledFaces);
+ cm->setDescription("creates ruled faces that contains a list of wires");
+ cm->setArgNames("{{list of edges},{list of edges},...}",NULL);
+ cm = cb->addMethod("addPlanarFace", &GModel::addPlanarFace);
+ cm->setDescription("creates a planar face that contains a list of wires");
+ cm->setArgNames("{{list of edges},{list of edges},...}",NULL);
+ cm = cb->addMethod("addPipe", &GModel::addPipe);
+ cm->setDescription("creates a pipe with a base and a wire for the direction");
+ cm->setArgNames("ge","{list of edges}",NULL);
cm = cb->addMethod("addCircleArcCenter", &GModel::addCircleArcCenter);
cm->setDescription("create a circle arc going from v1 to v2 with its center "
"at (x,y,z)");
@@ -1878,9 +1918,9 @@ void GModel::registerBindings(binding *b)
"one (tool). The third parameter tells if a new model has to "
"be created");
cm->setArgNames("tool","createNewGModel",NULL);
- // cm = cb->addMethod("glue", &GModel::glue);
- // cm->setDescription("glue the geometric model using geometric tolerance eps");
- // cm->setArgNames("eps",NULL);
+ cm = cb->addMethod("glue", &GModel::glue);
+ cm->setDescription("glue the geometric model using geometric tolerance eps");
+ cm->setArgNames("eps",NULL);
cm = cb->addMethod("setAsCurrent", &GModel::setAsCurrent);
cm->setDescription("set the model as the current (active) one");
cm = cb->setConstructor<GModel>();
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 4c2fb7463d..8d48df7e76 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -376,6 +376,10 @@ class GModel
GEntity *revolve(GEntity *e, std::vector<double> p1, std::vector<double> p2,
double angle);
GEntity *extrude(GEntity *e, std::vector<double> p1, std::vector<double> p2);
+ GEntity *addPipe(GEntity *e, std::vector<GEdge *> edges);
+ void addRuledFaces (std::vector<std::vector<GEdge *> > edges);
+ GFace* addFace (std::vector<GEdge *> edges, std::vector< std::vector<double > > points);
+ GFace* addPlanarFace (std::vector<std::vector<GEdge *> > edges);
// create solid geometry primitives using the factory
GEntity *addSphere(double cx, double cy, double cz, double radius);
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 559590e8af..ac016c4fb7 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -11,11 +11,17 @@
#include "GModelIO_OCC.h"
#include "OCCVertex.h"
#include "OCCEdge.h"
+#include "OCCFace.h"
#include "OCCRegion.h"
#include "BRepBuilderAPI_MakeVertex.hxx"
+#include "BRepOffsetAPI_MakePipe.hxx"
#include "BRepBuilderAPI_MakeEdge.hxx"
#include "BRepPrimAPI_MakeRevol.hxx"
#include "BRepPrimAPI_MakePrism.hxx"
+#include "BRepBuilderAPI_MakeWire.hxx"
+#include "BRepOffsetAPI_ThruSections.hxx"
+#include "BRepOffsetAPI_MakeFilling.hxx"
+#include "BRepBuilderAPI_MakeFace.hxx"
#include <GC_MakeArcOfCircle.hxx>
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>
@@ -170,8 +176,8 @@ GEdge *OCCFactory::addNURBS(GModel *gm, GVertex *start, GVertex *end,
}
const int degree = totKnots - nbControlPoints - 1;
- printf("creation of a nurbs of degree %d with %d control points\n",
- degree,nbControlPoints);
+ Msg::Debug("creation of a nurbs of degree %d with %d control points",
+ degree,nbControlPoints);
int index = 1;
ctrlPoints.SetValue(index++, gp_Pnt(start->x(), start->y(), start->z()));
@@ -544,4 +550,143 @@ void OCCFactory::rotate(GModel *gm, std::vector<double> p1,std::vector<double> p
gm->_occ_internals->buildGModel(gm);
}
+
+std::vector<GFace *> OCCFactory::addRuledFaces (GModel *gm, std::vector< std::vector<GEdge *> > wires){
+
+ std::vector<GFace*> faces;
+ Standard_Boolean anIsSolid = Standard_False;
+ Standard_Boolean anIsRuled = Standard_True;
+ BRepOffsetAPI_ThruSections aGenerator (anIsSolid,anIsRuled);
+
+ for (int i=0;i<wires.size();i++) {
+ BRepBuilderAPI_MakeWire wire_maker;
+ for (int j=0;j<wires[i].size();j++) {
+ GEdge *ge = wires[i][j];
+ printf("edge %d\n",ge->tag());
+ OCCEdge *occe = dynamic_cast<OCCEdge*>(ge);
+ if (occe){
+ wire_maker.Add(occe->getTopoDS_Edge());
+ }
+ }
+ aGenerator.AddWire (wire_maker.Wire());
+ }
+
+ aGenerator.CheckCompatibility (Standard_False);
+
+ aGenerator.Build();
+
+ TopoDS_Shape aResult = aGenerator.Shape();
+
+ TopExp_Explorer exp2;
+ for(exp2.Init(TopoDS::Shell(aResult), TopAbs_FACE); exp2.More(); exp2.Next()){
+ TopoDS_Face face = TopoDS::Face(exp2.Current());
+ GFace *ret = gm->_occ_internals->addFaceToModel(gm, face);
+ faces.push_back(ret);
+ }
+ return faces;
+}
+
+GFace * OCCFactory::addFace (GModel *gm, std::vector<GEdge *> edges, std::vector< std::vector<double > > points){
+
+ BRepOffsetAPI_MakeFilling aGenerator;
+
+ for (int i=0;i<edges.size();i++) {
+ GEdge *ge = edges[i];
+ OCCEdge *occe = dynamic_cast<OCCEdge*>(ge);
+ if (occe){
+ aGenerator.Add(occe->getTopoDS_Edge(),GeomAbs_C0);
+ }
+ }
+ for (int i=0;i<points.size();i++) {
+ gp_Pnt aPnt (points[i][0],points[i][1],points[i][2]);
+ aGenerator.Add(aPnt);
+ }
+
+ aGenerator.Build();
+
+ TopoDS_Shape aResult = aGenerator.Shape();
+ return gm->_occ_internals->addFaceToModel(gm, TopoDS::Face(aResult));
+}
+
+extern void computeMeanPlane(const std::vector<SPoint3> &points, mean_plane &meanPlane);
+
+GFace * OCCFactory::addPlanarFace (GModel *gm, std::vector< std::vector<GEdge *> > wires){
+
+ std::set<GVertex*> verts;
+ for (int i=0;i<wires.size();i++) {
+ for (int j=0;j<wires[i].size();j++) {
+ GEdge *ge = wires[i][j];
+ verts.insert(ge->getBeginVertex());
+ verts.insert(ge->getEndVertex());
+ }
+ }
+ std::vector<SPoint3> points;
+ std::set<GVertex*>::iterator it = verts.begin();
+ for ( ; it != verts.end(); ++it){
+ points.push_back(SPoint3((*it)->x(),(*it)->y(),(*it)->z()));
+ }
+ mean_plane meanPlane;
+ computeMeanPlane(points, meanPlane);
+
+ gp_Pln aPlane (meanPlane.a,meanPlane.b,meanPlane.c,meanPlane.d);
+ BRepBuilderAPI_MakeFace aGenerator (aPlane);
+
+ for (int i=0;i<wires.size();i++) {
+ BRepBuilderAPI_MakeWire wire_maker;
+ for (int j=0;j<wires[i].size();j++) {
+ GEdge *ge = wires[i][j];
+ OCCEdge *occe = dynamic_cast<OCCEdge*>(ge);
+ if (occe){
+ wire_maker.Add(occe->getTopoDS_Edge());
+ }
+ }
+ TopoDS_Wire myWire = wire_maker.Wire();
+ if (i)myWire.Reverse();
+ aGenerator.Add (myWire);
+ }
+
+ aGenerator.Build();
+
+ TopoDS_Shape aResult = aGenerator.Shape();
+ return gm->_occ_internals->addFaceToModel(gm, TopoDS::Face(aResult));
+}
+
+GEntity * OCCFactory::addPipe (GModel *gm, GEntity *base, std::vector<GEdge *> wire){
+
+ BRepBuilderAPI_MakeWire wire_maker;
+ for (int j=0;j<wire.size();j++) {
+ GEdge *ge = wire[j];
+ OCCEdge *occe = dynamic_cast<OCCEdge*>(ge);
+ if (occe){
+ wire_maker.Add(occe->getTopoDS_Edge());
+ }
+ }
+ TopoDS_Wire myWire = wire_maker.Wire();
+
+
+ GEntity *ret;
+ if (base->cast2Vertex()){
+ OCCVertex *occv = dynamic_cast<OCCVertex*>(base);
+ BRepOffsetAPI_MakePipe myNiceLittlePipe (myWire,occv->getShape());
+ TopoDS_Edge result = TopoDS::Edge(myNiceLittlePipe.Shape());
+ ret = gm->_occ_internals->addEdgeToModel(gm, result);
+ }
+ if (base->cast2Edge()){
+ OCCEdge *occe = dynamic_cast<OCCEdge*>(base);
+ BRepOffsetAPI_MakePipe myNiceLittlePipe (myWire,occe->getTopoDS_Edge());
+ TopoDS_Face result = TopoDS::Face(myNiceLittlePipe.Shape());
+ ret = gm->_occ_internals->addFaceToModel(gm, result);
+ }
+ if (base->cast2Face()){
+ OCCFace *occf = dynamic_cast<OCCFace*>(base);
+ BRepOffsetAPI_MakePipe myNiceLittlePipe (myWire,occf->getTopoDS_Face());
+ TopoDS_Solid result = TopoDS::Solid(myNiceLittlePipe.Shape());
+ ret = gm->_occ_internals->addRegionToModel(gm, result);
+ }
+ return ret;
+}
+
+
+
+
#endif
diff --git a/Geo/GModelFactory.h b/Geo/GModelFactory.h
index cad242d024..ae3f340e0e 100644
--- a/Geo/GModelFactory.h
+++ b/Geo/GModelFactory.h
@@ -25,8 +25,10 @@ class GModelFactory {
// brep primitives
enum arcCreationMethod {THREE_POINTS=1, CENTER_START_END=2};
enum splineType {BEZIER=1, BSPLINE=2};
+ // vertex primitive
virtual GVertex *addVertex(GModel *gm, double x, double y, double z,
double lc) = 0;
+ // edge primitives
virtual GEdge *addLine(GModel *, GVertex *v1, GVertex *v2) = 0;
virtual GEdge *addCircleArc(GModel *gm, const arcCreationMethod &method,
GVertex *start, GVertex *end,
@@ -40,10 +42,19 @@ class GModelFactory {
std::vector<double> knots,
std::vector<double> weights,
std::vector<int> multiplicity) = 0;
+ // faces primitives
+ // this one tries to build a model face with one single list
+ // of faces. If boundaries are co-planar, then it's a plane,
+ // otherwise, we tru ruled, sweep or other kind of surfaces
+ virtual std::vector<GFace *> addRuledFaces (GModel *gm, std::vector<std::vector<GEdge *> > edges) = 0;
+ virtual GFace * addFace (GModel *gm, std::vector<GEdge *> edges, std::vector< std::vector<double > > points) = 0;
+ virtual GFace * addPlanarFace (GModel *gm, std::vector<std::vector<GEdge *> > edges) = 0;
+ // sweep stuff
virtual GEntity *revolve(GModel *gm, GEntity*, std::vector<double> p1,
std::vector<double> p2, double angle) = 0;
virtual GEntity *extrude(GModel *gm, GEntity*, std::vector<double> p1,
std::vector<double> p2) = 0;
+ virtual GEntity* addPipe (GModel *gm, GEntity *base, std::vector<GEdge *> wire) = 0;
// solid primitives
virtual GEntity *addSphere(GModel *gm, double cx, double cy, double cz,
@@ -99,9 +110,13 @@ class OCCFactory : public GModelFactory {
std::vector<double> p2, double angle);
GEntity *extrude (GModel *gm, GEntity*,std::vector<double> p1,
std::vector<double> p2);
+ GEntity* addPipe (GModel *gm, GEntity *base, std::vector<GEdge *> wire);
GEntity *addSphere(GModel *gm,double cx, double cy, double cz, double radius);
GEntity *addCylinder(GModel *gm,std::vector<double> p1, std::vector<double> p2,
double radius);
+ std::vector<GFace *> addRuledFaces (GModel *gm, std::vector<std::vector<GEdge *> > edges);
+ GFace * addFace (GModel *gm, std::vector<GEdge *> edges, std::vector< std::vector<double > > points);
+ GFace * addPlanarFace (GModel *gm, std::vector<std::vector<GEdge *> > edges);
GEntity *addTorus(GModel *gm,std::vector<double> p1, std::vector<double> p2,
double radius1, double radius2);
GEntity *addBlock(GModel *gm,std::vector<double> p1, std::vector<double> p2);
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 05dac36895..2072897494 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -36,7 +36,6 @@ void OCC_Internals::buildShapeFromLists(TopoDS_Shape _shape)
BRep_Builder B;
TopoDS_Compound C;
B.MakeCompound(C);
- // B.Add(C,_shape);
TopTools_ListOfShape theList;
addSimpleShapes(_shape, theList);
@@ -45,9 +44,12 @@ void OCC_Internals::buildShapeFromLists(TopoDS_Shape _shape)
for(int i = 1; i <= vmap.Extent(); i++) B.Add(C, vmap(i));
for(int i = 1; i <= emap.Extent(); i++) B.Add(C, emap(i));
+ for(int i = 1; i <= wmap.Extent(); i++) B.Add(C, wmap(i));
for(int i = 1; i <= fmap.Extent(); i++) B.Add(C, fmap(i));
+ for(int i = 1; i <= shmap.Extent(); i++) B.Add(C, shmap(i));
for(int i = 1; i <= somap.Extent(); i++) B.Add(C, somap(i));
shape = C;
+
}
const TopoDS_Shape *OCC_Internals::lookupInLists(TopoDS_Shape _shape)
@@ -413,7 +415,15 @@ void OCC_Internals::loadBREP(const char *fn)
void OCC_Internals::writeBREP(const char *fn)
{
- BRepTools::Write(shape, (char*)fn);
+ std::ofstream myFile;
+ myFile.open (fn);
+ try {
+ BRepTools::Write(shape, myFile);
+ }
+ catch(Standard_Failure &err){
+ Msg::Error("%s", err.GetMessageString());
+ }
+ myFile.close ();
}
void OCC_Internals::loadSTEP(const char *fn)
@@ -472,6 +482,8 @@ GVertex *OCC_Internals::addVertexToModel(GModel *model, TopoDS_Vertex vertex)
GVertex *gv = getOCCVertexByNativePtr(model, vertex);
if (gv) return gv;
addShapeToLists(vertex);
+ buildShapeFromLists(vertex);
+ // buildLists();
buildGModel(model);
return getOCCVertexByNativePtr (model,vertex);
}
@@ -481,6 +493,8 @@ GEdge *OCC_Internals::addEdgeToModel(GModel *model, TopoDS_Edge edge)
GEdge *ge = getOCCEdgeByNativePtr(model, edge);
if (ge) return ge;
addShapeToLists(edge);
+ buildShapeFromLists(edge);
+ // buildLists();
buildGModel(model);
return getOCCEdgeByNativePtr(model,edge);
}
@@ -489,6 +503,8 @@ GFace* OCC_Internals::addFaceToModel(GModel *model, TopoDS_Face face){
GFace *gf = getOCCFaceByNativePtr(model,face);
if (gf) return gf;
addShapeToLists(face);
+ buildShapeFromLists(face);
+ // buildLists();
buildGModel(model);
return getOCCFaceByNativePtr(model,face);
}
@@ -498,6 +514,8 @@ GRegion* OCC_Internals::addRegionToModel(GModel *model, TopoDS_Solid region){
GRegion *gr = getOCCRegionByNativePtr(model,region);
if (gr) return gr;
addShapeToLists(region);
+ buildShapeFromLists(region);
+ //buildLists();
// buildLists();
// TEST
buildGModel(model);
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index 9bd1539211..fa0c64346b 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -265,6 +265,11 @@ bool GRegion::edgeConnected(GRegion *r) const
void GRegion::replaceFaces (std::list<GFace*> &new_faces)
{
replaceFacesInternal (new_faces);
+ if (l_faces.size() != new_faces.size()){
+ Msg::Fatal("impossible to replace faces in region %d (%d vs %d)",tag(),
+ l_faces.size(),new_faces.size());
+ }
+
std::list<GFace*>::iterator it = l_faces.begin();
std::list<GFace*>::iterator it2 = new_faces.begin();
std::list<int>::iterator it3 = l_dirs.begin();
diff --git a/Geo/OCCRegion.cpp b/Geo/OCCRegion.cpp
index 946da65ffe..5a1dff017d 100644
--- a/Geo/OCCRegion.cpp
+++ b/Geo/OCCRegion.cpp
@@ -23,6 +23,7 @@ OCCRegion::OCCRegion(GModel *m, TopoDS_Solid _s, int num)
void OCCRegion::setup()
{
+ l_faces.clear();
TopExp_Explorer exp2, exp3;
for(exp2.Init(s, TopAbs_SHELL); exp2.More(); exp2.Next()){
TopoDS_Shape shell = exp2.Current();
@@ -57,8 +58,6 @@ bool FaceHaveDifferentOrientations(const TopoDS_Face& aFR,
void OCCRegion::replaceFacesInternal(std::list<GFace*> &new_faces)
{
// we simply replace old faces by new faces in the structure
- Standard_Integer aNbS;
- TopTools_IndexedMapOfShape aMS;
TopExp_Explorer aExpS, aExpF;
BRep_Builder aBB;
TopoDS_Compound aCmp;
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index a7f3f43e68..ac501538e7 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -38,6 +38,32 @@ template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,Fu
}
}
+template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,
+ FunctionSpaceBase &shapeFcts,
+ FunctionSpaceBase &testFcts,
+ Iterator itbegin,
+ Iterator itend,
+ QuadratureBase &integrator,
+ Assembler &assembler) // non symmetric
+{
+ fullMatrix<typename Assembler::dataMat> localMatrix;
+ std::vector<Dof> R;
+ std::vector<Dof> C;
+ for (Iterator it = itbegin;it!=itend; ++it)
+ {
+ MElement *e = *it;
+ R.clear();
+ C.clear();
+ IntPt *GP;
+ int npts=integrator.getIntPoints(e,&GP);
+ term.get(e,npts,GP,localMatrix);
+ shapeFcts.getKeys(e,R);
+ testFcts.getKeys(e,C);
+ assembler.assemble(R, C, localMatrix);
+ }
+}
+
+
template<class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,MElement *e,QuadratureBase &integrator,Assembler &assembler) // symmetric
{
fullMatrix<typename Assembler::dataMat> localMatrix;
diff --git a/benchmarks/boolean/square1.lua b/benchmarks/boolean/square1.lua
index e9619fc833..a5a277a971 100644
--- a/benchmarks/boolean/square1.lua
+++ b/benchmarks/boolean/square1.lua
@@ -1,10 +1,36 @@
g = GModel()
-v1 = g:addVertex(0, 0, 0, 1)
-v2 = g:addVertex(1, 0, 0, 1)
-v3 = g:addVertex(1, 1, 0, 1)
-v4 = g:addVertex(0, 1, 0, 1)
+v1 = g:addVertex(0, 0, 0, .1)
+v2 = g:addVertex(1, 0, 0, .1)
+v3 = g:addVertex(1, 1, 0, .1)
+v4 = g:addVertex(0, 1, 0, .1)
e1 = g:addLine(v1, v2)
e2 = g:addLine(v2, v3)
e3 = g:addLine(v3, v4)
e4 = g:addLine(v4, v1)
+v11 = g:addVertex(.4, .4, 0, .1)
+v12 = g:addVertex(.5, .4, 0, .1)
+v13 = g:addVertex(.5, .5, 0, .1)
+v14 = g:addVertex(.4, .5, 0, .1)
+e11 = g:addLine(v11, v12)
+e12 = g:addLine(v12, v13)
+e13 = g:addLine(v13, v14)
+e14 = g:addLine(v14, v11)
+f = g:addPlanarFace ({{e1,e2,e3,e4},{e11,e12,e13,e14}})
+
+v100 = g:addVertex(0, 0, 0, .1)
+v200 = g:addVertex(0, 0, 1, .1)
+v300 = g:addVertex(0, .1, .33, .1)
+v400 = g:addVertex(.1, .1, .66, .1)
+line = g:addBezier(v100,v200,{{v300:x(),v300:y(),v300:z()},{v400:x(),v400:y(),v400:z()}});
+g:addPipe (f, {line})
+
+g:glue(1.e-9);
+
+myTool = GModel();
+myTool:addSphere(0.2,0.2,0.1,.52012);
+
+--g:addSphere(1,1.3,1,.3);
+g:computeDifference(myTool,0);
+
+g:setAsCurrent();
--
GitLab