diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 4f8685ff42167732df7a0a5118931e7b7611ad52..514533cc4e539888164f5efaf9042d15a863b2c2 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -10,7 +10,8 @@ set(SRC
GVertex.cpp GEdge.cpp GFace.cpp GRegion.cpp
GEdgeLoop.cpp GEdgeCompound.cpp GFaceCompound.cpp
GRegionCompound.cpp GRbf.cpp
- gmshVertex.cpp gmshEdge.cpp gmshFace.cpp gmshRegion.cpp gmshSurface.cpp
+ gmshVertex.cpp gmshEdge.cpp gmshFace.cpp gmshRegion.cpp
+ gmshCurve.cpp gmshSurface.cpp
OCCVertex.cpp OCCEdge.cpp OCCFace.cpp OCCRegion.cpp
discreteEdge.cpp discreteFace.cpp discreteRegion.cpp
fourierEdge.cpp fourierFace.cpp fourierProjectionFace.cpp
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 9f1620c42c4767d1f9236568eb0062fad2c2b448..238c63b82b0ba06cffb37b9dc35c7e5b06eebf76 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -967,7 +967,6 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[
GPoint pnt = point(initialGuess[0], initialGuess[1]);
SPoint3 spnt(pnt.x(), pnt.y(), pnt.z());
double dist = queryPoint.distance(spnt);
- printf("dist = %12.5E\n",dist);
}
// FILE *F = Fopen ("hop.pos","w");
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 096918d28027447ca7a0511c908eefb6732ab34f..d2e04c5aebe241c1dc59ac3abdbc57829f2152b8 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2652,6 +2652,15 @@ GEdge *GModel::addBezier(GVertex *start, GVertex *end,
return 0;
}
+GEdge *GModel::addBSpline(GVertex *start, GVertex *end,
+ std::vector<std::vector<double> > points)
+{
+ if(_factory)
+ return _factory->addSpline(this, GModelFactory::BSPLINE, start, end,
+ points);
+ return 0;
+}
+
GEdge *GModel::addNURBS(GVertex *start, GVertex *end,
std::vector<std::vector<double> > points,
std::vector<double> knots,
diff --git a/Geo/GModel.h b/Geo/GModel.h
index f026f38b6ce8ee25f0d002eeffc46a0a7abf6339..8fc5068a07e53c392aca034c9fd9992ecabf6486 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -489,6 +489,7 @@ class GModel
GEdge *addCircleArcCenter(GVertex *start, GVertex *center, GVertex *end);
GEdge *addCircleArc3Points(double x, double y, double z, GVertex *start, GVertex *end);
GEdge *addBezier(GVertex *start, GVertex *end, std::vector<std::vector<double> > points);
+ GEdge *addBSpline(GVertex *start, GVertex *end, std::vector<std::vector<double> > points);
GEdge *addNURBS(GVertex *start, GVertex *end,
std::vector<std::vector<double> > points,
std::vector<double> knots,
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 19b29e28063707280589e51142c591ed82b5dafb..4a7cbeae41e38b6269649a5b54553c73b10b6482 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -20,6 +20,7 @@
#include "Geo.h"
#include "GmshDefines.h"
+
GVertex *GeoFactory::addVertex(GModel *gm, double x, double y, double z, double lc)
{
int num = gm->getMaxElementaryNumber(0) + 1;
@@ -481,6 +482,7 @@ std::vector<GEntity*> GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e,
#include <TopoDS_Wire.hxx>
#include <TopoDS_Compound.hxx>
#include <TopTools_ListOfShape.hxx>
+#include <GeomAPI_PointsToBSpline.hxx>
#include "OCC_Connect.h"
#if defined(HAVE_SALOME)
#include "Partition_Spliter.hxx"
@@ -581,18 +583,33 @@ GEdge *OCCFactory::addSpline(GModel *gm, const splineType &type,
TColgp_Array1OfPnt ctrlPoints(1, nbControlPoints + 2);
int index = 1;
ctrlPoints.SetValue(index++, gp_Pnt(start->x(), start->y(), start->z()));
+
+ // printf("%d %d %d %d\n",points.size(),points[0].size(),points[1].size(),points[2].size());
for (int i = 0; i < nbControlPoints; i++) {
+ // printf("%g %g %g\n",points[i][0],points[i][1],points[i][2]);
gp_Pnt aP(points[i][0],points[i][1],points[i][2]);
ctrlPoints.SetValue(index++, aP);
}
ctrlPoints.SetValue(index++, gp_Pnt(end->x(), end->y(), end->z()));
if (type == BEZIER) {
+ if (nbControlPoints >= 20)Msg::Fatal("OCC Bezier Curves have a maximum degree of 20");
Handle(Geom_BezierCurve) Bez = new Geom_BezierCurve(ctrlPoints);
if (occv1 && occv2)
occEdge = BRepBuilderAPI_MakeEdge(Bez,occv1->getShape(),occv2->getShape()).Edge();
else
occEdge = BRepBuilderAPI_MakeEdge(Bez).Edge();
}
+ else if (type == BSPLINE) {
+
+ Handle(Geom_BSplineCurve) Bez = GeomAPI_PointsToBSpline(ctrlPoints).Curve();
+
+ if (occv1 && occv2)
+ occEdge = BRepBuilderAPI_MakeEdge(Bez,occv1->getShape(),occv2->getShape()).Edge();
+ else
+ occEdge = BRepBuilderAPI_MakeEdge(Bez).Edge();
+ }
+
+
return gm->_occ_internals->addEdgeToModel(gm, occEdge);
}
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 42bb88ce9eba4c0d9893b2c78b1d564b5a034a9e..f6b8a077f84602ee40bf8fe105e44de3fef54bcf 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -692,10 +692,17 @@ GRegion* OCC_Internals::addRegionToModel(GModel *model, TopoDS_Solid region)
{
GRegion *gr = getOCCRegionByNativePtr(model, region);
if(gr) return gr;
- addShapeToLists(region);
+
buildShapeFromLists(region);
+ model->destroy();
+ buildLists();
buildGModel(model);
return getOCCRegionByNativePtr(model, region);
+
+ // addShapeToLists(region);
+ // buildShapeFromLists(region);
+ // buildGModel(model);
+ // return getOCCRegionByNativePtr(model, region);
}
/* I needed getGTagOfOCC*ByNativePtr whithin setPhysicalNumToEntitiesInBox */
diff --git a/Geo/MElementOctree.cpp b/Geo/MElementOctree.cpp
index ce8db23964fa307ab2ca2b4611364a0d605b6aa0..e65860ca16fbcb8d1e03650ceaa1815c5438cc66 100644
--- a/Geo/MElementOctree.cpp
+++ b/Geo/MElementOctree.cpp
@@ -9,7 +9,7 @@
#include "Octree.h"
#include "Context.h"
-static void MElementBB(void *a, double *min, double *max)
+void MElementBB(void *a, double *min, double *max)
{
MElement *e = (MElement*)a;
MVertex *v = e->getVertex(0);
@@ -47,7 +47,7 @@ static void MElementCentroid(void *a, double *x)
x[2] *= oc;
}
-static int MElementInEle(void *a, double *x)
+int MElementInEle(void *a, double *x)
{
MElement *e = (MElement*)a;
double uvw[3];
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 40fff3e7cf056e791050c2e8d00ebf1f19be8c3b..98dceccb732b641d583aa304098376eb85204618 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -452,17 +452,6 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
}
}
}
- /*
- if (p1.size() == 8 ||p2.size() == 8){
- for (int i=0;i<p1.size();i++){
- printf("p1[%d] = %g %g\n",i,p1[i].x(),p1[i].y());
- }
- for (int i=0;i<p2.size();i++){
- printf("p2[%d] = %g %g\n",i,p2[i].x(),p2[i].y());
- }
- printf("%d %d\n",imin,jmin);
- }
- */
param1 = p1[jmin];
param2 = p2[imin];
}
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 589d9a1668c2dc29322e178fc49e7617c963db54..23c59eba933d3fe5e5b4c72860df3653c75b7a22 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -93,18 +93,20 @@ SPoint2 OCCEdge::reparamOnFace(const GFace *face, double epar, int dir) const
pnt.Coord(u, v);
// sometimes OCC miserably fails ...
- GPoint p1 = point(epar);
- GPoint p2 = face->point(u, v);
- const double dx = p1.x()-p2.x();
- const double dy = p1.y()-p2.y();
- const double dz = p1.z()-p2.z();
- if(sqrt(dx * dx + dy * dy + dz * dz) > 1.e-2 * CTX::instance()->lc){
- Msg::Warning("Reparam on face was inaccurate for curve %d on surface %d at point %g",
- tag(), face->tag(), epar);
- Msg::Warning("On the face %d local (%g %g) global (%g %g %g)",
- face->tag(), u, v, p2.x(), p2.y(), p2.z());
- Msg::Warning("On the edge %d local (%g) global (%g %g %g)",
- tag(), epar, p1.x(), p1.y(), p1.z());
+ if (0){
+ GPoint p1 = point(epar);
+ GPoint p2 = face->point(u, v);
+ const double dx = p1.x()-p2.x();
+ const double dy = p1.y()-p2.y();
+ const double dz = p1.z()-p2.z();
+ if(sqrt(dx * dx + dy * dy + dz * dz) > 1.e-2 * CTX::instance()->lc){
+ Msg::Warning("Reparam on face was inaccurate for curve %d on surface %d at point %g",
+ tag(), face->tag(), epar);
+ Msg::Warning("On the face %d local (%g %g) global (%g %g %g)",
+ face->tag(), u, v, p2.x(), p2.y(), p2.z());
+ Msg::Warning("On the edge %d local (%g) global (%g %g %g)",
+ tag(), epar, p1.x(), p1.y(), p1.z());
+ }
}
return SPoint2(u, v);
}
diff --git a/Geo/boundaryLayersData.cpp b/Geo/boundaryLayersData.cpp
index 7d572a276ce97ad5d169396d3218e5474a28fb70..c6430d92619e39caebe8570a8b0d956ff82c46e9 100644
--- a/Geo/boundaryLayersData.cpp
+++ b/Geo/boundaryLayersData.cpp
@@ -17,10 +17,6 @@
#include "OS.h"
#include "BackgroundMesh.h"
-#if defined(HAVE_RTREE)
-#include "rtree.h"
-#endif
-
#if !defined(HAVE_MESH) || !defined(HAVE_ANN)
BoundaryLayerField* getBLField(GModel *gm){ return 0; }
@@ -152,23 +148,10 @@ const BoundaryLayerData & BoundaryLayerColumns::getColumn(MVertex *v, MFace f)
{
int N = getNbColumns(v) ;
if (N == 1) return getColumn(v, 0);
- if (isOnWedge (v)){
- GFace *gf = _inverse_classification[f];
- BoundaryLayerFanWedge3d w = getWedge(v);
- if (w.isLeft(gf))return getColumn(v, 0);
- if (w.isRight(gf))return getColumn(v, N-1);
- Msg::Error("Strange behavior for a wedge");
- }
- if (isCorner (v)){
- GFace *gf = _inverse_classification[f];
- BoundaryLayerFanCorner3d c = getCorner(v);
- int k = 0;
- for (unsigned int i=0;i<c._gf.size();i++){
- if (c._gf[i] == gf){
- return getColumn(v, k);
- }
- k += (c._fanSize - 1);
- }
+ GFace *gf = _inverse_classification[f];
+ for (int i=0;i<N;i++){
+ const BoundaryLayerData & c = getColumn(v, i);
+ if (std::find(c._joint.begin(),c._joint.end(),gf) != c._joint.end())return c;
}
static BoundaryLayerData error;
return error;
@@ -782,8 +765,6 @@ bool buildAdditionalPoints2D(GFace *gf)
}
// DEBUG STUFF
- _columns->filterPoints(gf,0.21);
-
char name[256];
sprintf(name,"points_face_%d.pos",gf->tag());
FILE *f = Fopen (name,"w");
@@ -835,6 +816,56 @@ static void createBLPointsInDir(GRegion *gr,
}
}
+
+static bool createWedgeBetweenTwoFaces(MVertex *myV,
+ SVector3 n1, SVector3 n2,
+ std::vector<SVector3> &shoot)
+{
+ double angle = angle_0_180 (n1,n2);
+ int fanSize = FANSIZE__; //angle / _treshold;
+ for (int i=-1; i<=fanSize; i++){
+
+ double ti = (double)(i+1)/ (fanSize+1);
+ double angle_t = ti * angle;
+ double cosA = cos (angle_t);
+ double cosAlpha = dot(n1,n2);
+
+ const double A = (1.- 2.*cosA*cosA) + cosAlpha*cosAlpha - 2 * cosAlpha*(1.-cosA*cosA);
+ const double B = -2*(1.-cosA*cosA)*(1-cosAlpha);
+ const double C = 1.-cosA*cosA;
+ double DELTA = B*B-4*A*C;
+ double t = 0.0;
+ if (A == 0.0){
+ t = -C / B;
+ }
+ else if (C != 0.0){
+ if (DELTA < 0){
+ Msg::Error("this should not happen DELTA = %12.5E",DELTA);
+ DELTA = 0;
+ }
+ const double t1 = (-B+sqrt(DELTA))/(2.*A);
+ const double t2 = (-B-sqrt(DELTA))/(2.*A);
+
+ SVector3 x1 (n1*(1.-t1) + n2 * t2);
+ SVector3 x2 (n1*(1.-t2) + n2 * t2);
+ double a1 = angle_0_180 (n1,x1);
+ double a2 = angle_0_180 (n1,x2);
+ if (fabs(a1 - angle_t) < fabs(a2 - angle_t))t = t1;
+ else t = t2;
+ }
+ SVector3 x (n1*(1.-t) + n2 * t);
+ x.normalize();
+ shoot.push_back(x);
+ }
+ return true;
+}
+
+void computeAllGEdgesThatAreCreatingFans (GRegion *gr, std::vector<GEdge*> &fans)
+{
+
+}
+
+
static void createColumnsBetweenFaces(GRegion *gr,
MVertex *myV,
BoundaryLayerField *blf,
@@ -870,6 +901,7 @@ static void createColumnsBetweenFaces(GRegion *gr,
// printf("vertex %d %d faces\n",myV->getNum(),count);
std::set<int> done;
+ std::vector< std::vector<GFace*> > joints;
for (int i=0;i<count;i++){
if (done.find(i) == done.end()){
SVector3 n1 = n[i];
@@ -893,91 +925,38 @@ static void createColumnsBetweenFaces(GRegion *gr,
if (joint.size()){
std::vector<MVertex*> _column;
std::vector<SMetric3> _metrics;
+ joints.push_back(joint);
avg.normalize();
createBLPointsInDir (gr,myV,blf,avg,_column,_metrics);
_columns->addColumn(avg,myV, _column, _metrics, joint);
}
- // printf("adding one column for %d faces\n",joint.size());
+ // printf("adding one column for %d fqaces\n",joint.size());
}
}
-}
-
-/*
-static bool createWedgeBetweenTwoFaces(bool testOnly,
- MVertex *myV,
- GFace *gf1, GFace *gf2,
- std::multimap<GFace*,MTriangle*> & _faces,
- std::map<MFace,SVector3,Less_Face> &_normals,
- double _treshold,
- std::vector<SVector3> &shoot)
-{
- SVector3 n1,n2;
- SPoint3 c1,c2;
- for (std::multimap<GFace*,MTriangle*>::iterator itm =
- _faces.lower_bound(gf1);
- itm != _faces.upper_bound(gf1); ++itm){
- n1 += _normals[itm->second->getFace(0)];
- c1 = itm->second->getFace(0).barycenter();
- }
- for (std::multimap<GFace*,MTriangle*>::iterator itm =
- _faces.lower_bound(gf2);
- itm != _faces.upper_bound(gf2); ++itm){
- n2 += _normals[itm->second->getFace(0)];
- c2 = itm->second->getFace(0).barycenter();
- }
- n1.normalize();
- n2.normalize();
- // FIXME WRONG FOR INTERNAL CORNERS !!!
- double angle = angle_0_180 (n1,n2);
- double sign = dot((n1-n2),(c1-c2));
- if (angle > _treshold && sign > 0){
- if(testOnly)return true;
- int fanSize = FANSIZE__; //angle / _treshold;
- for (int i=-1; i<=fanSize; i++){
-
- double ti = (double)(i+1)/ (fanSize+1);
- double angle_t = ti * angle;
- double cosA = cos (angle_t);
- double cosAlpha = dot(n1,n2);
-
- const double A = (1.- 2.*cosA*cosA) + cosAlpha*cosAlpha - 2 * cosAlpha*(1.-cosA*cosA);
- const double B = -2*(1.-cosA*cosA)*(1-cosAlpha);
- const double C = 1.-cosA*cosA;
- double DELTA = B*B-4*A*C;
- double t = 0.0;
- if (A == 0.0){
- t = -C / B;
- }
- else if (C != 0.0){
- if (DELTA < 0){
- Msg::Error("this should not happen DELTA = %12.5E",DELTA);
- DELTA = 0;
+ // create wedges
+ if (joints.size() > 1){
+ for (unsigned int I = 0 ; I < joints.size() ; I++){
+ const BoundaryLayerData & c0 = _columns->getColumn(myV, I);
+ for (unsigned int J = I+1 ; J < joints.size() ; J++){
+ const BoundaryLayerData & c1 = _columns->getColumn(myV, J);
+ std::vector<SVector3> shoot;
+ createWedgeBetweenTwoFaces(myV,c0._n,c1._n,shoot);
+ for (unsigned int i=1;i<shoot.size()-1;i++){
+ std::vector<MVertex*> _column;
+ std::vector<SMetric3> _metrics;
+ createBLPointsInDir (gr,myV,blf,shoot[i],_column,_metrics);
+ _columns->addColumn(shoot[i] , myV, _column, _metrics);
}
- const double t1 = (-B+sqrt(DELTA))/(2.*A);
- const double t2 = (-B-sqrt(DELTA))/(2.*A);
-
- SVector3 x1 (n1*(1.-t1) + n2 * t2);
- SVector3 x2 (n1*(1.-t2) + n2 * t2);
- double a1 = angle_0_180 (n1,x1);
- double a2 = angle_0_180 (n1,x2);
- if (fabs(a1 - angle_t) < fabs(a2 - angle_t))t = t1;
- else t = t2;
}
- SVector3 x (n1*(1.-t) + n2 * t);
- x.normalize();
- shoot.push_back(x);
}
- return true;
}
- else {
- if(testOnly)return false;
- SVector3 n = n1+n2;
- n.normalize();
- shoot.push_back(n);
- return false;
+ // we have a corner : in a first step, only add one single
+ // in the average direction
+ if (joints.size() > 2){
}
}
-*/
+
+
BoundaryLayerColumns *buildAdditionalPoints3D(GRegion *gr)
{
@@ -1079,6 +1058,8 @@ BoundaryLayerColumns *buildAdditionalPoints3D(GRegion *gr)
}
}
+ // we have a 3D boundary layer that connects a 2D boundary layer
+ // this is the tricky case
if (onSymmetryPlane){
for ( std::list<GFace*>::iterator itf = faces.begin(); itf!= faces.end() ; ++itf){
BoundaryLayerColumns* _face_columns = (*itf)->getColumns();
@@ -1087,10 +1068,70 @@ BoundaryLayerColumns *buildAdditionalPoints3D(GRegion *gr)
std::vector<GFace*> _joint;
_joint.insert(_joint.begin(),_allGFaces.begin(),_allGFaces.end());
const BoundaryLayerData & c = _face_columns->getColumn(*it,0);
- _columns->addColumn(_allDirs[0],*it, c._column, c._metrics, _joint);
+ _columns->addColumn(c._n,*it, c._column, c._metrics, _joint);
}
else if (N > 1){
- Msg::Error ("Impossible connexion between face and region BLs");
+ if (_allGFaces.size() != 2){
+ Msg::Fatal("cannot solve such a strange stuff in the BL");
+ }
+ // printf("%d columns\n",N);
+ std::set<GFace*>::iterator itff = _allGFaces.begin();
+ GFace *g1 = *itff ; ++itff; GFace *g2 = *itff;
+ int sense = 1;
+ std::vector<GFace*> _joint;
+
+ const BoundaryLayerFan *fan = _face_columns->getFan(*it);
+
+ if (fan){
+ MVertex *v11 = fan->_e1.getVertex(0);
+ MVertex *v12 = fan->_e1.getVertex(1);
+ std::list<GEdge*> l1 = g1->edges();
+ std::list<GEdge*> l2 = g2->edges();
+ if (v11 == *it){
+ if (v12->onWhat()->dim() == 1){
+ GEdge *ge1 = (GEdge*)v12->onWhat();
+ // printf("COUCOU %d %d %d\n",fan->sense,std::find(l1.begin(),l1.end(),ge1) != l1.end(),std::find(l2.begin(),l2.end(),ge1) != l2.end());
+ if (std::find(l1.begin(),l1.end(),ge1) != l1.end())sense = fan->sense;
+ else if (std::find(l2.begin(),l2.end(),ge1) != l2.end())sense = -fan->sense;
+ else printf("strange1 %d %d \n");
+ }
+ else Msg::Error("Cannot choose between directions in a BL (dim = %d)",v12->onWhat()->dim());
+ }
+ else {
+ if (v11->onWhat()->dim() == 1){
+ GEdge *ge1 = (GEdge*)v11->onWhat();
+ if (std::find(l1.begin(),l1.end(),ge1) != l1.end())sense = fan->sense;
+ else if (std::find(l2.begin(),l2.end(),ge1) != l2.end())sense = -fan->sense;
+ else printf("strange2 %d %d \n");
+ }
+ else Msg::Error("Cannot choose between directions in a BL");
+ }
+ }
+ else{
+ Msg::Error("No fan on the outgoing BL");
+ }
+ _joint.push_back(g1);
+ const BoundaryLayerData & c0 = _face_columns->getColumn(*it,sense==1 ? 0 : N-1);
+ _columns->addColumn(c0._n,*it, c0._column, c0._metrics,_joint);
+ _joint.clear();
+ _joint.push_back(g2);
+ const BoundaryLayerData & cN = _face_columns->getColumn(*it,sense==1 ? N-1 : 0);
+ _columns->addColumn(cN._n,*it, cN._column, cN._metrics,_joint);
+ // printf("%g %g %g --> %g %g %g\n",c0._n.x(),c0._n.y(),c0._n.z(),cN._n.x(),cN._n.y(),cN._n.z());
+ if (sense==1){
+ for (int k=1;k<N-1;k++){
+ const BoundaryLayerData & c = _face_columns->getColumn(*it,k);
+ _columns->addColumn(c._n,*it, c._column, c._metrics);
+ // printf("%g %g %g \n",c._n.x(),c._n.y(),c._n.z());
+ }
+ }
+ else {
+ for (int k=N-2;k>0;k--){
+ const BoundaryLayerData & c = _face_columns->getColumn(*it,k);
+ _columns->addColumn(c._n,*it, c._column, c._metrics);
+ // printf("%g %g %g \n",c._n.x(),c._n.y(),c._n.z());
+ }
+ }
}
}
}
@@ -1101,7 +1142,7 @@ BoundaryLayerColumns *buildAdditionalPoints3D(GRegion *gr)
// DEBUG STUFF
- FILE *f = fopen ("test3D.pos","w");
+ FILE *f = fopen ("POINTS3D.pos","w");
fprintf(f,"View \"\" {\n");
for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){
MVertex *v = *it;
@@ -1121,176 +1162,4 @@ BoundaryLayerColumns *buildAdditionalPoints3D(GRegion *gr)
return _columns;
}
-struct blPoint_wrapper
-{
- bool _tooclose;
- MVertex *_v;
- std::map<MVertex*,MVertex*> &_v2v;
- blPoint_wrapper (MVertex *v, std::map<MVertex*,MVertex*> &v2v)
- : _tooclose(false), _v(v), _v2v(v2v) {}
-};
-
-struct blPoint_rtree
-{
- MVertex *_v;
- double _size;
- blPoint_rtree (MVertex *v, double size) :
- _v(v), _size(size) {}
- bool inExclusionZone (MVertex *v){
- double d = _v->distance(v);
- //printf("d = %12.5E\n",d);
- if (d <= _size) return true;
- return false;
- }
- void minmax (double min[3], double max[3]){
- min[0] = _v->x() - _size;
- min[1] = _v->y() - _size;
- min[2] = _v->z() - _size;
- max[0] = _v->x() + _size;
- max[1] = _v->y() + _size;
- max[2] = _v->z() + _size;
- }
-};
-
-
-bool rtree_callback(blPoint_rtree *neighbour,void* point){
- blPoint_wrapper *w = static_cast<blPoint_wrapper*>(point);
-
- const MVertex *from_1 = w->_v2v[neighbour->_v];
- const MVertex *from_2 = w->_v2v[w->_v];
-
- // printf("%p %p\n",from_1,from_2);
-
- if (from_1 == from_2) {
- return true;
- }
-
- if (neighbour->inExclusionZone(w->_v)){
- w->_tooclose = true;
- return false;
- }
- return true;
-}
-
-#if defined(HAVE_RTREE)
-bool inExclusionZone_filter (MVertex* p,
- std::map <MVertex*, MVertex*> &v2v,
- RTree< blPoint_rtree *,double,3,double> &rtree){
- // should assert that the point is inside the domain
- {
- double u, v;
- p->getParameter(0,u);
- p->getParameter(1,v);
- if (!backgroundMesh::current()->inDomain(u,v,0)) return true;
- }
-
- blPoint_wrapper w (p,v2v);
- double _min[3] = {p->x()-1.e-1, p->y()-1.e-1,p->z()-1.e-1};
- double _max[3] = {p->x()+1.e-1, p->y()+1.e-1,p->z()+1.e-1};
- rtree.Search(_min,_max,rtree_callback,&w);
-
- return w._tooclose;
-}
-#endif
-
-
-void BoundaryLayerColumns::filterPoints(GEntity *ge, double factor)
-{
-#if defined(HAVE_RTREE)
- // return;
- // compute the element sizes
- std::map<MVertex*,double> sizes;
- if (ge->dim() == 2){
- backgroundMesh::set((GFace *)ge);
- std::list<GEdge*> edges = ge->edges();
- std::list<GEdge*>::iterator it = edges.begin();
- for ( ; it != edges.end() ; ++it){
- GEdge *ged = *it;
- for (unsigned int i=0;i<ged->lines.size();i++){
- MLine *e = ged->lines[i];
- MVertex *v0 = e->getVertex(0);
- MVertex *v1 = e->getVertex(1);
- double d = v0->distance(v1);
- std::map<MVertex*,double>::iterator it0 = sizes.find(v0);
- if (it0 == sizes.end()) sizes[v0] = d;
- else it0->second = std::max(d, it0->second);
- std::map<MVertex*,double>::iterator it1 = sizes.find(v1);
- if (it1 == sizes.end()) sizes[v1] = d;
- else it1->second = std::max(d, it1->second);
- }
- }
- }
- else {
- Msg::Fatal("code ce truc JF !");
- }
-
- // a RTREE data structure that enables to verify if
- // points are too close
- RTree<blPoint_rtree*,double,3,double> rtree;
- // stores the info "where the new vertex comes form"
- std::map <MVertex*, MVertex*> v2v;
-
- // compute maximum column size
- // initialize the RTREE with points on the boundary
- unsigned int MAXCOLSIZE = 0;
- BoundaryLayerColumns::iter it = _data.begin();
- for ( ; it != _data.end() ; ++it) {
- BoundaryLayerData & d = it->second;
- MAXCOLSIZE = MAXCOLSIZE > d._column.size() ? MAXCOLSIZE : d._column.size();
- MVertex * v = it->first;
- double largeMeshSize = factor*sizes[v];
- blPoint_rtree *p = new blPoint_rtree(v,largeMeshSize);
- double _min[3],_max[3];
- p->minmax (_min,_max);
- rtree.Insert(_min,_max,p);
- v2v[v] = v;
- for (unsigned int k = 0 ; k < d._column.size() ; k++)
- v2v[d._column[k]] = v;
- }
-
- // go layer by layer
- for (unsigned int LAYER = 0 ; LAYER < MAXCOLSIZE ; LAYER++){
- // store accepted points that will be inserted in the rtree
- // afterwards
- std::set<MVertex*> accepted;
- it = _data.begin();
- for ( ; it != _data.end() ; ++it) {
- MVertex * v = it->first;
- double largeMeshSize = sizes[v];
- BoundaryLayerData & d = it->second;
- // take the point if the number of layers is
- // large enough
- if (d._column.size() > LAYER){
- // check if the vertex in the column at position LAYER
- // isn't too close to another vertex
- MVertex *toCheck = d._column[LAYER];
- if (LAYER){
- double DD = toCheck->distance ( d._column[LAYER-1] );
- // do not allow to have elements that are stretched the
- // other way around !
- if (DD > largeMeshSize) largeMeshSize *= 100;
- largeMeshSize = std::max (largeMeshSize, DD);
- }
- largeMeshSize *= factor;
- bool exclude = inExclusionZone_filter (toCheck,v2v,rtree);
- if (!exclude){
- v2v [toCheck] = v;
- blPoint_rtree *p = new blPoint_rtree(toCheck,largeMeshSize);
- double _min[3],_max[3];
- p->minmax (_min,_max);
- rtree.Insert(_min,_max,p);
- }
- else {
- std::vector<MVertex*> newColumn;
- for (unsigned int k = 0 ; k < LAYER ; k++) newColumn.push_back(d._column[k]);
- for (unsigned int k = LAYER ; k < d._column.size() ; k++) delete d._column[k];
- d._column = newColumn;
- }
- }
- }
- }
-#else
- Msg::Warning ("Boundary Layer Points cannot be filtered without compiling gmsh with the rtree library");
-#endif
-}
#endif
diff --git a/Geo/boundaryLayersData.h b/Geo/boundaryLayersData.h
index ee02ed8516e7bfbdebfd1490af3e2fd075b8faa3..b626555c56dd88797c8b416fa9b6dee98f84df51 100644
--- a/Geo/boundaryLayersData.h
+++ b/Geo/boundaryLayersData.h
@@ -45,49 +45,6 @@ struct BoundaryLayerFan
BoundaryLayerFan(MEdge e1, MEdge e2 , bool s = true)
: _e1(e1),_e2(e2) , sense (s){}
};
-
-// wedge between 2 sets of continuous faces
-struct BoundaryLayerFanWedge3d
-{
- std::vector<GFace*> _gf1, _gf2;
- BoundaryLayerFanWedge3d(GFace *gf1=0, GFace *gf2=0)
- {
- if(gf1)_gf1.push_back(gf1);
- if(gf2)_gf2.push_back(gf2);
- }
- BoundaryLayerFanWedge3d(std::vector<GFace*> &gf1,
- std::vector<GFace*> &gf2):_gf1(gf1),_gf2(gf2)
- {
- }
- bool isLeft (const GFace *gf) const {
- for (unsigned int i=0;i<_gf1.size();i++)if (_gf1[i] == gf)return true;
- return false;
- }
- bool isRight (const GFace *gf) const{
- for (unsigned int i=0;i<_gf2.size();i++)if (_gf2[i] == gf)return true;
- return false;
- }
- bool isLeft (const std::vector<GFace *> &gf) const {
- for (unsigned int i=0;i<gf.size();i++)if (isLeft(gf[i]))return true;
- return false;
- }
- bool isRight (const std::vector<GFace *> &gf) const {
- for (unsigned int i=0;i<gf.size();i++)if (isRight(gf[i]))return true;
- return false;
- }
-
-};
-
-struct BoundaryLayerFanCorner3d
-{
- int _fanSize;
- std::vector<GFace *> _gf;
- BoundaryLayerFanCorner3d() : _fanSize(0){}
- BoundaryLayerFanCorner3d(int fanSize ,std::vector<GFace *> &gf)
- : _fanSize(fanSize), _gf(gf){}
-};
-
-
struct edgeColumn {
const BoundaryLayerData &_c1, &_c2;
edgeColumn(const BoundaryLayerData &c1, const BoundaryLayerData &c2)
@@ -111,9 +68,11 @@ struct faceColumn {
class BoundaryLayerColumns
{
std::map<MVertex*, BoundaryLayerFan> _fans;
- std::map<MVertex*, BoundaryLayerFanWedge3d> _wedges;
- std::map<MVertex*, BoundaryLayerFanCorner3d> _corners;
public:
+ // Element columns
+ std::map<MElement*,MElement*> _toFirst;
+ std::map<MElement*,std::vector<MElement*> > _elemColumns;
+
std::map<MFace, GFace*, Less_Face> _inverse_classification;
std::multimap<MVertex*, BoundaryLayerData> _data;
size_t size () const {return _data.size();}
@@ -145,35 +104,15 @@ public:
{
_fans.insert(std::make_pair(v,BoundaryLayerFan(e1,e2,s)));
}
- inline void addWedge(MVertex *v, GFace *gf1, GFace *gf2)
- {
- _wedges.insert(std::make_pair(v,BoundaryLayerFanWedge3d(gf1,gf2)));
- }
- inline void addWedge(MVertex *v, std::vector<GFace *>&gf1, std::vector<GFace *> &gf2)
- {
- _wedges.insert(std::make_pair(v,BoundaryLayerFanWedge3d(gf1,gf2)));
- }
- inline void addCorner(MVertex *v, int fanSize, std::vector<GFace *> &gfs)
- {
- _corners.insert(std::make_pair(v,BoundaryLayerFanCorner3d(fanSize, gfs)));
- }
- inline bool isCorner (MVertex* v) const{
- return _corners.find(v) != _corners.end();
- }
- inline bool isOnWedge (MVertex* v) const{
- return _wedges.find(v) != _wedges.end();
- }
- inline BoundaryLayerFanWedge3d getWedge (MVertex* v) {
- std::map<MVertex*, BoundaryLayerFanWedge3d>::iterator it = _wedges.find(v);
- return it->second;
- }
- inline BoundaryLayerFanCorner3d getCorner (MVertex* v) {
- std::map<MVertex*, BoundaryLayerFanCorner3d>::iterator it = _corners.find(v);
- return it->second;
+ inline const BoundaryLayerFan *getFan(MVertex *v) const{
+ std::map<MVertex*,BoundaryLayerFan>::const_iterator it = _fans.find(v);
+ if (it != _fans.end()){
+ return &it->second;
+ }
+ return 0;
}
-
const BoundaryLayerData &getColumn(MVertex *v, MFace f);
inline const BoundaryLayerData &getColumn(MVertex *v, MEdge e)
@@ -205,7 +144,6 @@ public:
static BoundaryLayerData error;
return error;
}
- void filterPoints(GEntity *ge, double factor);
};
BoundaryLayerField* getBLField(GModel *gm);
@@ -213,5 +151,4 @@ bool buildAdditionalPoints2D (GFace *gf ) ;
BoundaryLayerColumns * buildAdditionalPoints3D (GRegion *gr) ;
void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3]);
-
#endif
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index bd9cde08ae46c0cac77f2bfda9dc04568c4ade96..e13f9c2100f6ba388b73f178b01af485485a0869 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -4,6 +4,7 @@
# bugs and problems to the public mailing list <gmsh@geuz.org>.
set(SRC
+ filterElements.cpp
Generator.cpp
meshGEdge.cpp
meshGEdgeExtruded.cpp
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 4342f40e2fd2275f865970c35132bbe5f1a5056e..ea1153ff9acb6192fae49c8be62e741531658e6d 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -802,7 +802,7 @@ void GenerateMesh(GModel *m, int ask)
if(CTX::instance()->mesh.hoOptimize < 0){
ElasticAnalogy(GModel::current(), CTX::instance()->mesh.hoThresholdMin, false);
}
- else{
+ else{
OptHomParameters p;
p.nbLayers = CTX::instance()->mesh.hoNLayers;
p.BARRIER_MIN = CTX::instance()->mesh.hoThresholdMin;
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index b9f3e0255aa74e45490a4e05ba39ecd3e5c5ccec..b5f405bc9e88e5ef141340593a76a1a6313b3ebd 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -224,9 +224,30 @@ static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
bool reparamOK = true;
if(!linear){
reparamOK = reparamMeshEdgeOnFace(v0, v1, gf, p0, p1);
- if(reparamOK)
- computeEquidistantParameters(gf, p0[0], p1[0], p0[1], p1[1], nPts + 2,
- US, VS);
+ if(reparamOK) {
+ if (nPts >= 30)computeEquidistantParameters(gf, p0[0], p1[0], p0[1], p1[1], nPts + 2,
+ US, VS);
+ else {
+ US[0] = p0[0];
+ VS[0] = p0[1];
+ US[nPts+1] = p1[0];
+ VS[nPts+1] = p1[1];
+ for(int j = 0; j < nPts; j++){
+ const double t = (double)(j + 1) / (nPts + 1);
+ SPoint3 pc = edge.interpolate(t);
+ SPoint2 guess = p0 * (1.-t) + p1 * t;
+ GPoint gp = gf->closestPoint(pc, guess);
+ if(gp.succeeded()){
+ US[j+1] = gp.u();
+ VS[j+1] = gp.v();
+ }
+ else{
+ US[j+1] = guess.x();
+ VS[j+1] = guess.y();
+ }
+ }
+ }
+ }
}
std::vector<MVertex*> temp;
for(int j = 0; j < nPts; j++){
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index a7e6b12dbfe002d6a394748c6dd3e2a312e12cb5..5358debc2a9b523f767dc0a723f89585105d20c7 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -44,6 +44,7 @@
#include "multiscalePartition.h"
#include "meshGFaceLloyd.h"
#include "boundaryLayersData.h"
+#include "filterElements.h"
inline double myAngle(const SVector3 &a, const SVector3 &b, const SVector3 &d)
{
@@ -593,46 +594,11 @@ static void addOrRemove(MVertex *v1, MVertex *v2, std::set<MEdge,Less_Edge> & be
else bedges.erase(it);
}
-void filterOverlappingElements(int dim, std::vector<MElement*> &e,
- std::vector<MElement*> &eout,
- std::vector<MElement*> &einter)
-{
- eout.clear();
- MElementOctree octree (e);
- for (unsigned int i = 0; i < e.size(); ++i){
- MElement *el = e[i];
- bool intersection = false;
- for (int j=0;j<el->getNumVertices();++j){
- MVertex *v = el->getVertex(j);
- std::vector<MElement *> inters = octree.findAll(v->x(), v->y(), v->z(), dim);
- std::vector<MElement *> inters2;
- for (unsigned int k = 0; k < inters.size(); k++){
- bool found = false;
- for (int l = 0; l < inters[k]->getNumVertices(); l++){
- if (inters[k]->getVertex(l) == v)found = true;
- }
- if (!found)inters2.push_back(inters[k]);
- }
- if (inters2.size() >= 1 ){
- intersection = true;
- }
- }
- if (intersection){
- // printf("intersection found\n");
- einter.push_back(el);
- }
- else {
- eout.push_back(el);
- }
- }
-}
-
void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
{
if (!buildAdditionalPoints2D (gf))return;
BoundaryLayerColumns* _columns = gf->getColumns();
-
std::set<MEdge,Less_Edge> bedges;
std::vector<MQuadrangle*> blQuads;
@@ -650,8 +616,8 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
MVertex *v2 = (*ite)->lines[i]->getVertex(1);
MEdge dv(v1,v2);
addOrRemove(v1,v2,bedges);
-
for (unsigned int SIDE = 0 ; SIDE < _columns->_normals.count(dv); SIDE ++){
+ std::vector<MElement*> myCol;
edgeColumn ec = _columns->getColumns(v1, v2, SIDE);
const BoundaryLayerData & c1 = ec._c1;
const BoundaryLayerData & c2 = ec._c2;
@@ -672,7 +638,10 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
//avoid convergent errors
if (dv2.length() < 0.03 * dv.length())break;
- blQuads.push_back(new MQuadrangle(v11,v21,v22,v12));
+ MQuadrangle *qq = new MQuadrangle(v11,v21,v22,v12);
+ myCol.push_back(qq);
+ qq->setPartition(l);
+ blQuads.push_back(qq);
fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
v11->x(),v11->y(),v11->z(),
v12->x(),v12->y(),v12->z(),
@@ -680,7 +649,8 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
v21->x(),v21->y(),v21->z());
}
// int M = std::max(c1._column.size(),c2._column.size());
-
+ for (unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0];
+ _columns->_elemColumns[myCol[0]] = myCol;
}
}
++ite;
@@ -692,6 +662,7 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
MVertex *v = itf->first;
int nbCol = _columns->getNbColumns(v);
+ std::vector<MElement*> myCol;
for (int i=0;i<nbCol-1;i++){
const BoundaryLayerData & c1 = _columns->getColumn(v,i);
const BoundaryLayerData & c2 = _columns->getColumn(v,i+1);
@@ -709,7 +680,10 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
v12 = c2._column[l-1];
}
if (v11 != v12){
- blQuads.push_back(new MQuadrangle(v11,v12,v22,v21));
+ MQuadrangle *qq = new MQuadrangle(v11,v12,v22,v21);
+ myCol.push_back(qq);
+ qq->setPartition(l);
+ blQuads.push_back(qq);
fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
v11->x(),v11->y(),v11->z(),
v12->x(),v12->y(),v12->z(),
@@ -717,7 +691,10 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
v21->x(),v21->y(),v21->z());
}
else {
- blTris.push_back(new MTriangle(v,v22,v21));
+ MTriangle *qq = new MTriangle(v,v22,v21);
+ myCol.push_back(qq);
+ qq->setPartition(l);
+ blTris.push_back(qq);
fprintf(ff2,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
v->x(),v->y(),v->z(),
v22->x(),v22->y(),v22->z(),
@@ -725,19 +702,16 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
}
}
}
+ for (unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0];
+ _columns->_elemColumns[myCol[0]] = myCol;
}
}
fprintf(ff2,"};\n");
fclose(ff2);
- std::vector<MElement*> els,newels,oldels;
- for (unsigned int i = 0; i < blQuads.size();i++) els.push_back(blQuads[i]);
- filterOverlappingElements (2,els,newels,oldels);
- blQuads.clear();
- for (unsigned int i = 0; i < newels.size(); i++)
- blQuads.push_back((MQuadrangle*)newels[i]);
- for (unsigned int i = 0; i < oldels.size(); i++) delete oldels[i];
+
+ filterOverlappingElements (blTris,blQuads,_columns->_elemColumns,_columns->_toFirst);
for (unsigned int i = 0; i < blQuads.size();i++){
addOrRemove(blQuads[i]->getVertex(0),blQuads[i]->getVertex(1),bedges);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index e22a13eb84b1f7a13e05b6bbe48415965747cd74..f09154f43a79ca413c77a7f8272c75a7cd541f6e 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1054,7 +1054,7 @@ void bowyerWatson(GFace *gf, int MAXPNT,
}
#endif
transferDataStructure(gf, AllTris, DATA);
- removeThreeTrianglesNodes(gf);
+ // removeThreeTrianglesNodes(gf);
}
/*
@@ -1332,14 +1332,18 @@ void bowyerWatsonFrontal(GFace *gf,
*/
}
- // char name[245];
- // sprintf(name,"frontal%d-real.pos", gf->tag());
- // _printTris (name, AllTris, Us, Vs,false);
- // sprintf(name,"frontal%d-param.pos", gf->tag());
- // _printTris (name, AllTris, Us, Vs,true);
+
nbSwaps = edgeSwapPass(gf, AllTris, SWCR_QUAL, DATA);
+ /*
+ char name[245];
+ sprintf(name,"frontal%d-real.pos", gf->tag());
+ _printTris (name, AllTris.begin(), AllTris.end(), DATA,false);
+ sprintf(name,"frontal%d-param.pos", gf->tag());
+ _printTris (name, AllTris.begin(), AllTris.end(), DATA,true);
+ */
transferDataStructure(gf, AllTris, DATA);
- removeThreeTrianglesNodes(gf);
+ // removeThreeTrianglesNodes(gf);
+
// in case of boundary layer meshing
#if defined(HAVE_ANN)
{
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index bdb26f58899d15256e2fc0d8ca570b3634c95971..45e5e7073fd82452107fc95375b1294a9a34d30b 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -10,6 +10,7 @@
#include "MQuadrangle.h"
#include "STensor3.h"
#include "GEntity.h"
+#include "MFace.h"
#include <list>
#include <set>
#include <map>
@@ -127,7 +128,8 @@ class compareTri3Ptr
{
if(a->getRadius() > b->getRadius()) return true;
if(a->getRadius() < b->getRadius()) return false;
- return a<b;
+ Less_Face lf;
+ return lf(a->tri()->getFace(0), b->tri()->getFace(0));
}
};
@@ -159,11 +161,12 @@ struct edgeXface
{
v[0] = t1->tri()->getVertex(iFac == 0 ? 2 : iFac-1);
v[1] = t1->tri()->getVertex(iFac);
- if(v[0]->getNum() > v[1]->getNum()){
- MVertex *tmp = v[0];
- v[0] = v[1];
- v[1] = tmp;
- }
+ if (v[0]->getNum() > v[1]->getNum())
+ {
+ MVertex *tmp = v[0];
+ v[0] = v[1];
+ v[1] = tmp;
+ }
}
inline bool operator < ( const edgeXface &other) const
{
@@ -174,8 +177,7 @@ struct edgeXface
}
inline bool operator == ( const edgeXface &other) const
{
- if(v[0]->getNum() == other.v[0]->getNum() &&
- v[1]->getNum() == other.v[1]->getNum()) return true;
+ if(v[0]->getNum() == other.v[0]->getNum() && v[1]->getNum() == other.v[1]->getNum()) return true;
return false;
}
};
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 796873e3d7c0529f402018745bb57b942867cbd4..b61d1e3cc7ab8ca287783fb887601f7f3b1e0aa6 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -33,6 +33,7 @@
#include "Levy3D.h"
#include "directions3D.h"
#include "discreteFace.h"
+#include "filterElements.h"
#if defined(HAVE_ANN)
#include "ANN/ANN.h"
@@ -667,6 +668,108 @@ bool AssociateElementsToModelRegionWithBoundaryLayers (GRegion *gr,
}
+static int getWedge (BoundaryLayerColumns* _columns, MVertex *v1, MVertex *v2,
+ int indicesVert1 [], int indicesVert2 [])
+{
+ int N1 = _columns->getNbColumns(v1) ;
+ int N2 = _columns->getNbColumns(v2) ;
+ int fanSize = 4;
+ int NW1 = 0;
+ int NW2 = 0;
+ for (int i=0;i<N1;i++){
+ const BoundaryLayerData & c1 = _columns->getColumn(v1,i);
+ if (c1._joint.size())NW1++;
+ }
+ for (int i=0;i<N2;i++){
+ const BoundaryLayerData & c2 = _columns->getColumn(v2,i);
+ if (c2._joint.size())NW2++;
+ }
+
+
+
+
+ std::map<int,int> one2two;
+ for (int i=0;i<NW1;i++){
+ const BoundaryLayerData & c1 = _columns->getColumn(v1,i);
+ for (int j=0;j<NW2;j++){
+ const BoundaryLayerData & c2 = _columns->getColumn(v2,j);
+ for (unsigned int k=0;k<c2._joint.size();k++){
+ if (std::find(c1._joint.begin(),c1._joint.end(),c2._joint[k]) !=
+ c1._joint.end()) {
+ one2two[i] = j;
+ }
+ }
+ }
+ }
+
+ // printf("%d %d %d %d \n",N1,N2,NW1,NW2);
+ // for (std::map<int,int>::iterator it = one2two.begin(); it != one2two.end(); it++){
+ // printf("one2two[%d] = %d\n",it->first,it->second);
+ // }
+ if (one2two.size() != 2)return 0;
+
+ int vert1Start,vert1End;
+ int vert2Start,vert2End;
+ std::map<int,int>::iterator it = one2two.begin();
+ vert1Start = it->first;
+ vert2Start = it->second;
+ ++it;
+ vert1End = it->first;
+ vert2End = it->second;
+
+
+ int INDEX1, count = 0;
+ for (int i=0;i<NW1;i++){
+ for (int j=i+1;j<NW1;j++){
+ if ((vert1Start == i && vert1End == j) ||
+ (vert1Start == j && vert1End == i))
+ {
+ INDEX1 = count;
+ }
+ count++;
+ }
+ }
+ int INDEX2;
+ count = 0;
+ for (int i=0;i<NW2;i++){
+ for (int j=i+1;j<NW2;j++){
+ if ((vert2Start == i && vert2End == j) ||
+ (vert2Start == i && vert2End == j))
+ {
+ INDEX2 = count;
+ }
+ count++;
+ }
+ }
+
+ int indexVert1Start = NW1 + fanSize * ( 0 + INDEX1);
+ int indexVert1End = NW1 + fanSize * ( 1 + INDEX1);
+
+ int indexVert2Start = NW2 + fanSize * ( 0 + INDEX2);
+ int indexVert2End = NW2 + fanSize * ( 1 + INDEX2);
+
+ indicesVert1[0] = vert1Start;
+ int k=1;
+ for (int i=indexVert1Start;i< indexVert1End;++i)indicesVert1[k++] = i;
+ indicesVert1[fanSize+1] = vert1End;
+
+ indicesVert2[0] = vert2Start;
+ k = 1;
+ if (indexVert2End > indexVert2Start){
+ for (int i=indexVert2Start;i< indexVert2End;++i)indicesVert2[k++] = i;
+ }
+ else {
+ for (int i=indexVert2Start-1;i<= indexVert2End;--i)indicesVert2[k++] = i;
+ }
+ indicesVert2[fanSize+1] = vert2End;
+
+
+ // printf("%d %d %d %d %d %d %d %d\n",vert1Start,vert1End,vert2Start,vert2End,indexVert1Start,indexVert1End,indexVert2Start,indexVert2End);
+ // return 0;
+
+ return fanSize + 2;
+}
+
static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
{
if (getBLField(gr->model())) insertVerticesInRegion(gr,-1);
@@ -681,8 +784,6 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
std::list<GFace*> embedded_faces = gr->embeddedFaces();
faces.insert(faces.begin(), embedded_faces.begin(),embedded_faces.end());
- FILE *ff2 = fopen ("tato3D.pos","w");
- fprintf(ff2,"View \" \"{\n");
std::set<MVertex*> verts;
std::list<GFace*>::iterator itf = faces.begin();
@@ -700,6 +801,7 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
const BoundaryLayerData & c3 = fc._c3;
int N = std::min(c1._column.size(),std::min(c2._column.size(),c3._column.size()));
// printf("%d Layers\n",N);
+ std::vector<MElement*> myCol;
for (int l=0;l < N ;++l){
MVertex *v11,*v12,*v13,*v21,*v22,*v23;
v21 = c1._column[l];
@@ -716,15 +818,15 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
v13 = c3._column[l-1];
}
// printf("coucoucouc %p %p %p %p %p %p\n",v11,v12,v13,v21,v22,v23);
- blPrisms.push_back(new MPrism(v11,v12,v13,v21,v22,v23));
- fprintf(ff2,"SI (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1,1,1};\n",
- v11->x(),v11->y(),v11->z(),
- v12->x(),v12->y(),v12->z(),
- v13->x(),v13->y(),v13->z(),
- v21->x(),v21->y(),v21->z(),
- v22->x(),v22->y(),v22->z(),
- v23->x(),v23->y(),v23->z());
- // printf("done %d\n",l);
+ MPrism *prism = new MPrism(v11,v12,v13,v21,v22,v23);
+ // store the layer the element belongs
+ prism->setPartition(l+1);
+ blPrisms.push_back(prism);
+ myCol.push_back(prism);
+ }
+ if (!myCol.empty()){
+ for (unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0];
+ _columns->_elemColumns[myCol[0]] = myCol;
}
}
}
@@ -748,167 +850,77 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
MEdge e = *ite;
MVertex *v1 = e.getVertex(0);
MVertex *v2 = e.getVertex(1);
- bool onWedge1 = _columns->isOnWedge(v1);
- bool onWedge2 = _columns->isOnWedge(v2);
- bool onCorner1 = _columns->isCorner(v1);
- bool onCorner2 = _columns->isCorner(v2);
- if ((onWedge1 || onCorner1) && (onWedge2 || onCorner2)){
- int N1=0,N2=0;
- std::vector<GFace *>gfs1,gfs2;
- BoundaryLayerFanWedge3d w1,w2;
- BoundaryLayerFanCorner3d c1,c2;
- if (onWedge1) {
- N1 = _columns->getNbColumns(v1);
- w1 = _columns->getWedge(v1);
- gfs1 = w1._gf1;
- gfs2 = w1._gf2;
- }
- else {
- c1 = _columns->getCorner(v1);
- N1 = c1._fanSize;
- }
- if (onWedge2) {
- N2 = _columns->getNbColumns(v2);
- w2 = _columns->getWedge(v2);
- gfs1 = w2._gf1;
- gfs2 = w2._gf2;
- }
- else {
- c2 = _columns->getCorner(v2);
- N2 = c2._fanSize;
- }
-
- // have to go from gf1 to gf2 to be aware of the sense!!!
- if (N1 == N2){
- for (int i = 0; i < N1 - 1; i++){
-
- unsigned int i11=i, i12=i+1, i21=i,i22=i+1;
-
- if (onWedge1){
- // if (w1.isLeft(gfs1)){i11=i;i12=i+1;}
- // else if (w1.isRight(gfs1)){i11=N1-1-i;i12=N1-2-i;}
- // printf("%d %d %d\n",w1.isLeft(gfs1),gfs1[0] == w1._gf1[0],gfs1.size());
- // printf("%d %d %d\n",w1.isRight(gfs1),gfs1[0] == w1._gf2[0],gfs1.size());
- if (gfs1[0] == w1._gf1[0]){i11=i;i12=i+1;}
- else if (gfs1[0] == w1._gf2[0]){i11=N1-1-i;i12=N1-2-i;}
- }
- else {
- /*
-
-
-
- | 0 --> column 0
- |
- | fan with N1 columns
- gf.size |
- ---------- + ---------- 1 --> column N1-1
- |
- |
- |
- | 2 --> column 2*(N1-1)
-
- 0 1 2 3 --> 0 -> 4-1
- 3 4 5 6 --> 4-1 -> 2 (4-1)
- 6 7 8 0
-
- 1 --> 0 ==>
-
- */
- int K1 = -1, K2 = -1;
- for (unsigned int s=0;s < c1._gf.size();s++){
- if (w2.isLeft(c1._gf[s]))K1 = s;
- if (w2.isRight(c1._gf[s]))K2 = s;
- }
- if (K1==-1) Msg::Error("K1 = -1");
- if (K2==-1) Msg::Error("K2 = -1");
- if (K1+1==K2){ i11=K1*(N1-1)+i;i12=i11+1; }
- else if (K2+1 == K1){ i11=K1*(N1-1)-i;i12=i11-1; }
- else if (K2 == 0 && K1 == (int)c1._gf.size() - 1){ i11=K1*(N1-1)+i;i12=i11+1; }
- else if (K1 == 0 && K2 == (int)c1._gf.size() - 1){ i11=(K2+1)*(N1-1)-i;i12=i11-1; }
- else Msg::Error("KROUPOUK 1 %d %d !", K1, K2);
- if (i12 == (c1._gf.size()) * (N1-1))i12=0;
- if (i11 == (c1._gf.size()) * (N1-1))i11=0;
- }
- if (onWedge2){
- if (w2.isLeft(gfs1)){ i21=i;i22=i+1; }
- else if (w2.isRight(gfs1)){ i21=N1-1-i;i22=N1-2-i; }
- }
- else {
- int K1 = -1, K2 = -1;
- for (unsigned int s=0;s<c2._gf.size();s++){
- if (w1.isLeft(c2._gf[s]))K1 = s;
- if (w1.isRight(c2._gf[s]))K2 = s;
- }
- if (K1+1 == K2){ i21=K1*(N1-1)+i;i22=i21+1; }
- else if (K2+1 == K1){ i21=K1*(N1-1)-i;i22=i21-1; }
- else if (K2 == 0 && K1 == (int)c2._gf.size()-1 ){ i21=K1*(N1-1)+i;i22=i21+1; }
- else if (K1 == 0 && K2 == (int)c2._gf.size()-1){ i21=(K2+1)*(N1-1)-i;i22=i21-1; }
- else Msg::Error("KROUPOUK 2 %d %d !",K1,K2);
- if (i22 == (c2._gf.size()) * (N1-1))i22=0;
- if (i21 == (c2._gf.size()) * (N1-1))i21=0;
- }
-
- BoundaryLayerData c11 = _columns->getColumn(v1,i11);
- BoundaryLayerData c12 = _columns->getColumn(v1,i12);
- BoundaryLayerData c21 = _columns->getColumn(v2,i21);
- BoundaryLayerData c22 = _columns->getColumn(v2,i22);
- int N = std::min(c11._column.size(),
- std::min(c12._column.size(),
- std::min(c21._column.size(), c22._column.size())));
- for (int l=0;l < N ;++l){
- MVertex *v11,*v12,*v13,*v14;
- MVertex *v21,*v22,*v23,*v24;
- v21 = c11._column[l];
- v22 = c12._column[l];
- v23 = c22._column[l];
- v24 = c21._column[l];
- if (l == 0){
- v11 = v12 = v1;
- v13 = v14 = v2;
- }
- else {
- v11 = c11._column[l-1];
- v12 = c12._column[l-1];
- v13 = c22._column[l-1];
- v14 = c21._column[l-1];
- }
-
- if (l == 0){
- blPrisms.push_back(new MPrism(v12,v21,v22,v13,v24,v23));
- fprintf(ff2,"SI (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){3,3,3,3,3,3};\n",
- v12->x(),v12->y(),v12->z(),
- v21->x(),v21->y(),v21->z(),
- v22->x(),v22->y(),v22->z(),
- v13->x(),v13->y(),v13->z(),
- v24->x(),v24->y(),v24->z(),
- v23->x(),v23->y(),v23->z());
- }
- else {
- blHexes.push_back(new MHexahedron(v11,v12,v13,v14,v21,v22,v23,v24));
- fprintf(ff2,"SH (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2,2,2,2,2,2};\n",
- v11->x(),v11->y(),v11->z(),
- v12->x(),v12->y(),v12->z(),
- v13->x(),v13->y(),v13->z(),
- v14->x(),v14->y(),v14->z(),
- v21->x(),v21->y(),v21->z(),
- v22->x(),v22->y(),v22->z(),
- v23->x(),v23->y(),v23->z(),
- v24->x(),v24->y(),v24->z());
- }
- }
+ int indices1[256];
+ int indices2[256];
+ int NbW = getWedge (_columns, v1, v2, indices1,indices2);
+ for (int i=0;i<NbW-1;i++){
+ int i11 = indices1[i];
+ int i12 = indices1[i+1];
+ int i21 = indices2[i];
+ int i22 = indices2[i+1];
+ BoundaryLayerData c11 = _columns->getColumn(v1,i11);
+ BoundaryLayerData c12 = _columns->getColumn(v1,i12);
+ BoundaryLayerData c21 = _columns->getColumn(v2,i21);
+ BoundaryLayerData c22 = _columns->getColumn(v2,i22);
+ int N = std::min(c11._column.size(),
+ std::min(c12._column.size(),
+ std::min(c21._column.size(), c22._column.size())));
+ std::vector<MElement*> myCol;
+ for (int l=0;l < N ;++l){
+ MVertex *v11,*v12,*v13,*v14;
+ MVertex *v21,*v22,*v23,*v24;
+ v21 = c11._column[l];
+ v22 = c12._column[l];
+ v23 = c22._column[l];
+ v24 = c21._column[l];
+ if (l == 0){
+ v11 = v12 = v1;
+ v13 = v14 = v2;
+ }
+ else {
+ v11 = c11._column[l-1];
+ v12 = c12._column[l-1];
+ v13 = c22._column[l-1];
+ v14 = c21._column[l-1];
+ }
+
+ if (l == 0){
+ MPrism *prism = new MPrism(v12,v21,v22,v13,v24,v23);
+ // store the layer the element belongs
+ prism->setPartition(l+1);
+ myCol.push_back(prism);
+
+ blPrisms.push_back(prism);
+ }
+ else {
+ MHexahedron *hex = new MHexahedron(v11,v12,v13,v14,v21,v22,v23,v24);
+ // store the layer the element belongs
+ myCol.push_back(hex);
+ hex->setPartition(l+1);
+ blHexes.push_back(hex);
}
}
- else{
- Msg::Error("cannot create 3D BL FAN %d %d -- %d %d %d %d",
- N1,N2,onWedge1,onWedge1,onCorner1,onCorner2);
+ if (!myCol.empty()){
+ for (unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0];
+ _columns->_elemColumns[myCol[0]] = myCol;
}
}
++ite;
}
- fprintf(ff2,"};\n");
- fclose(ff2);
-
+ filterOverlappingElements (blPrisms,blHexes,_columns->_elemColumns,_columns->_toFirst);
+ {
+ FILE *ff2 = fopen ("tato3D.pos","w");
+ fprintf(ff2,"View \" \"{\n");
+ for (unsigned int i = 0; i < blPrisms.size();i++){
+ blPrisms[i]->writePOS(ff2,1,0,0,0,0,0);
+ }
+ for (unsigned int i = 0; i < blHexes.size();i++){
+ blHexes[i]->writePOS(ff2,1,0,0,0,0,0);
+ }
+ fprintf(ff2,"};\n");
+ fclose(ff2);
+ }
for (unsigned int i = 0; i < blPrisms.size();i++){
for (unsigned int j=0;j<5;j++)
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index e4f78f93a1e821ce0e1759873d97a2fd7e39abf0..6e726cf0997a1b3f5d4627f56e4f117a971ba7d7 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -27,6 +27,7 @@ set(SRC
GaussQuadraturePyr.cpp
GaussLegendreSimplex.cpp
GaussJacobi1D.cpp
+ geodesics.cpp
robustPredicates.cpp
mathEvaluator.cpp
Iso.cpp
diff --git a/benchmarks/2d/Square-02.geo b/benchmarks/2d/Square-02.geo
index 70c16d58ad5b2437d44faef176bdefb25e685eaa..b2612898bd2e83282c9fa358ca3ca28626c1817d 100644
--- a/benchmarks/2d/Square-02.geo
+++ b/benchmarks/2d/Square-02.geo
@@ -13,4 +13,4 @@ Line(3) = {1,4};
Line(4) = {4,3};
Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};
-Recombine Surface {6};
+//Recombine Surface {6};
diff --git a/benchmarks/2d/square.geo b/benchmarks/2d/square.geo
index 6c186bf2f90bf7915d17de352e0192800338b7fb..f3050fb2548065f88789884592b0ed4ba884daa7 100644
--- a/benchmarks/2d/square.geo
+++ b/benchmarks/2d/square.geo
@@ -1,4 +1,4 @@
-lc=1;
+lc=1000;
Point(1) = {0, 0, 0,lc*.1};
Point(2) = {0, 10, 0,lc};
Point(3) = {10, 10, 0,lc};
@@ -11,6 +11,9 @@ Line(4) = {1, 2};
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(10) = {5};
+Physical Line("wall")={1,2,3,4};
+Physical Surface("air")={10};
+
//----------------------
//Compound Line(10)={1,2,3,4};
@@ -24,4 +27,4 @@ Plane Surface(10) = {5};
-Recombine Surface {10};
+//Recombine Surface {10};