diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index abf1d2764e87e33fd913620cdd30651536ac4923..e9e5bdfff12f4f21ec267565cdc98073be216d7d 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -2001,12 +2001,13 @@ double GFaceCompound::locCurvature(MTriangle *t, double u, double v) const
SPoint2 GFaceCompound::parFromVertex(MVertex *v) const
{
double U=0.0,V=0.0;
- if(v->onWhat()->dim()==2){// && meshStatistics.status == GFace::DONE) {
+ if(v->onWhat()->dim()==2){
v->getParameter(0, U);
v->getParameter(1, V);
}
if (v->onWhat()->dim()==1 ||
- (v->onWhat()->dim()==2 && U == -1 && V==-1)){ //if MFaceVertex created on edge in bunin
+ (v->onWhat()->dim()==2 &&
+ U == -1 && V==-1)){ //if MFaceVertex created on edge in bunin
SPoint2 sp = getCoordinates(v);
U = sp.x();
V = sp.y();
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index e5a1e4c254ee4edcdeef76caf51b4d2692e7243a..3a166fdc20216fdc95d39e3d36ecc53a23f2b4dc 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1268,6 +1268,19 @@ static void _associateEntityWithElementVertices(GEntity *ge, std::vector<T*> &el
}
}
+void GModel:: _createGeometryOfDiscreteEntities() {
+ for(eiter it = firstEdge(); it != lastEdge(); ++it){
+ if((*it)->geomType() == GEntity::DiscreteCurve) {
+ discreteEdge *de = dynamic_cast<discreteEdge*> (*it);
+ if (!de)Msg::Fatal("no fun");
+ de->createGeometry();
+ }
+ }
+ for(fiter it = firstFace(); it != lastFace(); ++it){
+ }
+}
+
+
void GModel::_associateEntityWithMeshVertices()
{
// loop on regions, then on faces, edges and vertices and store the
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 7fc00f5e61ba5f7667d6dcc4f6b73869d6d68160..27faded9600f6072f6c42d37d2a6c68947a3de63 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -119,6 +119,9 @@ class GModel
// the fly if needed
void _storeElementsInEntities(std::map<int, std::vector<MElement*> > &map);
+ // Discrete Entities have to have their mesh moved to a geometry container
+ void _createGeometryOfDiscreteEntities();
+
// loop over all vertices connected to elements and associate
// geometrical entity
void _associateEntityWithMeshVertices();
diff --git a/Geo/GModelIO_MSH.cpp b/Geo/GModelIO_MSH.cpp
index 34656bf3b09e3b0dd8589fe862c5788f231aa852..1115b54deaef1d71d0cefe1c9a32e8852dc4c0c1 100644
--- a/Geo/GModelIO_MSH.cpp
+++ b/Geo/GModelIO_MSH.cpp
@@ -449,6 +449,7 @@ int GModel::readMSH(const std::string &name)
else
_storeVerticesInEntities(_vertexMapCache);
+ _createGeometryOfDiscreteEntities() ;
fclose(fp);
return postpro ? 2 : 1;
diff --git a/Geo/GModelIO_MSH2.cpp b/Geo/GModelIO_MSH2.cpp
index 2ea895cd8f2587f98d061572087c54b59fe9ec19..a450488c72e944f68f717b9eb028b66f78c77c56 100644
--- a/Geo/GModelIO_MSH2.cpp
+++ b/Geo/GModelIO_MSH2.cpp
@@ -689,6 +689,7 @@ int GModel::_readMSH2(const std::string &name)
for(int i = 0; i < 4; i++)
_storePhysicalTagsInEntities(i, physicals[i]);
+ _createGeometryOfDiscreteEntities() ;
fclose(fp);
return postpro ? 2 : 1;
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 793c30ece00b2f9431f1d97da5ff2c67e542043a..5ccf30023d59d9ba6e58672a9a996f11478f19aa 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -271,6 +271,7 @@ void discreteEdge::parametrize(std::map<GFace*, std::map<MVertex*, MVertex*,
std::map<GRegion*, std::map<MVertex*, MVertex*,
std::less<MVertex*> > > ®ion2Vert)
{
+ return;
if (_pars.empty()){
for (unsigned int i = 0; i < lines.size() + 1; i++){
_pars.push_back(i);
@@ -391,8 +392,7 @@ bool discreteEdge::getLocalParameter(const double &t, int &iLine,
double tmin = _pars[iLine];
double tmax = _pars[iLine+1];
if (t >= tmin && t <= tmax){
- tLoc = _orientation[iLine] ? (t-tmin)/(tmax-tmin) :
- 1 - (t-tmin)/(tmax-tmin);
+ tLoc = (t-tmin)/(tmax-tmin);
return true;
}
}
@@ -406,8 +406,8 @@ GPoint discreteEdge::point(double par) const
if(!getLocalParameter(par, iEdge, tLoc)) return GPoint();
double x, y, z;
- MVertex *vB = lines[iEdge]->getVertex(0);
- MVertex *vE = lines[iEdge]->getVertex(1);
+ MVertex *vB = discrete_lines[iEdge]->getVertex(0);
+ MVertex *vE = discrete_lines[iEdge]->getVertex(1);
//linear Lagrange mesh
x = vB->x() + tLoc * (vE->x()- vB->x());
@@ -422,14 +422,14 @@ SVector3 discreteEdge::firstDer(double par) const
int iEdge;
if(!getLocalParameter(par, iEdge, tLoc)) return SVector3();
- MVertex *vB = lines[iEdge]->getVertex(0);
- MVertex *vE = lines[iEdge]->getVertex(1);
+ MVertex *vB = discrete_lines[iEdge]->getVertex(0);
+ MVertex *vE = discrete_lines[iEdge]->getVertex(1);
double dx, dy, dz;
- double dt = 1.0;
- dx = (vE->x() - vB->x()) / dt;
- dy = (vE->y() - vB->y()) / dt;
- dz = (vE->z() - vB->z()) / dt;
+ // double dt = 1.0;
+ dx = (vE->x() - vB->x());// / dt;
+ dy = (vE->y() - vB->y());// / dt;
+ dz = (vE->z() - vB->z());// / dt;
SVector3 der(dx, dy, dz);
return der;
@@ -470,36 +470,43 @@ double discreteEdge::curvatures(const double par, SVector3 *dirMax, SVector3 *di
Range<double> discreteEdge::parBounds(int i) const
{
- return Range<double>(0, 1);
+ return Range<double>(0, (double)discrete_lines.size());
}
void discreteEdge::createGeometry(){
if (discrete_lines.empty()){
+ Msg::Info("creating the geometry of discrete curve %d",tag());
createTopo();
// copy the mesh
std::map<MVertex*,MVertex*> v2v;
for (unsigned int i = 0; i < mesh_vertices.size(); i++){
- MVertex *v = new MVertex(mesh_vertices[i]->x(),mesh_vertices[i]->y(),mesh_vertices[i]->z());
- v2v[mesh_vertices[i]] = v;
- discrete_vertices.push_back(v);
+ MVertex *v = new MEdgeVertex(mesh_vertices[i]->x(),mesh_vertices[i]->y(),mesh_vertices[i]->z(),this,(double)(i+1));
+ v2v [mesh_vertices[i]] = v;
}
+
+ std::vector<MLine*> _temp;
for (unsigned int i = 0; i < lines.size(); i++){
MVertex *v0 = lines[i]->getVertex(0);
MVertex *v1 = lines[i]->getVertex(1);
- discrete_lines.push_back(new MLine(v0,v1));
- }
- // compute parameters
- double L = 0.0;
- _pars.push_back(L);
- for (unsigned int i = 0; i < discrete_lines.size(); i++){
- MVertex *v0 = discrete_lines[i]->getVertex(0);
- MVertex *v1 = discrete_lines[i]->getVertex(1);
- L += distance(v0,v1);
- _pars.push_back(L);
+ MVertex *v00 = (v2v.find(v0) == v2v.end()) ? v0 : v2v[v0];
+ MVertex *v01 = (v2v.find(v1) == v2v.end()) ? v1 : v2v[v1];
+ if (_orientation[i] == 1){
+ discrete_lines.push_back(new MLine(v00,v01));
+ }
+ else{
+ discrete_lines.push_back(new MLine(v01,v00));
+ }
}
- for (unsigned int i = 0; i < _pars.size(); i++){
- _pars[i] /= L;
+
+
+ // compute parameters and recompute the vertices
+ _pars.push_back(0.0);
+ for (unsigned int i = 1; i < discrete_lines.size(); i++){
+ _pars.push_back((double)i);
+ MVertex *newv = discrete_lines[i]->getVertex(0);
+ discrete_vertices.push_back(newv);
}
+ _pars.push_back((double)discrete_lines.size());
}
}
@@ -508,8 +515,8 @@ void discreteEdge::mesh(bool verbose){
#if defined(HAVE_MESH)
// copy the mesh into geometrical entities
// FIXME
- return;
- createGeometry();
+ // return;
+ // createGeometry();
meshGEdge mesher;
mesher(this);
#endif
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index 6a09a39555eb8de914f19a31531d8cc526eb7a9b..7ece5ad825eda60dbd0a3c9877642798e8c7138b 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -42,6 +42,7 @@ class discreteEdge : public GEdge {
void computeNormals () const;
virtual void mesh(bool) ;
void writeGEO(FILE *fp);
+ int minimumDrawSegments() const {return 2*_pars.size();}
};
#endif
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 46aa13a4105d573b8f9e0e3520d658662d941bbd..3cd561be92ae489a6cdfc204d9edfef9e27670ef 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -139,7 +139,7 @@ void discreteFace::secondDer(const SPoint2 ¶m,
}
// FIXME PAB ----> really create an atlas !!!!!!!!!!!!!!!
-void discreteFace::createAtlas() {
+void discreteFace::createGeometry() {
if (!_atlas.empty())return;
// parametrization is done here !!!
discreteDiskFace *df = new discreteDiskFace (this, triangles);
@@ -160,7 +160,8 @@ void discreteFace::gatherMeshes() {
void discreteFace::mesh(bool verbose) {
#if defined(HAVE_MESH)
- createAtlas();
+ return;
+ createGeometry();
for (unsigned int i=0;i<_atlas.size();i++)
_atlas[i]->mesh(verbose);
gatherMeshes();
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index f8a91f1a91e8517547354969539568f171be301a..0fd0ac93f42e91e36629542f5d7e44fca2f2a23b 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -34,7 +34,7 @@ class discreteFace : public GFace {
void setBoundEdges(GModel *gm, std::vector<int> tagEdges);
void findEdges(std::map<MEdge, std::vector<int>, Less_Edge > &map_edges);
void writeGEO(FILE *fp);
- void createAtlas();
+ void createGeometry();
void gatherMeshes();
virtual void mesh (bool verbose);
std::vector<discreteDiskFace*> _atlas;
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index 3e3263e7bc99d8ab3c32e15ac848de7e7265ceb1..2071835ab26a31a2582672e39ca15b50e96ad092 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -118,8 +118,13 @@ std::string gmshEdge::getAdditionalInfoString()
int gmshEdge::minimumMeshSegments () const
{
int np;
- if(geomType() == Line)
+ if(geomType() == Line){
np = GEdge::minimumMeshSegments();
+ // FIXME FOR QUADS
+ if(List_Nbr(c->Control_Points) > 2){
+ np = 3*(List_Nbr(c->Control_Points))+1;
+ }
+ }
else if(geomType() == Circle || geomType() == Ellipse)
np = (int)(fabs(c->Circle.t1 - c->Circle.t2) *
(double)CTX::instance()->mesh.minCircPoints / M_PI) - 1;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 3f3f359396df0bf5b867733f62726677f0f7da07..b0b0e363e3427a86a5accfe06e2b31546a2814bc 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -53,7 +53,6 @@ void backgroundMesh::unset()
_current = 0;
}
-double backgroundMesh::sizeFactor = 1.0;
backgroundMesh::backgroundMesh(GFace *_gf, bool cfd)
#if defined(HAVE_ANN)
: _octree(0), uv_kdtree(0), nodes(0), angle_nodes(0), angle_kdtree(0)
@@ -503,21 +502,21 @@ void backgroundMesh::updateSizes(GFace *_gf)
MVertex *v = _2Dto3D[itv->first];
double lc;
if (v->onWhat()->dim() == 0){
- lc = sizeFactor * BGM_MeshSize(v->onWhat(), 0,0,v->x(),v->y(),v->z());
+ lc = BGM_MeshSize(v->onWhat(), 0,0,v->x(),v->y(),v->z());
}
else if (v->onWhat()->dim() == 1){
double u;
v->getParameter(0, u);
- lc = sizeFactor * BGM_MeshSize(v->onWhat(), u, 0, v->x(), v->y(), v->z());
+ lc = BGM_MeshSize(v->onWhat(), u, 0, v->x(), v->y(), v->z());
}
else{
reparamMeshVertexOnFace(v, _gf, p);
- lc = sizeFactor * BGM_MeshSize(_gf, p.x(), p.y(), v->x(), v->y(), v->z());
+ lc = BGM_MeshSize(_gf, p.x(), p.y(), v->x(), v->y(), v->z());
}
// printf("2D -- %g %g 3D -- %g %g\n",p.x(),p.y(),v->x(),v->y());
itv->second = std::min(lc,itv->second);
- itv->second = std::max(itv->second, sizeFactor * CTX::instance()->mesh.lcMin);
- itv->second = std::min(itv->second, sizeFactor * CTX::instance()->mesh.lcMax);
+ itv->second = std::max(itv->second, CTX::instance()->mesh.lcMin);
+ itv->second = std::min(itv->second, CTX::instance()->mesh.lcMax);
}
// do not allow large variations in the size field
// (Int. J. Numer. Meth. Engng. 43, 1143-1165 (1998) MESH GRADATION
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index 7ad287d6a8a6ce14aeb860bd53001aea840b19bf..0fcd71011ab2bd4ff365b86b25d35d6c22c99922 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -40,7 +40,6 @@ struct crossField2d
class backgroundMesh : public simpleFunction<double>
{
- static double sizeFactor;
MElementOctree *_octree;
std::vector<MVertex*> _vertices;
std::vector<MElement*> _triangles;
@@ -66,7 +65,6 @@ class backgroundMesh : public simpleFunction<double>
static void setCrossFieldsByDistance(GFace *);
static void unset();
static backgroundMesh *current () { return _current; }
- static void setSizeFactor (double s) {sizeFactor = s;}
void propagate1dMesh(GFace *);
void propagateCrossField(GFace *, simpleFunction<double> *);
void propagateCrossFieldHJ(GFace *);
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 652b7d351faa7df64a395ab261e4352df00b202b..8f04b245277c5fc6ce2583268aeb01be450eb3f2 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -5,6 +5,7 @@
#include "GmshConfig.h"
#include "GModel.h"
+#include "discreteEdge.h"
#include "meshGEdge.h"
#include "GEdge.h"
#include "MLine.h"
@@ -74,6 +75,13 @@ static double smoothPrimitive(GEdge *ge, double alpha,
return Points[Points.size() - 1].p;
}
+
+static double F_LcB(GEdge *ge, double t)
+{
+ GPoint p = ge->point(t);
+ return BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
+}
+
static double F_Lc(GEdge *ge, double t)
{
GPoint p = ge->point(t);
@@ -297,7 +305,11 @@ void copyMesh(GEdge *from, GEdge *to, int direction)
void deMeshGEdge::operator() (GEdge *ge)
{
- if(ge->geomType() == GEntity::DiscreteCurve) return;
+ if(ge->geomType() == GEntity::DiscreteCurve) {
+ // FIXME : NOTHING SPECIAL TO DO
+ return;
+ }
+
ge->deleteMesh();
ge->meshStatistics.status = GEdge::PENDING;
ge->correspondingVertices.clear();
@@ -327,6 +339,65 @@ static int increaseN (int N)
return N;
}
+// ensure not to have points that are too close to each other.
+// can be caused by a coarse 1D mesh or by a noisy curve
+static void filterPoints (GEdge*ge) {
+ if (ge->mesh_vertices.empty())return;
+ if(ge->meshAttributes.method == MESH_TRANSFINITE)return;
+ //if (ge->mesh_vertices.size() <=3)return;
+ bool forceOdd = false;
+ if((ge->meshAttributes.method != MESH_TRANSFINITE ||
+ CTX::instance()->mesh.flexibleTransfinite) &&
+ CTX::instance()->mesh.algoRecombine != 0){
+ if(CTX::instance()->mesh.recombineAll){
+ forceOdd = true;
+ }
+ }
+
+ MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
+ std::vector<std::pair<double, MVertex*> > lengths;
+ for (unsigned int i=0;i<ge->mesh_vertices.size();i++){
+ MEdgeVertex *v = dynamic_cast<MEdgeVertex*> (ge->mesh_vertices[i]);
+ if (!v)Msg::Fatal("in 1D mesh");
+ double d = distance (v,v0);
+ double t;
+ v->getParameter(0,t);
+ if (i != 0){
+ double t0;
+ v0->getParameter(0,t0);
+ t=0.5*(t+t0);
+ }
+ double lc = F_LcB(ge, t);
+ // double lc = v->getLc();
+ if (d < lc * .3) {
+ lengths.push_back(std::make_pair(lc/d,v));
+ }
+ else
+ v0=v;
+ }
+ std::sort(lengths.begin(),lengths.end());
+ int last = lengths.size();
+ if (forceOdd) {
+ while (last %2 != 0)last--;
+ }
+ /*
+ if (CTX::instance()->mesh.algoRecombine == 2){
+ if (last < 4)last = 0;
+ while (last %4 != 0)last--;
+ }
+ else {
+ while (last %2 != 0)last--;
+ }
+ }
+ */
+
+ for (int i=0;i<last;i++){
+ std::vector<MVertex*>::iterator it = std::find(ge->mesh_vertices.begin(),ge->mesh_vertices.end(),lengths[i].second);
+ ge->mesh_vertices.erase(it);
+ delete lengths[i].second;
+ }
+}
+
void meshGEdge::operator() (GEdge *ge)
{
#if defined(HAVE_ANN)
@@ -499,6 +570,10 @@ void meshGEdge::operator() (GEdge *ge)
mesh_vertices.resize(NUMP - 1);
}
+ // printf("%ld ----> ", ge->mesh_vertices.size());
+ filterPoints (ge);
+ // printf("%ld \n", ge->mesh_vertices.size());
+
for(unsigned int i = 0; i < mesh_vertices.size() + 1; i++){
MVertex *v0 = (i == 0) ?
ge->getBeginVertex()->mesh_vertices[0] : mesh_vertices[i - 1];
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index d42f09ff280d32cbb75e42c8c722b50548bea67a..38d856778bf13558654f2928009ba7cbc7847765 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -98,6 +98,9 @@ public:
MVertex *v1 = (*ite)->lines[i]->getVertex(0);
MVertex *v2 = (*ite)->lines[i]->getVertex(1);
MVertex *v3 = (*ite)->lines[i+1]->getVertex(1);
+ v2->x() = 0.5*(v1->x()+v3->x());
+ v2->y() = 0.5*(v1->y()+v3->y());
+ v2->z() = 0.5*(v1->z()+v3->z());
temp.push_back(new MLine(v1,v3));
if (v1->onWhat() == *ite) (*ite)->mesh_vertices.push_back(v1);
_middle[MEdge(v1,v3)] = v2;
@@ -108,8 +111,8 @@ public:
}
++ite;
}
+ CTX::instance()->mesh.lcFactor *=2.0;
}
- backgroundMesh::setSizeFactor(2.0);
}
void subdivide ()
{
@@ -129,7 +132,10 @@ public:
m[j] = p1;
if (it == _middle.end() && it2 == eds.end()){
GPoint gp = _gf->point((p1+p2)*0.5);
- v[j] = new MFaceVertex (gp.x(),gp.y(),gp.z(),_gf,gp.u(),gp.v());
+ double XX = 0.5*(E.getVertex(0)->x()+E.getVertex(1)->x());
+ double YY = 0.5*(E.getVertex(0)->y()+E.getVertex(1)->y());
+ double ZZ = 0.5*(E.getVertex(0)->z()+E.getVertex(1)->z());
+ v[j] = new MFaceVertex (XX,YY,ZZ,_gf,gp.u(),gp.v());
_gf->mesh_vertices.push_back(v[j]);
eds[E] = v[j];
}
@@ -142,7 +148,10 @@ public:
}
}
GPoint gp = _gf->point((m[0]+m[1]+m[2])*(1./3.));
- MFaceVertex *vmid = new MFaceVertex (gp.x(),gp.y(),gp.z(),_gf,gp.u(),gp.v());
+ double XX = (v[0]->x()+v[1]->x()+v[2]->x())*(1./3.);
+ double YY = (v[0]->y()+v[1]->y()+v[2]->y())*(1./3.);
+ double ZZ = (v[0]->z()+v[1]->z()+v[2]->z())*(1./3.);
+ MFaceVertex *vmid = new MFaceVertex (XX,YY,ZZ,_gf,gp.u(),gp.v());
_gf->mesh_vertices.push_back(vmid);
qnew.push_back(new MQuadrangle(_gf->triangles[i]->getVertex(0),v[0],vmid,v[2]));
qnew.push_back(new MQuadrangle(_gf->triangles[i]->getVertex(1),v[1],vmid,v[0]));
@@ -162,7 +171,10 @@ public:
m[j] = p1;
if (it == _middle.end() && it2 == eds.end()){
GPoint gp = _gf->point((p1+p2)*0.5);
- v[j] = new MFaceVertex (gp.x(),gp.y(),gp.z(),_gf,gp.u(),gp.v());
+ double XX = 0.5*(E.getVertex(0)->x()+E.getVertex(1)->x());
+ double YY = 0.5*(E.getVertex(0)->y()+E.getVertex(1)->y());
+ double ZZ = 0.5*(E.getVertex(0)->z()+E.getVertex(1)->z());
+ v[j] = new MFaceVertex (XX,YY,ZZ,_gf,gp.u(),gp.v());
_gf->mesh_vertices.push_back(v[j]);
eds[E] = v[j];
}
@@ -175,7 +187,11 @@ public:
}
}
GPoint gp = _gf->point((m[0]+m[1]+m[2]+m[3])*0.25);
- MVertex *vmid = new MFaceVertex (gp.x(),gp.y(),gp.z(),_gf,gp.u(),gp.v());
+ // FIXME : NOT EXACTLY CORRECT, BUT THAT'S THE PLACE WE WANT THE POINT TO RESIDE
+ double XX = 0.25*(v[0]->x()+v[1]->x()+v[2]->x()+v[3]->x());
+ double YY = 0.25*(v[0]->y()+v[1]->y()+v[2]->y()+v[3]->y());
+ double ZZ = 0.25*(v[0]->z()+v[1]->z()+v[2]->z()+v[3]->z());
+ MVertex *vmid = new MFaceVertex (XX,YY,ZZ,_gf,gp.u(),gp.v());
_gf->mesh_vertices.push_back(vmid);
qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(0),v[0],vmid,v[3]));
qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(1),v[1],vmid,v[0]));
@@ -188,17 +204,17 @@ public:
}
void finish ()
{
- backgroundMesh::setSizeFactor(1.0);
if((CTX::instance()->mesh.recombineAll || _gf->meshAttributes.recombine) &&
CTX::instance()->mesh.algoRecombine == 2){
// recombine the elements on the half mesh
- recombineIntoQuads(_gf,true,true,.1);
+ CTX::instance()->mesh.lcFactor /=2.0;
+ recombineIntoQuads(_gf,false,true,.1);
Msg::Info("subdividing");
- // _gf->model()->writeMSH("hop1.msh");
+ _gf->model()->writeMSH("hop1.msh");
subdivide();
- // _gf->model()->writeMSH("hop2.msh");
+ _gf->model()->writeMSH("hop2.msh");
restore();
- recombineIntoQuads(_gf,true,true,0.1);
+ recombineIntoQuads(_gf,false,true,0.1);
computeElementShapes(_gf,
_gf->meshStatistics.worst_element_shape,
_gf->meshStatistics.average_element_shape,
@@ -2367,7 +2383,14 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
void deMeshGFace::operator() (GFace *gf)
{
- if(gf->geomType() == GEntity::DiscreteSurface) return;
+ if(gf->geomType() == GEntity::DiscreteSurface) {
+ discreteFace *df = dynamic_cast<discreteFace*> (gf);
+ // copy the mesh and parametrize it to have a "geometry"
+ // of the discrete edge
+ // FIXME
+ return;
+ df->createGeometry();
+ }
gf->deleteMesh();
gf->meshStatistics.status = GFace::PENDING;
gf->meshStatistics.nbTriangle = gf->meshStatistics.nbEdge = 0;
@@ -2386,23 +2409,7 @@ void meshGFace::operator() (GFace *gf, bool print)
return;
}
- //-----------------------------------------------------------
- // FIXME PAB & JFR FIRST IMPLEMENTATION OF THE COMPOUND REPLACEMENT
- // if(gf->geomType() == GEntity::DiscreteSurface) {
- // discreteFace *df = dynamic_cast<discreteFace*> (gf);
- // if (df) {
- // df->createAtlas();
- // for (unsigned int i=0;i<df->_atlas.size();i++){
- // (*this)(df->_atlas[i],print);
- // }
- // df->gatherMeshes();
- // }
- // gf->meshStatistics.status = GFace::DONE;
- // return;
- // }
- //-----------------------------------------------------------
- //-----------------------------------------------------------
-
+ // if(gf->geomType() == GEntity::DiscreteFace) return;
if(gf->geomType() == GEntity::ProjectionFace) return;
if(gf->meshAttributes.method == MESH_NONE) return;
if(CTX::instance()->mesh.meshOnlyVisible && !gf->getVisibility()) return;
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 6500ecc58a3c703455d9bff2247651309f192688..85ca7063d1e921b79d186c7d561c08c6afc6eaf7 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1558,7 +1558,7 @@ void buildBackGroundMesh (GFace *gf,
int CurvControl = CTX::instance()->mesh.lcFromCurvature;
CTX::instance()->mesh.lcFromCurvature = 0;
// Do a background mesh
- bowyerWatson(gf,4000, equivalence,parametricCoordinates);
+ bowyerWatson(gf,40000, equivalence,parametricCoordinates);
// Re-enable curv control if asked
CTX::instance()->mesh.lcFromCurvature = CurvControl;
// apply this to the BGM
diff --git a/tutorial/t13.geo b/tutorial/t13.geo
index e102d9852f962476f4cd2a873d0e7beff87c0501..feebcb7c322fa9f6f7b25acdac2f474e8ac65187 100644
--- a/tutorial/t13.geo
+++ b/tutorial/t13.geo
@@ -39,6 +39,7 @@ Volume(1) = {1};
Physical Surface(1) = {s : s + #ss[]-1};
Physical Volume(1) = 1;
+Physical Line(1) = {26 ... 37};
// element size imposed by a size field
Field[1] = MathEval;