From bebbb3d8b6fefdcc8167d299c5045d099051c376 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 22 Nov 2013 10:22:03 +0000
Subject: [PATCH] high order points better positionned initially boundary
layers in 3D are advancing (filtering intersecting layers using rtrees)
randomness of 2D meshes is now close to be ok (we still have problems with
OCC ...)
---
Geo/CMakeLists.txt | 3 +-
Geo/GFace.cpp | 1 -
Geo/GModel.cpp | 9 +
Geo/GModel.h | 1 +
Geo/GModelFactory.cpp | 17 ++
Geo/GModelIO_OCC.cpp | 9 +-
Geo/MElementOctree.cpp | 4 +-
Geo/MVertex.cpp | 11 -
Geo/OCCEdge.cpp | 26 +-
Geo/boundaryLayersData.cpp | 411 ++++++++++------------------
Geo/boundaryLayersData.h | 83 +-----
Mesh/CMakeLists.txt | 1 +
Mesh/Generator.cpp | 2 +-
Mesh/HighOrder.cpp | 27 +-
Mesh/meshGFace.cpp | 68 ++---
Mesh/meshGFaceDelaunayInsertion.cpp | 18 +-
Mesh/meshGFaceDelaunayInsertion.h | 18 +-
Mesh/meshGRegion.cpp | 342 ++++++++++++-----------
Numeric/CMakeLists.txt | 1 +
benchmarks/2d/Square-02.geo | 2 +-
benchmarks/2d/square.geo | 7 +-
21 files changed, 455 insertions(+), 606 deletions(-)
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 4f8685ff42..514533cc4e 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 9f1620c42c..238c63b82b 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 096918d280..d2e04c5aeb 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 f026f38b6c..8fc5068a07 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 19b29e2806..4a7cbeae41 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 42bb88ce9e..f6b8a077f8 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 ce8db23964..e65860ca16 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 40fff3e7cf..98dceccb73 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 589d9a1668..23c59eba93 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 7d572a276c..c6430d9261 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 ee02ed8516..b626555c56 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 bd9cde08ae..e13f9c2100 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 4342f40e2f..ea1153ff9a 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 b9f3e0255a..b5f405bc9e 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 a7e6b12dbf..5358debc2a 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 e22a13eb84..f09154f43a 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 bdb26f5889..45e5e7073f 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 796873e3d7..b61d1e3cc7 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 e4f78f93a1..6e726cf099 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 70c16d58ad..b2612898bd 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 6c186bf2f9..f3050fb254 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};
--
GitLab