diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 0d0a618a6607c52f2abae9da1a663038969bed92..0fd1b5c576fe62d22cecd5b76bcdc7bdddc188b8 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -1,4 +1,4 @@
-// $Id: CommandLine.cpp,v 1.108 2007-11-26 14:34:09 remacle Exp $
+// $Id: CommandLine.cpp,v 1.109 2007-12-03 15:17:39 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -204,6 +204,10 @@ void Get_Options(int argc, char *argv[])
CTX.batch = 3;
i++;
}
+ else if(!strcmp(argv[i] + 1, "4")) {
+ CTX.batch = 4;
+ i++;
+ }
else if(!strcmp(argv[i] + 1, "pid")) {
fprintf(stdout, "%d\n", GetProcessId());
fflush(stdout);
diff --git a/Fltk/Main.cpp b/Fltk/Main.cpp
index d8f3a96e4252ac523d29f34cc7b9fd714f4a4c4d..8d1d2f829f5950c98fc0d0811d5d5a44759529d5 100644
--- a/Fltk/Main.cpp
+++ b/Fltk/Main.cpp
@@ -1,4 +1,4 @@
-// $Id: Main.cpp,v 1.112 2007-09-26 20:51:58 geuzaine Exp $
+// $Id: Main.cpp,v 1.113 2007-12-03 15:17:39 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -125,7 +125,11 @@ int main(int argc, char *argv[])
else
Msg(GERROR, "Invalid background mesh (no view)");
}
- if(CTX.batch > 0) {
+ if(CTX.batch == 4) {
+ AdaptMesh();
+ CreateOutputFile(CTX.output_filename, CTX.mesh.format);
+ }
+ else if(CTX.batch > 0) {
GenerateMesh(CTX.batch);
CreateOutputFile(CTX.output_filename, CTX.mesh.format);
}
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 3299ace23fe74ac69d480caf724fdb68c9d45579..e7e1e1b983d0b0dfa1b6b82d79c4c10b18dab2e2 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.23 2007-10-03 19:40:41 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.24 2007-12-03 15:17:40 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -20,6 +20,7 @@
// Please report all bugs and problems to <gmsh@geuz.org>.
#include <map>
+#include <cstring>
#include <string>
#include "Message.h"
diff --git a/Geo/GVertex.cpp b/Geo/GVertex.cpp
index ee938919107c6fed7291ef054d810f9eac019a02..6ce753e97114044141fa2eb8f2ad56ba9f8e4fb6 100644
--- a/Geo/GVertex.cpp
+++ b/Geo/GVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: GVertex.cpp,v 1.13 2007-09-03 12:00:28 geuzaine Exp $
+// $Id: GVertex.cpp,v 1.14 2007-12-03 15:17:40 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -20,6 +20,7 @@
// Please report all bugs and problems to <gmsh@geuz.org>.
#include <algorithm>
+#include <cstring>
#include "GVertex.h"
#include "GFace.h"
#include "Message.h"
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 1b4630561e5475162585cfc193f9544f488b137f..a514effaf605cd72a02d0c6833239de27925dad6 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -304,7 +304,7 @@ class MLineN : public MLine {
class MTriangle : public MElement {
protected:
MVertex *_v[3];
- virtual void jac(int order, MVertex *v[], double u, double v, double jac[2][2]);
+ virtual void jac(int order, MVertex *verts[], double u, double v, double jac[2][2]);
public :
MTriangle(MVertex *v0, MVertex *v1, MVertex *v2, int num=0, int part=0)
: MElement(num, part)
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 9302a3e2a5b911f0f596cfe91cb4f20b882cc612..0c7532683e6739bae793cdbd913d0ab6fee2440d 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: MVertex.cpp,v 1.16 2007-09-21 16:22:51 geuzaine Exp $
+// $Id: MVertex.cpp,v 1.17 2007-12-03 15:17:40 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -19,7 +19,8 @@
//
// Please report all bugs and problems to <gmsh@geuz.org>.
-#include <math.h>
+#include <cstring>
+#include <cmath>
#include "MVertex.h"
#include "GEdge.h"
#include "GFace.h"
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 77779abe43ff25424b5d256e37465268498092e3..95223af5c84936f8e29b7a2f392da5d131c68c59 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCEdge.cpp,v 1.27 2007-11-26 14:34:09 remacle Exp $
+// $Id: OCCEdge.cpp,v 1.28 2007-12-03 15:17:40 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -190,14 +190,15 @@ GEntity::GeomType OCCEdge::geomType() const
int OCCEdge::minimumMeshSegments() const
{
- int np;
- if(geomType() == Line)
- np= GEdge::minimumMeshSegments();
- else
+ int np;
+ if(geomType() == Line)
+ np= GEdge::minimumMeshSegments();
+ else
np=CTX.mesh.min_curv_points - 1;
- // if the edge is closed, ensure that at least 3 points are generated
- if (getBeginVertex() == getEndVertex()) np=std::max(2,np);
+ // if the edge is closed, ensure that at least 3 points are generated in the
+ // 1D mesh (4 segments, one of which is degenerated)
+ if (getBeginVertex() == getEndVertex()) np=std::max(4,np);
return np;
}
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index a64bac22a729cfd6a51e24dee73eb9e0402527ef..c6b86fa80ba00e23060e640799b0383a6cdfa92c 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.127 2007-11-28 14:18:10 remacle Exp $
+// $Id: Generator.cpp,v 1.128 2007-12-03 15:17:40 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -331,6 +331,36 @@ void OptimizeMesh(GModel *m)
Msg(STATUS1, "Mesh");
}
+void AdaptMesh()
+{
+ Msg(STATUS1, "Adapting the 3D Mesh...");
+ double t1 = Cpu();
+
+ if(CTX.threads_lock) {
+ Msg(INFO, "I'm busy! Ask me that later...");
+ return;
+ }
+
+ CTX.threads_lock = 1;
+
+ GModel *m = GModel::current();
+
+ std::for_each(m->firstRegion(), m->lastRegion(), adaptMeshGRegion());
+ std::for_each(m->firstRegion(), m->lastRegion(), adaptMeshGRegion());
+ std::for_each(m->firstRegion(), m->lastRegion(), adaptMeshGRegion());
+ std::for_each(m->firstRegion(), m->lastRegion(), adaptMeshGRegion());
+ std::for_each(m->firstRegion(), m->lastRegion(), adaptMeshGRegion());
+ std::for_each(m->firstRegion(), m->lastRegion(), adaptMeshGRegion());
+ std::for_each(m->firstRegion(), m->lastRegion(), adaptMeshGRegion());
+ std::for_each(m->firstRegion(), m->lastRegion(), adaptMeshGRegion());
+ std::for_each(m->firstRegion(), m->lastRegion(), adaptMeshGRegion());
+ std::for_each(m->firstRegion(), m->lastRegion(), adaptMeshGRegion());
+
+ double t2 = Cpu();
+ Msg(INFO, "Mesh Adaptation complete (%g s)", t2 - t1);
+ Msg(STATUS1, "Mesh");
+}
+
void GenerateMesh(int ask)
{
if(CTX.threads_lock) {
diff --git a/Mesh/Generator.h b/Mesh/Generator.h
index a39d30ff9ed0ab5a2ab8b69449e3eed31d358061..4bcdc68156146a176417007a88d69362960c91e7 100644
--- a/Mesh/Generator.h
+++ b/Mesh/Generator.h
@@ -23,6 +23,7 @@
class GModel;
void GetStatistics(double stat[50], double quality[3][100]=0);
+void AdaptMesh();
void GenerateMesh(int dimension);
void OptimizeMesh(GModel *m);
void OptimizeMeshNetgen(GModel *m);
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index ae3f8a89c6ef54eceb2c669e54a0560a866367c3..5d307bc3f6a41e01144d2d16e60bea4cb6557bb8 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -20,6 +20,7 @@
//
// Please report all bugs and problems to <gmsh@geuz.org>.
+class GModel;
#include "MElement.h"
#include <list>
#include <set>
@@ -132,5 +133,4 @@ struct edgeXface
};
-
#endif
diff --git a/Mesh/meshGRegion.h b/Mesh/meshGRegion.h
index 6744f4bf6f4e035da2a509e70193fd37a040a685..6f3fb49c7bfe205f396ea61558d6f492431d1237 100644
--- a/Mesh/meshGRegion.h
+++ b/Mesh/meshGRegion.h
@@ -56,6 +56,13 @@ class deMeshGRegion {
void operator () (GRegion *);
};
+// adapt the mesh of a region
+class adaptMeshGRegion {
+ public :
+ void operator () (GRegion *);
+};
+
+
void MeshDelaunayVolume(std::vector<GRegion*> &delaunay);
int MeshTransfiniteVolume(GRegion *gr);
int SubdivideExtrudedMesh(GModel *m);
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 16c8064b82a1fa062da6cf7a90cb1ef71431d512..42dffd8cad07a94903a749bd4fe40a6ede67a400 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.24 2007-11-28 14:18:10 remacle Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.25 2007-12-03 15:17:40 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -25,6 +25,7 @@
#include "meshGRegionDelaunayInsertion.h"
#include "GModel.h"
#include "GRegion.h"
+#include "meshGRegion.h"
#include "Numeric.h"
#include "Message.h"
#include <set>
@@ -324,7 +325,6 @@ static void setLcs ( MTetrahedron *t, std::map<MVertex*,double> &vSizes)
}
-typedef std::multimap<MVertex*,std::pair<MTriangle*,GFace*> > fs_cont ;
GFace* findInFaceSearchStructure ( MVertex *p1,MVertex *p2,MVertex *p3, const fs_cont&search )
@@ -432,9 +432,10 @@ void recur_classify ( MTet4 *t ,
}
}
-//template <class CONTAINER, class DATA>
-void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
+void adaptMeshGRegion::operator () (GRegion *gr)
{
+ const gmshQualityMeasure4Tet qm = QMTET_2;
+
typedef std::list<MTet4 *> CONTAINER ;
CONTAINER allTets;
for (int i=0;i<gr->tetrahedra.size();i++){
@@ -444,8 +445,180 @@ void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
connectTets(allTets.begin(),allTets.end());
+ double t1 = Cpu();
+ std::vector<MTet4*> illegals;
+ const int nbRanges=10;
+ int quality_ranges [nbRanges];
+ {
+ double totalVolumeb = 0.0;
+ double worst = 1.0;
+ double avg = 0;
+ int count=0;
+ for (int i=0;i<nbRanges;i++)quality_ranges[i] = 0;
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
+ {
+ if (!(*it)->isDeleted()){
+ double vol = fabs((*it)->tet()->getVolume());
+ double qual = (*it)->getQuality();
+ worst = std::min(qual,worst);
+ avg+=qual;
+ count++;
+ totalVolumeb+=vol;
+ for (int i=0;i<nbRanges;i++){
+ double low = (double)i/(nbRanges);
+ double high = (double)(i+1)/(nbRanges);
+ if (qual >= low && qual < high)quality_ranges[i]++;
+ }
+ }
+ }
+ Msg(INFO,"Opti : START with %12.5E QBAD %12.5E QAVG %12.5E",totalVolumeb,worst,avg/count);
+ for (int i=0;i<nbRanges;i++){
+ double low = (double)i/(nbRanges);
+ double high = (double)(i+1)/(nbRanges);
+ Msg(INFO,"Opti : %3.2f < QUAL < %3.2f : %9d elements ",low,high,quality_ranges[i]);
+ }
+ }
+
+ double qMin = 0.5;
+ double sliverLimit = 0.2;
+ int nbESwap=0, nbFSwap=0, nbReloc=0, nbCollapse = 0;
+
+ while (1){
+ std::vector<MTet4*> newTets;
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it){
+ if (!(*it)->isDeleted()){
+ for (int i=0;i<4;i++){
+ for (int j=0;j<4;j++){
+ if (gmshCollapseVertex(newTets,*it,i,j,QMTET_2)){nbCollapse++;i=j=10;}
+ }
+ }
+ }
+ }
+
+ printf("nbCollapses = %d\n",nbCollapse);
+
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it){
+ if (!(*it)->isDeleted()){
+ double vol;
+ double qq = (*it)->getQuality();
+ if (qq < qMin){
+ for (int i=0;i<4;i++){
+ if (gmshFaceSwap(newTets,*it,i,qm)){nbFSwap++;break;}
+ }
+ }
+ }
+ }
+
+ illegals.clear();
+ for (int i=0;i<nbRanges;i++)quality_ranges[i] = 0;
+
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
+ {
+ if (!(*it)->isDeleted()){
+ double qq = (*it)->getQuality();
+ if (qq < qMin)
+ for (int i=0;i<6;i++){
+ if (gmshEdgeSwap(newTets,*it,i,qm)) {nbESwap++;break;}
+ }
+ if (!(*it)->isDeleted()){
+ if (qq < sliverLimit)illegals.push_back(*it);
+ for (int i=0;i<nbRanges;i++){
+ double low = (double)i/(nbRanges);
+ double high = (double)(i+1)/(nbRanges);
+ if (qq >= low && qq < high)quality_ranges[i]++;
+ }
+ }
+ }
+ }
+
+ if (!newTets.size())break;
+
+ // add all the new tets in the container
+ for (int i=0;i<newTets.size();i++){
+ if (!newTets[i]->isDeleted()){
+ allTets.push_back(newTets[i]);
+ }
+ else{
+ delete newTets[i]->tet();
+ delete newTets[i];
+ }
+ }
+
+ // relocate vertices
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
+ {
+ if (!(*it)->isDeleted()){
+ double qq = (*it)->getQuality();
+ if (qq < qMin)
+ for (int i=0;i<4;i++){
+ if (gmshSmoothVertex(*it,i,qm))nbReloc++;
+ }
+ }
+ }
+
+ double totalVolumeb = 0.0;
+ double worst = 1.0;
+ double avg = 0;
+ int count=0;
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
+ {
+ if (!(*it)->isDeleted()){
+ double vol = fabs((*it)->tet()->getVolume());
+ double qual = (*it)->getQuality();
+ worst = std::min(qual,worst);
+ avg+=qual;
+ count++;
+ totalVolumeb+=vol;
+ }
+ }
+ double t2 = Cpu();
+ Msg(INFO,"Opti : (%d,%d,%d) = %12.5E QBAD %12.5E QAVG %12.5E (%8.3f sec)",nbESwap,nbFSwap,nbReloc,totalVolumeb,worst,avg/count,t2-t1);
+ break;
+ }
+
+ int nbSlivers = 0;
+ for (int i=0;i<illegals.size();i++)
+ if (!(illegals[i]->isDeleted()))nbSlivers ++;
+
+ if (nbSlivers){
+ Msg(INFO,"Opti : %d illegal tets are still in the mesh, trying to remove them",nbSlivers);
+ }
+ else{
+ Msg(INFO,"Opti : no illegal tets in the mesh ;-)",nbSlivers);
+ }
+
+ for (int i=0;i<nbRanges;i++){
+ double low = (double)i/(nbRanges);
+ double high = (double)(i+1)/(nbRanges);
+ Msg(INFO,"Opti : %3.2f < QUAL < %3.2f : %9d elements ",low,high,quality_ranges[i]);
+ }
+
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
+ {
+ if (!(*it)->isDeleted()){
+ gr->tetrahedra.push_back((*it)->tet());
+ delete *it;
+ }
+ else{
+ delete (*it)->tet();
+ delete *it;
+ }
+ }
+}
+
+//template <class CONTAINER, class DATA>
+void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
+{
+ typedef std::list<MTet4 *> CONTAINER ;
+ CONTAINER allTets;
+ for (int i=0;i<gr->tetrahedra.size();i++){
+ allTets.push_back(new MTet4(gr->tetrahedra[i],qm));
+ }
+ gr->tetrahedra.clear();
+
+ connectTets(allTets.begin(),allTets.end());
double t1 = Cpu();
std::vector<MTet4*> illegals;
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index d10f0c91dea908b26f800db6e62f5fdf48e5e29e..1c6469b5357e188f9abe373e21115bdaa225fdcf 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -24,10 +24,13 @@
#include "qualityMeasures.h"
#include <list>
#include <set>
+#include <map>
#include <stack>
//#define _GMSH_PRE_ALLOCATE_STRATEGY_ 1
class GRegion;
+class GFace;
+class GModel;
class MTet4Factory;
@@ -232,5 +235,8 @@ private:
};
void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm);
+typedef std::multimap<MVertex*,std::pair<MTriangle*,GFace*> > fs_cont ;
+GFace* findInFaceSearchStructure ( MVertex *p1,MVertex *p2,MVertex *p3, const fs_cont&search );
+bool buildFaceSearchStructure ( GModel *model , fs_cont&search );
#endif
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index e0017a26b59799b14cf863b30bf9b99312a4496c..c2066916a715dcfac70098e979be01c92457f982 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -1,5 +1,6 @@
#include "meshGRegionLocalMeshMod.h"
#include "GEntity.h"
+#include "Message.h"
static int edges[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
static int efaces[6][2] = {{0,2},{0,1},{1,2},{0,3},{2,3},{1,3}};
@@ -13,6 +14,35 @@ static int vFac[4][3] = {{0,1,2},{0,2,3},{0,1,3},{1,2,3}};
// all tets that share this edge and all vertices that are
// forming the outer ring of the cavity
// we return true if the cavity is closed and false if it is open
+
+void computeNeighboringTetsOfACavity (const std::vector<MTet4*> &cavity,std::vector<MTet4*> &outside)
+{
+ outside.clear();
+ for (int i=0;i<cavity.size();i++){
+ for (int j=0;j<4;j++){
+ MTet4 * neigh = cavity[i]->getNeigh(j);
+ if (neigh)
+ {
+ bool found = false;
+ for (int k=0;k<outside.size();k++){
+ if(outside[k]==neigh){
+ found=true;
+ break;
+ }
+ }
+ if (!found){
+ for (int k=0;k<cavity.size() ;k++){
+ if(cavity[k] ==neigh){
+ found=true;
+ }
+ }
+ }
+ if(!found)outside.push_back(neigh);
+ }
+ }
+ }
+}
+
bool gmshBuildEdgeCavity ( MTet4 *t,
int iLocalEdge,
MVertex **v1,MVertex **v2,
@@ -21,7 +51,6 @@ bool gmshBuildEdgeCavity ( MTet4 *t,
std::vector<MVertex*> &ring ){
cavity.clear();
ring.clear();
- outside.clear();
*v1 = t->tet()->getVertex(edges[iLocalEdge][0]);
*v2 = t->tet()->getVertex(edges[iLocalEdge][1]);
@@ -68,31 +97,7 @@ bool gmshBuildEdgeCavity ( MTet4 *t,
throw;
}
}
-
- for (int i=0;i<cavity.size();i++){
- for (int j=0;j<4;j++){
- MTet4 * neigh = cavity[i]->getNeigh(j);
- if (neigh)
- {
- bool found = false;
- for (int k=0;k<outside.size();k++){
- if(outside[k]==neigh){
- found=true;
- break;
- }
- }
- if (!found){
- for (int k=0;k<cavity.size() ;k++){
- if(cavity[k] ==neigh){
- found=true;
- }
- }
- }
- if(!found)outside.push_back(neigh);
- }
- }
- }
-
+ computeNeighboringTetsOfACavity (cavity,outside);
return true;
}
@@ -319,6 +324,47 @@ bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
return true;
}
+bool gmshEdgeSplit (std::vector<MTet4 *> &newTets,
+ MTet4 *tet,
+ MVertex *newVertex,
+ int iLocalEdge,
+ const gmshQualityMeasure4Tet &cr){
+ std::vector<MTet4*> cavity;
+ std::vector<MTet4*> outside;
+ std::vector<MVertex*> ring;
+ MVertex *v1,*v2;
+
+ bool closed = gmshBuildEdgeCavity ( tet,iLocalEdge,&v1,&v2,cavity,outside,ring);
+ if (!closed)return false;
+
+ for (int j=0;j<ring.size();j++){
+ MVertex *pv1 = ring[j];
+ MVertex *pv2 = ring[(j+1)%ring.size()];
+ MTetrahedron *tr1 = new MTetrahedron ( pv1,
+ pv2,
+ newVertex,
+ v1);
+ MTetrahedron *tr2 = new MTetrahedron ( newVertex,
+ pv2,
+ pv1,
+ v2);
+ MTet4 *t41 = new MTet4 ( tr1 , cr) ;
+ MTet4 *t42 = new MTet4 ( tr2 , cr);
+ t41->setOnWhat(cavity[0]->onWhat());
+ t42->setOnWhat(cavity[0]->onWhat());
+ outside.push_back(t41);
+ outside.push_back(t42);
+ newTets.push_back(t41);
+ newTets.push_back(t42);
+ }
+
+ for (int i=0;i<cavity.size();i++)cavity[i]->setDeleted(true);
+
+ connectTets ( outside );
+
+ return true;
+}
+
// swap a face i.e. remove a face shared by 2 tets
bool gmshFaceSwap (std::vector<MTet4 *> &newTets,
MTet4 *t1,
@@ -413,7 +459,7 @@ void gmshBuildVertexCavity_recur ( MTet4 *t,
std::vector<MTet4*> &cavity){
// if (recur > 20)printf("oufti %d\n",recur);
if(t->isDeleted()){
- printf("a deleted triangle is a neighbor of a non deleted triangle\n");
+ Msg(FATAL,"a deleted triangle is a neighbor of a non deleted triangle");
}
int iV=-1;
for (int i=0; i<4; i++){
@@ -423,7 +469,7 @@ void gmshBuildVertexCavity_recur ( MTet4 *t,
}
}
if (iV==-1){
- printf("this tet does not contain a given vertex\n");
+ Msg(FATAL,"trying to build a cavity of tets for a vertex that does not belong to this tet");
}
for (int i=0; i<3; i++){
MTet4 *neigh = t->getNeigh(vFac[iV][i]);
@@ -444,6 +490,79 @@ void gmshBuildVertexCavity_recur ( MTet4 *t,
}
+bool gmshCollapseVertex ( std::vector<MTet4 *> &newTets,
+ MTet4 *t,
+ int iVertex,
+ int iTarget,
+ const gmshQualityMeasure4Tet &cr){
+
+ if(t->isDeleted())throw;
+
+ MVertex *v = t->tet()->getVertex(iVertex);
+ MVertex *tg = t->tet()->getVertex(iTarget);
+
+ if(v->onWhat()->dim() < 3)return false;
+ if(tg->onWhat()->dim() < 3)return false;
+
+ std::vector<MTet4*> cavity_v;
+ std::vector<MTet4*> outside;
+ cavity_v.push_back(t);
+ gmshBuildVertexCavity_recur (t,v,cavity_v);
+
+ std::vector<MTet4*> toDelete;
+ std::vector<MTet4*> toUpdate;
+ double volume = 0;
+ for (int i=0;i<cavity_v.size();i++){
+ bool found = false;
+ volume+=fabs(cavity_v[i]->tet()->getVolume());
+ for (int j=0;j<4;j++){
+ if (cavity_v[i]->tet()->getVertex(j) == tg)found=true;
+ }
+ if (found) toDelete.push_back(cavity_v[i]);
+ else toUpdate.push_back(cavity_v[i]);
+ }
+
+ double x = v->x();
+ double y = v->y();
+ double z = v->z();
+ v->x() = tg->x();
+ v->y() = tg->y();
+ v->z() = tg->z();
+
+ double volume_update=0;
+ for (int i=0;i<toUpdate.size();i++){
+ volume_update+=fabs(toUpdate[i]->tet()->getVolume());
+ }
+
+ // printf("%12.5E %12.5E\n",volume,volume_update);
+
+ if (fabs(volume-volume_update) > 1.e-10 * volume)
+ {
+ v->x() = x;
+ v->y() = y;
+ v->z() = z;
+ return false;
+ }
+ // ok we collapse
+ computeNeighboringTetsOfACavity (cavity_v,outside);
+ for (int i=0;i<toUpdate.size();i++){
+ MTetrahedron *tr1 = new MTetrahedron ( toUpdate[i]->tet()->getVertex(0) == v ? tg : toUpdate[i]->tet()->getVertex(0),
+ toUpdate[i]->tet()->getVertex(1) == v ? tg : toUpdate[i]->tet()->getVertex(1),
+ toUpdate[i]->tet()->getVertex(2) == v ? tg : toUpdate[i]->tet()->getVertex(2),
+ toUpdate[i]->tet()->getVertex(3) == v ? tg : toUpdate[i]->tet()->getVertex(3));
+ MTet4 *t41 = new MTet4 ( tr1 , cr) ;
+ t41->setOnWhat(cavity_v[0]->onWhat());
+ outside.push_back(t41);
+ newTets.push_back(t41);
+ }
+ for (int i=0;i<cavity_v.size();i++)cavity_v[i]->setDeleted(true);
+
+ connectTets ( outside );
+
+ return true;
+}
+
+
bool gmshSmoothVertex ( MTet4 *t,
int iVertex,
const gmshQualityMeasure4Tet &cr){
@@ -517,3 +636,303 @@ bool gmshSmoothVertex ( MTet4 *t,
}
+// Edge split sets ...
+// Here, we only build a list of tets that are a subdivision
+// of the given
+//static int edges[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
+
+// the resulting triangles are 012 and 230 in vector res
+// void splitPrism (std::vector<MTet4 *> &newTets,
+// std::vector<MVertex *> &steinerPoints,
+// MVertex* v1,
+// MVertex* v2,
+// MVertex* v3,
+// MVertex* v4,
+// MVertex* v5,
+// MVertex* v6,
+// const gmshQualityMeasure4Tet &cr,
+// MVertex *res[4],
+// fs_cont &search,
+// GFace *fakeFace){
+// }
+
+// void splitQuadFace (MVertex* newv1,
+// MVertex* newv2,
+// MVertex* v21,
+// MVertex* v11,
+// const gmshQualityMeasure4Tet &cr,
+// MVertex *res[4],
+// fs_cont &search,
+// GFace *fakeFace){
+// GFace* gfound = findInFaceSearchStructure (newv1,newv2,v11,search);
+// if (gfound){
+// res[0] = newv2;
+// res[1] = newv1;
+// res[2] = v11;
+// res[3] = v21;
+// }
+// else{
+// GFace* gfound = findInFaceSearchStructure (newv1,newv2,v21,search);
+// if (gfound){
+// res[0] = newv1;
+// res[1] = newv2;
+// res[2] = v21;
+// res[3] = v11;
+// }
+// // choose the best configuration
+// else{
+// double q1 = std::min(qmTri(newv1,newv2,v11,cr),qmTet(newv2,v11,v21,cr));
+// double q2 = std::min(qmTet(newv1,newv2,v21,cr),qmTet(newv1,v11,v21,cr));
+// if (q1 > q2){
+// res[0] = newv2;
+// res[1] = newv1;
+// res[2] = v11;
+// res[3] = v21;
+// MVertex *p1 = std::min(std::min(newv1,newv2),v11);
+// MVertex *p2 = std::min(std::min(newv2,v11),v21);
+// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p1, std::pair<MTriangle*,GFace*>(new MTriangle(newv1,newv2,v11),fakeFace)));
+// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p2, std::pair<MTriangle*,GFace*>(new MTriangle(newv2,v11,v21),fakeFace)));
+// }
+// else{
+// res[0] = newv1;
+// res[1] = newv2;
+// res[2] = v21;
+// res[3] = v11;
+// MVertex *p1 = std::min(std::min(newv1,newv2),v21);
+// MVertex *p2 = std::min(std::min(newv2,v11),v21);
+// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p1, std::pair<MTriangle*,GFace*>(new MTriangle(newv1,newv2,v21),fakeFace)));
+// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p2, std::pair<MTriangle*,GFace*>(new MTriangle(newv1,v11,v21),fakeFace)));
+// }
+// }
+// }
+// }
+// void splitQuadFace (std::vector<MTet4 *> &newTets,
+// std::vector<MVertex *> &steinerPoints,
+// MVertex* newv1,
+// MVertex* newv2,
+// MVertex* v21,
+// MVertex* v11,
+// MVertex* other,
+// const gmshQualityMeasure4Tet &cr,
+// fs_cont &search,
+// GFace *fakeFace){
+// MTetrahedron *tr3,*tr4;
+// // now we look if there exists a face with the same vertices
+// GFace* gfound = findInFaceSearchStructure (newv1,newv2,v11,search);
+// if (gfound){
+// tr3 = new MTetrahedron ( newv1,newv2,v11,other);
+// tr4 = new MTetrahedron ( newv2,v11,v21,other);
+// }
+// else{
+// GFace* gfound = findInFaceSearchStructure (newv1,newv2,v21,search);
+// if (gfound){
+// tr3 = new MTetrahedron ( newv1,newv2,v21,other);
+// tr4 = new MTetrahedron ( newv1,v11,v21,other);
+// }
+// // choose the best configuration
+// else{
+// double vol;
+// double q1 = std::min(qmTet(newv1,newv2,v11,other,cr,vol),qmTet(newv2,v11,v21,other,cr,vol));
+// double q2 = std::min(qmTet(newv1,newv2,v21,other,cr,vol),qmTet(newv1,v11,v21,other,cr,vol));
+// if (q1 > q2){
+// tr3 = new MTetrahedron ( newv1,newv2,v11,other);
+// tr4 = new MTetrahedron ( newv2,v11,v21,other);
+// MVertex *p1 = std::min(std::min(newv1,newv2),v11);
+// MVertex *p2 = std::min(std::min(newv2,v11),v21);
+// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p1, std::pair<MTriangle*,GFace*>(new MTriangle(newv1,newv2,v11),fakeFace)));
+// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p2, std::pair<MTriangle*,GFace*>(new MTriangle(newv2,v11,v21),fakeFace)));
+// }
+// else{
+// tr3 = new MTetrahedron ( newv1,newv2,v21,other);
+// tr4 = new MTetrahedron ( newv1,v11,v21,other);
+// MVertex *p1 = std::min(std::min(newv1,newv2),v21);
+// MVertex *p2 = std::min(std::min(newv2,v11),v21);
+// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p1, std::pair<MTriangle*,GFace*>(new MTriangle(newv1,newv2,v21),fakeFace)));
+// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p2, std::pair<MTriangle*,GFace*>(new MTriangle(newv1,v11,v21),fakeFace)));
+// }
+// }
+// }
+// }
+
+// bool splitEdgesOfTet (std::vector<MTet4 *> &newTets,
+// std::vector<MVertex *> &steinerPoints,
+// MTet4 *t1,
+// int nbEdges,
+// int e[6],
+// MVertex* pts[6],
+// const gmshQualityMeasure4Tet &cr,
+// fs_cont &search,
+// GFace *fakeFace){
+// switch(nbEdges){
+// case 1 :
+// {
+// int iE = edges[0];
+// MVertex *newv = pts[0];
+// MVertex *v1 = t1->tet()->getVertex(e[iE][0]);
+// MVertex *v2 = t1->tet()->getVertex(e[iE][1]);
+// int oE = 5-iE;
+// MVertex *o1 = t1->tet()->getVertex(e[oE][0]);
+// MVertex *o2 = t1->tet()->getVertex(e[oE][1]);
+// MTetrahedron *tr1 = new MTetrahedron ( newv,v1,o1,o2);
+// MTetrahedron *tr2 = new MTetrahedron ( newv,v2,o2,o1);
+// MTet4 *t41 = new MTet4 (tr1,cr);
+// MTet4 *t42 = new MTet4 (tr2,cr);
+// newTets.push_back(t41);
+// newTets.push_back(t42);
+// }
+// break;
+// case 2 :
+// {
+// int iE1 = edges[0];
+// int iE2 = edges[1];
+
+// MVertex *newv1 = pts[0];
+// MVertex *newv2 = pts[1];
+
+// MVertex *v11 = t1->tet()->getVertex(e[iE1][0]);
+// MVertex *v12 = t1->tet()->getVertex(e[iE1][1]);
+// MVertex *v21 = t1->tet()->getVertex(e[iE2][0]);
+// MVertex *v22 = t1->tet()->getVertex(e[iE2][1]);
+
+// if (iE1 == 5-iE2){ // two opposite edges
+// MTetrahedron *tr1 = new MTetrahedron ( newv1,newv2,v11,v21);
+// MTetrahedron *tr2 = new MTetrahedron ( newv1,newv2,v21,v12);
+// MTetrahedron *tr3 = new MTetrahedron ( newv1,newv2,v12,v22);
+// MTetrahedron *tr4 = new MTetrahedron ( newv1,newv2,v22,v11);
+// MTet4 *t41 = new MTet4 (tr1,cr);
+// MTet4 *t42 = new MTet4 (tr2,cr);
+// MTet4 *t43 = new MTet4 (tr3,cr);
+// MTet4 *t44 = new MTet4 (tr4,cr);
+// newTets.push_back(t41);
+// newTets.push_back(t42);
+// newTets.push_back(t43);
+// newTets.push_back(t44);
+// }
+// else{ //two adjacend edges
+// MVertex *vsame,*other=0;
+// if (v11 == v21){vsame=v11; v11=v12; v21=v22;}
+// else if (v11 == v22){vsame=v11; v11=v12}
+// else if (v12 == v21){vsame=v12; v21=v22}
+// else if (v12 == v22){vsame=v12;}
+// else throw;
+// for (int i=0;i<4;i++){
+// if (vsame != t1->tet()->getVertex(i) &&
+// v11 != t1->tet()->getVertex(i) &&
+// v21 != t1->tet()->getVertex(i)){
+// other = t1->tet->getVertex(i);
+// break;
+// }
+// }
+// if (!other)throw;
+// MTetrahedron *tr1 = new MTetrahedron ( newv1,newv2,vsame,other);
+// splitQuadFace (newTets,steinerPoints,newv1,newv2,v21,v11,other,cr,search,fakeFace);
+// }
+// }
+// break;
+// case 3 :
+// {
+// int iE1 = edges[0];
+// int iE2 = edges[1];
+// int iE3 = edges[2];
+// MVertex *newv1;
+// MVertex *newv2;
+// MVertex *newv3;
+// std::sort(edges,edges+3);
+// if (iE1 == edges[0])newv1=pts[0];
+// else if (iE2 == edges[0])newv1=pts[1];
+// else if (iE3 == edges[0])newv1=pts[2];
+// if (iE1 == edges[1])newv2=pts[0];
+// else if (iE2 == edges[1])newv2=pts[1];
+// else if (iE3 == edges[1])newv2=pts[2];
+// if (iE1 == edges[2])newv3=pts[0];
+// else if (iE2 == edges[2])newv3=pts[1];
+// else if (iE3 == edges[2])newv3=pts[2];
+// iE1 = edges[0];
+// iE2 = edges[1];
+// iE3 = edges[2];
+// // edges are sorted and points correspond
+
+// mVertex *v[4];
+// // the 3 edges coincide at one vertex
+// int config;
+// if (iE1 == 0 && iE2 == 1 && iE3 == 2){
+// config = 1;
+// v[0] = t1->tet()->getVertex(0);
+// v[1] = t1->tet()->getVertex(1);
+// v[2] = t1->tet()->getVertex(2);
+// v[3] = t1->tet()->getVertex(3);
+// }
+// else if (iE1 == 0 && iE2 == 3 && iE3 == 4){
+// config = 1;
+// v[0] = t1->tet()->getVertex(1);
+// v[1] = t1->tet()->getVertex(0);
+// v[2] = t1->tet()->getVertex(2);
+// v[3] = t1->tet()->getVertex(3);
+// }
+// else if (iE1 == 1 && iE2 == 3 && iE3 == 5){
+// config = 1;
+// v[0] = t1->tet()->getVertex(2);
+// v[1] = t1->tet()->getVertex(0);
+// v[2] = t1->tet()->getVertex(1);
+// v[3] = t1->tet()->getVertex(3);
+// }
+// else if (iE1 == 2 && iE2 == 4 && iE3 == 5){
+// config = 1;
+// v[0] = t1->tet()->getVertex(3);
+// v[1] = t1->tet()->getVertex(0);
+// v[2] = t1->tet()->getVertex(1);
+// v[3] = t1->tet()->getVertex(2);
+// }
+// // three edges of the same face are cut
+// else if (iE1 == 0 && iE2 == 1 && iE3 == 3){
+// config = 2;
+// v[0] = t1->tet()->getVertex(3);
+// v[1] = t1->tet()->getVertex(1);
+// v[2] = t1->tet()->getVertex(0);
+// v[3] = t1->tet()->getVertex(2);
+// }
+// else if (iE1 == 0 && iE2 == 2 && iE3 == 4){
+// config = 2;
+// v[0] = t1->tet()->getVertex(2);
+// v[1] = t1->tet()->getVertex(1);
+// v[2] = t1->tet()->getVertex(0);
+// v[3] = t1->tet()->getVertex(3);
+// }
+// else if (iE1 == 1 && iE2 == 2 && iE3 == 5){
+// config = 2;
+// v[0] = t1->tet()->getVertex(1);
+// v[1] = t1->tet()->getVertex(2);
+// v[2] = t1->tet()->getVertex(0);
+// v[3] = t1->tet()->getVertex(3);
+// }
+// else if (iE1 == 3 && iE2 == 4 && iE3 == 5){
+// config = 2;
+// v[0] = t1->tet()->getVertex(0);
+// v[1] = t1->tet()->getVertex(2);
+// v[2] = t1->tet()->getVertex(1);
+// v[3] = t1->tet()->getVertex(3);
+// }
+// // the three edges for a kind of Z
+// else if (iE1 == 0 && iE2 == 3 && iE3 == 5){
+// config = 3;
+// }
+
+// if (config == 1){
+// }
+// else if (config == 2){
+// }
+// else{
+// throw; // for the moment.
+// }
+// }
+// default :
+// throw;
+// }
+// t1->setDeleted(true);
+
+
+// }
+
+
+
+// collapse
diff --git a/Mesh/meshGRegionLocalMeshMod.h b/Mesh/meshGRegionLocalMeshMod.h
index 7a4ef0dfc1369341c92a0083879928ca2f8085b9..1342d6038132bc052b0b920a7accd6e0e962c02f 100644
--- a/Mesh/meshGRegionLocalMeshMod.h
+++ b/Mesh/meshGRegionLocalMeshMod.h
@@ -22,19 +22,36 @@
#include "meshGRegionDelaunayInsertion.h"
#include "qualityMeasures.h"
-
+// local mesh modification operators. Those
+// operators only apply to the "bulk" of the
+// mesh and cannot be applied to boudnaries.
+// I'm working on it
bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
MTet4 *tet,
int iLocalEdge,
const gmshQualityMeasure4Tet &cr);
+
bool gmshFaceSwap (std::vector<MTet4 *> &newTets,
MTet4 *tet,
int iLocalFace,
const gmshQualityMeasure4Tet &cr);
+
bool gmshSmoothVertex ( MTet4 *t,
int iLocalVertex,
const gmshQualityMeasure4Tet &cr);
-
+
+bool gmshCollapseVertex ( std::vector<MTet4 *> &newTets,
+ MTet4 *t,
+ int iVertex,
+ int iTarget,
+ const gmshQualityMeasure4Tet &cr);
+
+bool gmshEdgeSplit (std::vector<MTet4 *> &newTets,
+ MTet4 *tet,
+ MVertex *newVertex,
+ int iLocalEdge,
+ const gmshQualityMeasure4Tet &cr);
+
#endif