new 3D mesher has optimization

parent a67ac8e8
Pipeline #162 passed with stage
in 9 minutes 1 second
......@@ -20,13 +20,17 @@
#endif
typedef unsigned char CHECKTYPE;
struct Tet;
struct Vert {
private :
double _x[3];
double _lc;
unsigned int _num;
Tet *_t;
public :
inline void setT(Tet*t) { _t = t;}
inline Tet* getT() const { return _t;}
inline unsigned int getNum () const {return _num;}
inline void setNum (unsigned int n) {_num=n;}
unsigned char _thread;
......@@ -40,7 +44,7 @@ public :
inline double &lc() {return _lc;}
inline operator double *() { return _x; }
Vert (double X=0, double Y=0, double Z=0, double lc=0, int num = 0)
: _num(num), _thread(0)
: _num(num),_t(NULL), _thread(0)
{
_x[0] = X; _x[1] = Y; _x[2] = Z; _lc = lc;
}
......@@ -223,6 +227,7 @@ struct Tet {
{
_modified=true;
V[0] = v0; V[1] = v1; V[2] = v2; V[3] = v3;
for (int i=0;i<4;i++)if (V[i])V[i]->setT(this);
// for (int i=0;i<4;i++)_copy[i] = *V[i];
return 1;
}
......@@ -232,6 +237,7 @@ struct Tet {
double val = robustPredicates::orient3d((double*)v0, (double*)v1,
(double*)v2, (double*)v3);
V[0] = v0; V[1] = v1; V[2] = v2; V[3] = v3;
for (int i=0;i<4;i++)if (V[i])V[i]->setT(this);
if (val > 0){
// for (int i=0;i<4;i++)_copy[i] = *V[i];
return 1;
......@@ -397,6 +403,7 @@ void delaunayTrgl(const unsigned int numThreads,
unsigned int Npts,
std::vector<Vert*> assignTo[],
tetContainer &allocator, edgeContainer *embedded = 0);
bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread);
void edgeSwapPass (int numThreads, tetContainer &allocator, edgeContainer &embeddedEdges);
void vertexRelocationPass (int numThreads, std::vector<Vert*> &v);
#endif
......@@ -142,20 +142,10 @@ void saturateEdge(Edge &e, std::vector<Vert*> &S, std::stack<IPT> &temp)
bgm = fields->get(fields->getBackgroundField());
}
if (0){
const double length = p1.distance(p2);
const double lc_avg = (e.first->lc() + e.second->lc())*.5;
const double l_adim_approx = length / lc_avg;
if (l_adim_approx > 10.){
SPoint3 pmid = (p1+p2)*0.5;
const double lc = GMSHSIZE(pmid, bgm, lc_avg);
S.push_back(new Vert(pmid.x(),pmid.y(),pmid.z(),lc));
return;
}
}
// double _a = Cpu();
// printf("%g %g \n",e.first->lc(), e.second->lc());
const double dN = adaptiveTrapezoidalRule(p1,p2,e.first->lc(), e.second->lc(),
_result, dl, temp, bgm);
// _C1 += Cpu() - _a;
......@@ -449,15 +439,25 @@ void computeMeshSizes (GRegion *gr, std::map<MVertex*, double> &s){
updateSize ((*it)->triangles[i]->getVertex((j+1)%3), s, d);
}
}
}
}
// there may be steiner points
for (unsigned int i=0;i<gr->tetrahedra.size();i++){
for (unsigned int j = 0; j < 4; j++){
MVertex *v = gr->tetrahedra[i]->getVertex(j);
if (s.find(v) == s.end()){
s[v] = 1.0;
}
}
}
}
extern double tetQuality (Tet* t, double *volume);
void edgeBasedRefinement(const int numThreads,
const int nptsatonce,
GRegion *gr)
{
// fill up old Datastructures
edgeContainer embeddedEdges (10000);
std::map<MVertex*, double> sizes;
......@@ -525,7 +525,6 @@ void edgeBasedRefinement(const int numThreads,
}
}
// do not allow to saturate boundary edges
{
for (unsigned int i=0;i< allocator.size(0);i++) {
......@@ -550,7 +549,6 @@ void edgeBasedRefinement(const int numThreads,
for (unsigned int i=0;i<_vertices.size();i++){
_filter.insert( _vertices[i] );
}
int iter = 1;
Msg::Info("------------------------------------- SATUR FILTR SORTH DELNY TIME TETS");
......@@ -561,10 +559,10 @@ void edgeBasedRefinement(const int numThreads,
double t1 = Cpu();
// C_COUNT = 0;
saturateEdges (ec, allocator, numThreads, add);
// printf("%d calls %d points added",C_COUNT,add.size());
// printf("%d points added",add.size());
double t2 = Cpu();
filterVertices (numThreads, _filter, add);
// printf(" %d points remain\n",add.size());
// printf(" %d points remain\n",add.size());
double t3 = Cpu();
if (add.empty())break;
// randomize vertices (EXTREMELY IMPORTANT FOR NOT DETERIORATING PERFORMANCE)
......@@ -572,6 +570,7 @@ void edgeBasedRefinement(const int numThreads,
// sort them using BRIO
std::vector<int> indices;
SortHilbert(add, indices);
// printf("%d indices\n",indices.size());
double t4 = Cpu();
delaunayTrgl (1,1,add.size(), &add,allocator,embeddedEdges.empty() ? NULL : &embeddedEdges);
double t5 = Cpu();
......@@ -588,12 +587,45 @@ void edgeBasedRefinement(const int numThreads,
}
}
std::vector<Vert*> vv;
for (int myThread=0; myThread < numThreads; myThread++)
for (unsigned int i=0;i<allocator.size(myThread);i++)
for (unsigned int j=0;j<4;j++)
if (allocator(myThread,i)->V[j])
allocator(myThread,i)->V[j]->setNum(0);
for (int myThread=0; myThread < numThreads; myThread++)
for (unsigned int i=0;i<allocator.size(myThread);i++)
for (unsigned int j=0;j<4;j++)
if (allocator(myThread,i)->V[j] && allocator(myThread,i)->V[j]->getNum() == 0){
allocator(myThread,i)->V[j]->setNum(1);
std::map<Vert*,MVertex*>::iterator it = _ma.find(allocator(myThread,i)->V[j]);
if (it == _ma.end())vv.push_back(allocator(myThread,i)->V[j]);
}
double t6 = Cpu();
Msg::Info("Optimizing");
edgeSwapPass (numThreads, allocator, embeddedEdges);
vertexRelocationPass (numThreads,vv);
edgeSwapPass (numThreads, allocator, embeddedEdges);
vertexRelocationPass (numThreads,vv);
edgeSwapPass (numThreads, allocator, embeddedEdges);
vertexRelocationPass (numThreads,vv);
double t7 = Cpu();
Msg::Info("Optimization done (%g seconds)",t7-t6);
int cat [10] = {0,0,0,0,0,0,0,0,0,0};
double MIN = 1.0;
for (unsigned int i=0; i< allocator.size(0);i++){
Tet *tt = allocator (0,i);
MVertex *mvs[4];
if (tt->V[0]){
double vol;
double q = tetQuality (tt,&vol);
cat [(int)(q*10)] ++;
MIN = std::min(MIN,q);
for (int j=0;j<4;j++){
Vert *v = tt->V[j];
std::map<Vert*,MVertex*>::iterator it = _ma.find(v);
......@@ -609,6 +641,11 @@ void edgeBasedRefinement(const int numThreads,
}
}
Msg::Info ("Min Tet Quality %22.15E",MIN);
for (int i=0;i<10;i++){
Msg::Info ("Tet Quality [%4.3f,%4.3f] %8d",0.1*i,0.1*(i+1),cat[i]);
}
if (Msg::GetVerbosity() == 99) {
std::map<Edge,double> _sizes;
for (unsigned int i=0; i< allocator.size(0);i++){
......@@ -638,4 +675,6 @@ void edgeBasedRefinement(const int numThreads,
Msg::Info("MESH EFFICIENCY : %22.15E %6d edges among %d are out of range",
exp (sum / _sizes.size()),nbBad,_sizes.size());
}
}
......@@ -107,7 +107,7 @@ bool tetgenmesh::reconstructmesh(void *p)
double t_start = Cpu();
std::vector<MVertex*> _vertices;
std::map<int,MVertex*> _extras;
// Get the set of vertices from GRegion.
{
std::set<MVertex*, MVertexLessThanNum> all;
......@@ -526,6 +526,7 @@ bool tetgenmesh::reconstructmesh(void *p)
clock_t t;
// printf("coucou2\n");
Msg::Info("Boundary Recovery...");
recoverboundary(t);
// printf("coucou3\n");
......@@ -651,7 +652,7 @@ bool tetgenmesh::reconstructmesh(void *p)
}
v->setIndex(pointmark(pointloop));
_gr->mesh_vertices.push_back(v);
_vertices.push_back(v);
_extras[pointmark(pointloop)] = v;
}
spivot(parentseg, parentsh);
if (parentsh.sh != NULL) {
......@@ -675,7 +676,7 @@ bool tetgenmesh::reconstructmesh(void *p)
}
v->setIndex(pointmark(pointloop));
_gr->mesh_vertices.push_back(v);
_vertices.push_back(v);
_extras[pointmark(pointloop)] = v;
}
}
// Record all the GFaces' tag at this segment.
......@@ -690,9 +691,8 @@ bool tetgenmesh::reconstructmesh(void *p)
// Create an interior mesh vertex.
MVertex *v = new MVertex(pointloop[0], pointloop[1], pointloop[2], _gr);
v->setIndex(pointmark(pointloop));
_gr->mesh_vertices.push_back(v);
_vertices.push_back(v);
}
_extras[pointmark(pointloop)] = v;
_gr->mesh_vertices.push_back(v); }
}
else if (pointtype(pointloop) == FREEFACETVERTEX) {
sdecode(point2sh(pointloop), parentsh);
......@@ -718,29 +718,34 @@ bool tetgenmesh::reconstructmesh(void *p)
}
v->setIndex(pointmark(pointloop));
_gr->mesh_vertices.push_back(v);
_vertices.push_back(v);
_extras[pointmark(pointloop)] = v;
}
else {
// Create a mesh vertex.
MVertex *v = new MVertex(pointloop[0], pointloop[1], pointloop[2], _gr);
v->setIndex(pointmark(pointloop));
_gr->mesh_vertices.push_back(v);
_vertices.push_back(v);
_extras[pointmark(pointloop)] = v;
}
}
else {
MVertex *v = new MVertex(pointloop[0], pointloop[1], pointloop[2], _gr);
v->setIndex(pointmark(pointloop));
_gr->mesh_vertices.push_back(v);
_vertices.push_back(v);
_extras[pointmark(pointloop)] = v;
}
}
pointloop = pointtraverse();
}
assert((int)_vertices.size() == points->items);
// assert((int)_vertices.size() == points->items);
}
if (!_extras.empty())Msg::Info("We add %d steiner points...",_extras.size());
if (l_edges.size() > 0) {
// There are Steiner points on segments!
face segloop;
// Re-create the segment mesh in the corresponding GEdges.
......@@ -756,7 +761,7 @@ bool tetgenmesh::reconstructmesh(void *p)
}
}
assert(ge != NULL);
Msg::Info("Steiner points exist on GEdge %d", ge->tag());
// Msg::Info("Steiner points exist on GEdge %d", ge->tag());
// Delete the old triangles.
for(unsigned int i = 0; i < ge->lines.size(); i++)
delete ge->lines[i];
......@@ -770,8 +775,11 @@ bool tetgenmesh::reconstructmesh(void *p)
if (shellmark(segloop) == etag) {
p[0] = sorg(segloop);
p[1] = sdest(segloop);
MVertex *v1 = _vertices[pointmark(p[0])];
MVertex *v2 = _vertices[pointmark(p[1])];
MVertex *v1 = pointmark(p[0]) >= _vertices.size() ? _extras[pointmark(p[0])] : _vertices[pointmark(p[0])];
MVertex *v2 = pointmark(p[1]) >= _vertices.size() ? _extras[pointmark(p[1])] : _vertices[pointmark(p[1])];
// MVertex *v2 = _vertices[pointmark(p[1])];
// printf("%d %d %d\n",pointmark(p[0]),pointmark(p[1]),_vertices.size());
// printf("%g %g %g %g\n",v1->x(),v1->y(),v2->x(),v2->y());
MLine *t = new MLine(v1, v2);
ge->lines.push_back(t);
}
......@@ -814,9 +822,15 @@ bool tetgenmesh::reconstructmesh(void *p)
p[0] = sorg(subloop);
p[1] = sdest(subloop);
p[2] = sapex(subloop);
MVertex *v1 = _vertices[pointmark(p[0])];
MVertex *v2 = _vertices[pointmark(p[1])];
MVertex *v3 = _vertices[pointmark(p[2])];
MVertex *v1 = pointmark(p[0]) >= _vertices.size() ? _extras[pointmark(p[0])] : _vertices[pointmark(p[0])];
MVertex *v2 = pointmark(p[1]) >= _vertices.size() ? _extras[pointmark(p[1])] : _vertices[pointmark(p[1])];
MVertex *v3 = pointmark(p[2]) >= _vertices.size() ? _extras[pointmark(p[2])] : _vertices[pointmark(p[2])];
if (pointmark(p[0]) >= _vertices.size())printf("F %d %d\n",v1->getIndex(),pointmark(p[0]));
if (pointmark(p[1]) >= _vertices.size())printf("F %d %d\n",v2->getIndex(),pointmark(p[1]));
if (pointmark(p[2]) >= _vertices.size())printf("F %d %d\n",v3->getIndex(),pointmark(p[2]));
// MVertex *v1 = _vertices[pointmark(p[0])];
// MVertex *v2 = _vertices[pointmark(p[1])];
// MVertex *v3 = _vertices[pointmark(p[2])];
MTriangle *t = new MTriangle(v1, v2, v3);
gf->triangles.push_back(t);
}
......@@ -837,10 +851,20 @@ bool tetgenmesh::reconstructmesh(void *p)
p[2] = apex(tetloop);
p[3] = oppo(tetloop);
MVertex *v1 = _vertices[pointmark(p[0])];
MVertex *v2 = _vertices[pointmark(p[1])];
MVertex *v3 = _vertices[pointmark(p[2])];
MVertex *v4 = _vertices[pointmark(p[3])];
MVertex *v1 = pointmark(p[0]) >= _vertices.size() ? _extras[pointmark(p[0])] : _vertices[pointmark(p[0])];
MVertex *v2 = pointmark(p[1]) >= _vertices.size() ? _extras[pointmark(p[1])] : _vertices[pointmark(p[1])];
MVertex *v3 = pointmark(p[2]) >= _vertices.size() ? _extras[pointmark(p[2])] : _vertices[pointmark(p[2])];
MVertex *v4 = pointmark(p[3]) >= _vertices.size() ? _extras[pointmark(p[3])] : _vertices[pointmark(p[3])];
// if (pointmark(p[0]) >= _vertices.size())printf("%d %d\n",v1->getIndex(),pointmark(p[0]));
// if (pointmark(p[1]) >= _vertices.size())printf("%d %d\n",v2->getIndex(),pointmark(p[1]));
// if (pointmark(p[2]) >= _vertices.size())printf("%d %d\n",v3->getIndex(),pointmark(p[2]));
// if (pointmark(p[3]) >= _vertices.size())printf("%d %d\n",v4->getIndex(),pointmark(p[3]));
// MVertex *v1 = _vertices[pointmark(p[0])];
// MVertex *v2 = _vertices[pointmark(p[1])];
// MVertex *v3 = _vertices[pointmark(p[2])];
// MVertex *v4 = _vertices[pointmark(p[3])];
MTetrahedron *t = new MTetrahedron(v1, v2, v3, v4);
_gr->tetrahedra.push_back(t);
tetloop.tet = tetrahedrontraverse();
......
......@@ -193,16 +193,18 @@ struct faceXtet{
}
};
void connectTets_vector2(std::vector<MTet4*> &t, std::vector<faceXtet> &conn)
template <class ITER>
void connectTets_vector2_templ(size_t _size, ITER beg, ITER end, std::vector<faceXtet> &conn)
{
conn.clear();
conn.reserve(4*t.size());
conn.reserve(4*_size);
// unsigned int k = 0;
// printf("COUCOU1 %d tets\n",t.size());
for (unsigned int i=0;i<t.size();i++){
if (!t[i]->isDeleted()){
for (ITER IT = beg; IT != end; ++IT){
MTet4* t = *IT;
if (!t->isDeleted()){
for (int j = 0; j < 4; j++){
conn.push_back(faceXtet(t[i], j));
conn.push_back(faceXtet(t, j));
}
}
}
......@@ -253,6 +255,8 @@ void connectTets(ITER beg, ITER end, std::set<MFace, Less_Face> *allEmbeddedFace
void connectTets(std::list<MTet4*> &l) { connectTets(l.begin(), l.end()); }
void connectTets(std::vector<MTet4*> &l) { connectTets(l.begin(), l.end()); }
void connectTets_vector2(std::list<MTet4*> &l, std::vector<faceXtet> &conn) { connectTets_vector2_templ(l.size(), l.begin(), l.end(), conn); }
void connectTets_vector2(std::vector<MTet4*> &l, std::vector<faceXtet> &conn) { connectTets_vector2_templ(l.size(), l.begin(), l.end(), conn); }
// Ensure the star-shapeness of the delaunay cavity
// We use the visibility criterion : the vertex should be visible
......@@ -865,7 +869,7 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
while (1){
// printf("coucou\n");
std::vector<MTet4*> newTets;
for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
/* for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
if (!(*it)->isDeleted()){
double qq = (*it)->getQuality();
if (qq < qMin){
......@@ -878,6 +882,7 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
}
}
}
*/
// printf("coucou\n");
illegals.clear();
......@@ -1188,6 +1193,11 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
}
gr->tetrahedra.clear();
{
std::vector<faceXtet> conn;
// connectTets_vector2_templ(allTets.size(), allTets.begin(), allTets.end() , conn);
}
// SLOW
connectTets(allTets.begin(), allTets.end());
// classify the tets on the right region
......
......@@ -128,10 +128,22 @@ void BuildSwapPattern3(SwapPattern *sc)
sc->trianguls = trgul ;
}
/*
0 1
+------------+
| |
| |
| |
+------------+
3 2
*/
void BuildSwapPattern4(SwapPattern *sc)
{
static int trgl[][3] =
{ {0,1,2}, {0,2,3}, {0,1,3}, {1,2,3} };
{ {0,1,2}, {2,3,0}, {1,2,3}, {3,0,1} };
static int trgul[][5] =
{ {0,1,-1,-1,-1}, {2,3,-1,-1,-1} };
......@@ -282,6 +294,7 @@ bool edgeSwap(std::vector<MTet4 *> &newTets,
// if there exist no swap that enhance the quality
if (best <= tetQualityRef) return false;
// does random swaps if (best < .01) return false;
// we have the best configuration, so we swap
// printf("outside size = %d\n",outside.size());
......
......@@ -60,6 +60,34 @@ static double objective_function (double xi, MVertex *ver, GFace *gf,
}
static double objective_function (double xi, MVertex *ver, GFace *gf,
SPoint3 &p1, SPoint3 &p2,
const std::vector<MElement*> &lt){
double x = ver->x();
double y = ver->y();
double z = ver->z();
SPoint3 p = p1 * (1.-xi) + p2 * xi;
double initialGuess[2]={0,0};
GPoint pp = gf->closestPoint(p,initialGuess);
ver->x() = pp.x();
ver->y() = pp.y();
ver->z() = pp.z();
double minQual = 1.0;
for(unsigned int i = 0; i < lt.size(); i++){
if (lt[i]->getNumVertices () == 4)
minQual = std::min((lt[i]->etaShapeMeasure()), minQual);
else
minQual = std::min(fabs(lt[i]->gammaShapeMeasure()), minQual);
}
ver->x() = x;
ver->y() = y;
ver->z() = z;
return minQual;
}
#define sqrt5 2.236067977499789696
......@@ -158,6 +186,55 @@ double Maximize_Quality_Golden_Section( MVertex *ver, GFace *gf,
return a;
}
double Maximize_Quality_Golden_Section( MVertex *ver, GFace *gf,
SPoint3 &p1, SPoint3 &p2,
const std::vector<MElement*> &lt ,
double tol, double &worst)
{
static const double lambda = 0.5 * (sqrt5 - 1.0);
static const double mu = 0.5 * (3.0 - sqrt5); // = 1 - lambda
double a = 0.0;
double b = 1.0;
worst = objective_function (0.0, ver, gf, p1, p2, lt );
if (worst > 0.5) return 0.0;
double x1 = b - lambda * (b - a);
double x2 = a + lambda * (b - a);
double fx1 = objective_function (x1, ver, gf, p1, p2, lt );
double fx2 = objective_function (x2, ver, gf, p1, p2, lt );
if (tol < 0.0)return fx1 > fx2 ? x1 : x2;
while ( ! Stopping_Rule( a, b , tol) ) {
// printf("GOLDEN : %g %g (%12.5E,%12.5E)\n",a,b,fa,fb);
if (fx1 < fx2) {
a = x1;
if ( Stopping_Rule( a, b , tol) ) break;
x1 = x2;
fx1 = fx2;
x2 = b - mu * (b - a);
fx2 = objective_function (x2, ver, gf, p1, p2, lt );
} else {
b = x2;
if ( Stopping_Rule( a, b , tol) ) break;
x2 = x1;
fx2 = fx1;
x1 = a + mu * (b - a);
fx1 = objective_function (x1, ver, gf, p1, p2, lt );
}
}
double final = objective_function (a, ver, gf, p1, p2, lt );
if (final < worst) return 0.0;
worst = final;
// printf("finally : %g %g (%12.5E,%12.5E)\n",a,b,fa,fb);
return a;
}
void _relocateVertexGolden(MVertex *ver,
const std::vector<MElement*> &lt, double relax , double tol)
{
......@@ -184,6 +261,37 @@ void _relocateVertexGolden(MVertex *ver,
ver->z() = (1.-xi) * ver->z() + xi * z/N;
}
// use real space + projection at the end
static double _relocateVertex2(GFace* gf, MVertex *ver,
const std::vector<MElement*> &lt,
double tol) {
SPoint3 p1(0,0,0);
int counter = 0;
for(unsigned int i = 0; i < lt.size(); i++){
for (int j=0;j<lt[i]->getNumVertices();j++){
MVertex* v = lt[i]->getVertex(j);
p1 += SPoint3(v->x(),v->y(),v->z());
counter++;
}
}
p1 *= 1./(double)counter;
SPoint3 p2(ver->x(),ver->y(),ver->z());
double worst;
double xi = Maximize_Quality_Golden_Section( ver, gf, p1, p2, lt , tol, worst);
SPoint3 p = p1*(1-xi) + p2*xi;
double initialGuess[2]={0,0};
GPoint pp = gf->closestPoint(p,initialGuess);
if (!pp.succeeded())return 2.0;
ver->x() = pp.x();
ver->y() = pp.y();
ver->z() = pp.z();
return worst;
}
static double _relocateVertex(GFace* gf, MVertex *ver,
const std::vector<MElement*> &lt,
double tol) {
......@@ -191,15 +299,19 @@ static double _relocateVertex(GFace* gf, MVertex *ver,
SPoint2 p1(0,0);
SPoint2 p2;
ver->getParameter(0,p2[0]);
ver->getParameter(1,p2[1]);
if (ver->getParameter(0,p2[0])){
ver->getParameter(1,p2[1]);
}
else {
return _relocateVertex2(gf,ver,lt,tol);
}
int counter=0;
for(unsigned int i = 0; i < lt.size(); i++){
for (int j=0;j<lt[i]->getNumVertices();j++){
MVertex* v = lt[i]->getVertex(j);
SPoint2 pp;
reparamMeshVertexOnFace(v, gf, pp);
reparamMeshVertexOnFace(v, gf, pp);
counter++;
if (v->onWhat()->dim() == 1) {
GEdge *ge = dynamic_cast<GEdge*> (v->onWhat());
......@@ -230,7 +342,7 @@ void getAllBoundaryLayerVertices (GFace *gf, std::set<MVertex*> &vs);
void RelocateVertices (GFace* gf, int niter, double tol) {
std::set<MVertex*> vs;
getAllBoundaryLayerVertices (gf, vs);
v2t_cont adj;
buildVertexToElement(gf->triangles, adj);
buildVertexToElement(gf->quadrangles, adj);
......
......@@ -691,11 +691,11 @@ void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed
// add the vertices as additional vertices in the
// surface mesh
char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag());
FILE *f = Fopen(ccc,"w");
if(f) fprintf(f, "View \"\"{\n");
// char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag());
// FILE *f = Fopen(ccc,"w");
// if(f) fprintf(f, "View \"\"{\n");
for (unsigned int i=0;i<vertices.size();i++){
if(f) vertices[i]->print(f,i);
// if(f) vertices[i]->print(f,i);
if(vertices[i]->_v->onWhat() == gf) {
packed.push_back(vertices[i]->_v);
metrics.push_back(vertices[i]->_meshMetric);
......@@ -704,10 +704,10 @@ void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed
}
delete vertices[i];
}
if(f){
fprintf(f,"};");
fclose(f);
}
// if(f){
// fprintf(f,"};");
// fclose(f);
// }
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment