diff --git a/Geo/GFace.h b/Geo/GFace.h
index f1ca7c23220da3a6ceb2d1266d82b1f7a514947e..111a53adf09d32b9d5a1754dadd6ba0e4df55fcc 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -89,6 +89,8 @@ public:
// add embedded vertices/edges
void addEmbeddedVertex(GVertex *v) { embedded_vertices.insert(v); }
void addEmbeddedEdge(GEdge *e) { embedded_edges.push_back(e); }
+ // FIXME DEBUG
+ void clearEmbeddedEdges() { embedded_edges.clear(); }
// delete the edge from the face (the edge is supposed to be a free
// edge in the face, not part of any edge loops--use with caution!)
@@ -327,7 +329,6 @@ public:
std::vector<SPoint2> stl_vertices_uv;
std::vector<SPoint3> stl_vertices_xyz;
std::vector<SVector3> stl_normals;
- std::vector<SVector3> stl_curvatures;
std::vector<int> stl_triangles;
// a vertex array containing a geometrical representation of the
@@ -362,8 +363,7 @@ public:
std::vector<SVector3> storage3; // sizes and directions storage
std::vector<double> storage4; // sizes and directions storage
- virtual bool reorder(const int elementType,
- const std::vector<std::size_t> &ordering);
+ virtual bool reorder(const int elementType, const std::vector<std::size_t> &ordering);
};
#endif
diff --git a/Geo/GModelParametrize.cpp b/Geo/GModelParametrize.cpp
deleted file mode 100644
index bbd18693916c155c76da270586591391f18b8123..0000000000000000000000000000000000000000
--- a/Geo/GModelParametrize.cpp
+++ /dev/null
@@ -1,393 +0,0 @@
-// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
-
-#include <vector>
-#include <set>
-#include <map>
-#include <stack>
-#include "GmshConfig.h"
-#include "GModel.h"
-#include "GFace.h"
-#include "discreteFace.h"
-#include "discreteEdge.h"
-#include "MTriangle.h"
-#include "MEdge.h"
-#include "GEdge.h"
-#include "MLine.h"
-#include "Context.h"
-#include "GmshMessage.h"
-#include "GModelParametrize.h"
-
-#if defined(HAVE_MESH)
-#include "meshPartition.h"
-#endif
-
-#if defined(HAVE_HXT)
-extern "C" {
-#include "hxt_mesh.h"
-#include "hxt_edge.h"
-#include "hxt_mean_values.h"
-}
-
-static HXTStatus gmsh2hxt(GFace *gf, HXTMesh **pm,
- std::map<MVertex *, int> &v2c,
- std::vector<MVertex *> &c2v)
-{
- HXTContext *context;
- hxtContextCreate(&context);
- HXTMesh *m;
- HXT_CHECK(hxtMeshCreate(context, &m));
- std::set<MVertex *> all;
- for(size_t i = 0; i < gf->triangles.size(); i++) {
- all.insert(gf->triangles[i]->getVertex(0));
- all.insert(gf->triangles[i]->getVertex(1));
- all.insert(gf->triangles[i]->getVertex(2));
- }
- m->vertices.num = m->vertices.size = all.size();
- HXT_CHECK(
- hxtAlignedMalloc(&m->vertices.coord, 4 * m->vertices.num * sizeof(double)));
-
- size_t count = 0;
- c2v.resize(all.size());
- for(std::set<MVertex *>::iterator it = all.begin(); it != all.end(); it++) {
- m->vertices.coord[4 * count + 0] = (*it)->x();
- m->vertices.coord[4 * count + 1] = (*it)->y();
- m->vertices.coord[4 * count + 2] = (*it)->z();
- m->vertices.coord[4 * count + 3] = 0.0;
- v2c[*it] = count;
- c2v[count++] = *it;
- }
- all.clear();
-
- m->triangles.num = m->triangles.size = gf->triangles.size();
- HXT_CHECK(hxtAlignedMalloc(&m->triangles.node,
- (m->triangles.num) * 3 * sizeof(uint32_t)));
- HXT_CHECK(hxtAlignedMalloc(&m->triangles.colors,
- (m->triangles.num) * sizeof(uint16_t)));
- for(size_t i = 0; i < gf->triangles.size(); i++) {
- m->triangles.node[3 * i + 0] = v2c[gf->triangles[i]->getVertex(0)];
- m->triangles.node[3 * i + 1] = v2c[gf->triangles[i]->getVertex(1)];
- m->triangles.node[3 * i + 2] = v2c[gf->triangles[i]->getVertex(2)];
- m->triangles.colors[i] = gf->tag();
- }
-
- m->lines.num = m->lines.size = 0;
-
- *pm = m;
- return HXT_STATUS_OK;
-}
-
-HXTStatus computeDiscreteCurvatures(GModel *gm, std::map<MVertex* , std::pair<SVector3, SVector3> > &C ){
- C.clear();
-
- for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
- HXTMesh *m;
- HXTEdges *edges;
- double *nodalCurvatures;
- double *crossField;
- std::map<MVertex *, int> v2c;
- std::vector<MVertex *> c2v;
- gmsh2hxt(*it, &m, v2c, c2v);
- HXT_CHECK(hxtEdgesCreate(m, &edges));
- HXT_CHECK(hxtCurvatureRusinkiewicz(m, &nodalCurvatures, &crossField, edges, false));
- for (size_t i=0;i<m->vertices.num;i++){
- MVertex *v = c2v[i];
- double *c = &nodalCurvatures[6*i];
- SVector3 cMax (c[0],c[1],c[2]);
- SVector3 cMin (c[3],c[4],c[5]);
- std::pair<SVector3, SVector3> out = std::make_pair(cMax,cMin);
- C.insert(std::make_pair(v,out));
- }
- HXT_CHECK(hxtEdgesDelete(&edges));
- HXT_CHECK(hxtFree(&nodalCurvatures));
- HXT_CHECK(hxtFree(&crossField));
- HXT_CHECK(hxtMeshDelete(&m));
- }
- return HXT_STATUS_OK;
-}
-
-
-void parametrizeAllGEdge(GModel *gm)
-{
- for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it) {
- discreteEdge *de = dynamic_cast<discreteEdge *>(*it);
- if (de) de->createGeometry();
- }
-}
-
-int parametrizeAllGFace(GModel *gm, std::map<MVertex* , std::pair<SVector3, SVector3> > *C)
-{
- int result = 0;
- for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
- discreteFace *df = dynamic_cast<discreteFace *>(*it);
- if(df) result += parametrizeGFace(df,C);
- }
- return result;
-}
-
-int parametrizeGFace(discreteFace *gf, std::map<MVertex* , std::pair<SVector3, SVector3> > *C)
-{
- int n = 1;
- if(gf->triangles.empty()) return 0;
- HXT_CHECK(hxtInitializeLinearSystems(&n, NULL));
- HXTMesh *m;
- HXTMeanValues *param;
- HXTEdges *edges;
- std::map<MVertex *, int> v2c;
- std::vector<MVertex *> c2v;
- gmsh2hxt(gf, &m, v2c, c2v);
- // double *nodalCurvatures;
- // double *crossField;
- // printf("A %d triangles\n",m->triangles.num);
- HXT_CHECK(hxtEdgesCreate(m, &edges));
- HXT_CHECK(hxtMeanValuesCreate(edges, ¶m));
- HXT_CHECK(hxtMeanValuesCompute(param));
- double *uvc = NULL;
- int nv, ne;
- HXT_CHECK(hxtMeanValuesGetData(param, NULL, NULL, &uvc, &nv, &ne,1));
- gf->stl_vertices_uv.clear();
- gf->stl_vertices_uv.resize(nv);
- gf->stl_vertices_xyz.clear();
- gf->stl_vertices_xyz.resize(nv);
- gf->stl_curvatures.clear();
- if (C) gf->stl_curvatures.resize(2*nv);
- gf->stl_normals.clear();
- gf->stl_normals.resize(nv);
-
- for(int iv = 0; iv < nv; iv++) {
- if (C){
- MVertex *v = c2v[iv];
- std::map<MVertex* , std::pair<SVector3, SVector3> >::iterator it = C->find(v);
- if (it == C->end())Msg::Error("POINT %d NOT FOUND FOR COMPUTING CURVATURES",v->getNum());
- // printf("%g %g %g\n",it->second.first.x(),it->second.first.y(),it->second.first.z());
- gf->stl_curvatures[2*iv] = it->second.first;
- gf->stl_curvatures[2*iv+1] = it->second.second;
- }
- gf->stl_vertices_uv[iv] = SPoint2(uvc[2*iv],uvc[2*iv+1]);
-
- gf->stl_vertices_xyz[iv] =
- SPoint3(m->vertices.coord[4 * iv + 0], m->vertices.coord[4 * iv + 1],
- m->vertices.coord[4 * iv + 2]);
- gf->stl_normals[iv] = SVector3(1, 0, 0);
- }
- gf->stl_triangles.clear();
- gf->stl_triangles.resize(3 * ne);
- for(int ie = 0; ie < ne; ie++) {
- gf->stl_triangles[3 * ie + 0] = m->triangles.node[3 * ie + 0];
- gf->stl_triangles[3 * ie + 1] = m->triangles.node[3 * ie + 1];
- gf->stl_triangles[3 * ie + 2] = m->triangles.node[3 * ie + 2];
- }
- // printf("A\n");
- gf->fillVertexArray(false);
- HXT_CHECK(hxtMeshDelete(&m));
- HXT_CHECK(hxtEdgesDelete(&edges));
- HXT_CHECK(hxtFree(&uvc));
- // printf("B\n");
- return 0;
-}
-#else
-int parametrizeGFace(GFace *gf, std::map<std::pair<MVertex*,GFace*> , std::pair<SVector3, SVector3> > *C)
-{
- Msg::Error("Gmsh should be compiled against HXT for being able to compute "
- "discrete parametrizations");
- return -1;
-}
-#endif
-
-int isTriangulationParametrizable(const std::vector<MTriangle *> &t, int Nmax,
- double ar)
-{
- int XX = (int)t.size();
- if(XX > Nmax) return XX / Nmax + 1;
- std::set<MVertex *> v;
- std::map<MEdge, int, Less_Edge> e;
- double surf = 0;
- for(size_t i = 0; i < t.size(); ++i) {
- surf += t[i]->getVolume();
- for(int j = 0; j < 3; j++) {
- v.insert(t[i]->getVertex(j));
- std::map<MEdge, int, Less_Edge>::iterator it = e.find(t[i]->getEdge(j));
- if(it == e.end())
- e[t[i]->getEdge(j)] = 1;
- else
- it->second++;
- }
- }
- std::map<MEdge, int, Less_Edge>::iterator it = e.begin();
- std::vector<MEdge> _bnd;
- for(; it != e.end(); ++it){
- if(it->second == 1) _bnd.push_back(it->first);
- }
-
- std::vector<std::vector<MVertex *> > vs;
- if(!SortEdgeConsecutive(_bnd, vs)) {
- // we have vertices adjacent to more than 2 model edges
- // FILE *f = fopen("bug.pos","w");
- // fprintf(f,"View\"\"{\n");
- // for (size_t i=0;i<t.size();i++){
- // fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
- // t[i]->getVertex(0)->x(),t[i]->getVertex(0)->y(),t[i]->getVertex(0)->z(),
- // t[i]->getVertex(1)->x(),t[i]->getVertex(1)->y(),t[i]->getVertex(1)->z(),
- // t[i]->getVertex(2)->x(),t[i]->getVertex(2)->y(),t[i]->getVertex(2)->z(),
- // t[i]->getVertex(0)->getNum(),t[i]->getVertex(1)->getNum(),t[i]->getVertex(2)->getNum());
- // }
- // fprintf(f,"};\n");
- // fclose(f);
- // getchar();
- return 2;
- }
- double lmax = 0;
- for(size_t i = 0; i < vs.size(); i++) {
- double li = 0;
- for(size_t j = 1; j < vs[i].size(); j++) {
- li += distance(vs[i][j], vs[i][j - 1]);
- }
- lmax = std::max(li, lmax);
- }
-
- if(ar * lmax * lmax < 2 * M_PI * surf) { return 4; }
-
- int poincare =
- t.size() - (2 * (v.size() - 1) - _bnd.size() + 2 * (vs.size() - 1));
- return poincare <= 0 ? 1 : 2;
-}
-
-void makeMLinesUnique(std::vector<MLine *> &v)
-{
- std::sort(v.begin(), v.end(), compareMLinePtr());
- v.erase(std::unique(v.begin(), v.end(), equalMLinePtr()), v.end());
-}
-
-void computeEdgeCut(GModel *gm, std::vector<MLine *> &cut,
- int max_elems_per_cut)
-{
- GModel m;
-
- // ----------------------------------------------------------------------------------
- // STUPID FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
- std::vector<MTriangle*> aa;
- for(size_t i = 0; i < (*it)->triangles.size(); i++){
- if ((*it)->triangles[i]->getVertex(0) == (*it)->triangles[i]->getVertex(1) ||
- (*it)->triangles[i]->getVertex(0) == (*it)->triangles[i]->getVertex(2) ||
- (*it)->triangles[i]->getVertex(1) == (*it)->triangles[i]->getVertex(2)){
- // printf("TRIANGLE BAD %d %d %d\n",(*it)->triangles[i]->getVertex(0)->getNum()
- // ,(*it)->triangles[i]->getVertex(1)->getNum()
- // ,(*it)->triangles[i]->getVertex(2)->getNum());
- }
- else aa.push_back((*it)->triangles[i]);
- }
- (*it)->triangles = aa;
- }
- // END STUPID FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- // ----------------------------------------------------------------------------------
-
- for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
- int part = 0;
- if((*it)->triangles.empty()) continue;
- std::vector<MVertex *> verts = (*it)->mesh_vertices;
- std::map<MTriangle *, int> global;
- std::map<MEdge, int, Less_Edge> cuts;
- std::stack<std::vector<MTriangle *> > partitions;
- partitions.push((*it)->triangles);
- (*it)->triangles.clear();
-
- while(!partitions.empty()) {
- (*it)->triangles = partitions.top();
- (*it)->mesh_vertices.clear();
- std::set<MVertex *> vs;
- for(size_t i = 0; i < (*it)->triangles.size(); ++i) {
- for(size_t j = 0; j < 3; ++j)
- vs.insert((*it)->triangles[i]->getVertex(j));
- }
- (*it)->mesh_vertices.insert((*it)->mesh_vertices.begin(), vs.begin(),
- vs.end());
- partitions.pop();
- int np =
- isTriangulationParametrizable((*it)->triangles, max_elems_per_cut, 5.0);
- if(np == 1) {
- for(size_t i = 0; i < (*it)->triangles.size(); i++)
- global[(*it)->triangles[i]] = part;
- part++;
- }
- else {
-#if defined(HAVE_MESH)
- int *p = new int[(*it)->triangles.size()];
- if(!PartitionFace(*it, np, p)) {
- std::vector<MTriangle *> t[1000];
- for(size_t i = 0; i < (*it)->triangles.size(); i++)
- t[(*it)->triangles[i]->getPartition()].push_back(
- (*it)->triangles[i]);
- for(int i = 0; i < np; i++) { partitions.push(t[i]); }
- }
- delete[] p;
-#else
- Msg::Error("Partitioning surface requires Mesh module");
-#endif
- }
- }
- (*it)->triangles.clear();
- std::map<MTriangle *, int>::iterator it2 = global.begin();
- for(; it2 != global.end(); ++it2) {
- MTriangle *t = it2->first;
- (*it)->triangles.push_back(t);
- for(int i = 0; i < 3; i++) {
- MEdge ed = t->getEdge(i);
- std::map<MEdge, int, Less_Edge>::iterator it3 = cuts.find(ed);
- if(it3 == cuts.end())
- cuts[ed] = it2->second;
- else {
- if(it3->second != it2->second)
- cut.push_back(new MLine(ed.getVertex(0), ed.getVertex(1)));
- }
- }
- }
- (*it)->mesh_vertices = verts;
- }
- makeMLinesUnique(cut);
-}
-
-void computeNonManifoldEdges(GModel *gm, std::vector<MLine *> &cut,
- bool addBoundary)
-{
- std::map<MEdge, int, Less_Edge> m;
- for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
- for(size_t i = 0; i < (*it)->triangles.size(); i++) {
- for(int j = 0; j < 3; j++) {
- std::map<MEdge, int, Less_Edge>::iterator it2 =
- m.find((*it)->triangles[i]->getEdge(j));
- if(it2 == m.end())
- m[(*it)->triangles[i]->getEdge(j)] = 1;
- else
- it2->second++;
- }
- }
- }
- {
- int countNM = 0, countBND = 0;
- std::map<MEdge, int, Less_Edge>::iterator it = m.begin();
- for(; it != m.end(); ++it) {
- if(it->second > 2) {
- cut.push_back(
- new MLine(it->first.getVertex(0), it->first.getVertex(1)));
- countNM++;
- }
- if(addBoundary && it->second == 1) {
- cut.push_back(
- new MLine(it->first.getVertex(0), it->first.getVertex(1)));
- countBND++;
- }
- }
- if(countNM + countBND > 0)
- Msg::Info(
- "Model has %d non manifold mesh edges and %d boundary mesh edges",
- countNM, countBND);
- }
- makeMLinesUnique(cut);
-}
-
-// todo...
-void detectPlanarFaces (GModel *gm, std::vector<MLine *> &cut){
-}
diff --git a/Geo/GModelParametrize.h b/Geo/GModelParametrize.h
deleted file mode 100644
index 22b93922d80fb8043222ff7b60f7046cf226a7c2..0000000000000000000000000000000000000000
--- a/Geo/GModelParametrize.h
+++ /dev/null
@@ -1,21 +0,0 @@
-// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
-
-#ifndef GMODEL_PARAMETRIZE_H
-#define GMODEL_PARAMETRIZE_H
-
-int isTriangulationParametrizable(const std::vector<MTriangle *> &t, int Nmax,
- double ar);
-int isTriangulationParametrizable(const std::vector<MTriangle *> &t);
-void computeEdgeCut(GModel *gm, std::vector<MLine *> &cut,
- int max_elems_per_cut);
-void computeNonManifoldEdges(GModel *gm, std::vector<MLine *> &cut,
- bool addBoundary);
-int parametrizeAllGFace(GModel *gm, std::map<MVertex* , std::pair<SVector3, SVector3> > *C=NULL);
-void parametrizeAllGEdge(GModel *gm);
-int parametrizeGFace(discreteFace *gf, std::map<MVertex* , std::pair<SVector3, SVector3> > *C=NULL);
-HXTStatus computeDiscreteCurvatures(GModel *gm, std::map<MVertex* , std::pair<SVector3, SVector3> > &C );
-
-#endif
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index f8ee62b17876f0bed476e70b430859c2e6859b76..8c5ac7de56fdade3de3acd379f70601049dbab58 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1362,12 +1362,12 @@ static inline bool getOrderedNeighboringVertices (BDS_Point *p,
if (p->iD == CHECK){
for (size_t i = 0 ; i < ts.size() ; i++){
BDS_Point *pts[4];
- ts[i]->getNodes(pts);
- //printf("TR %d : %d %d %d\n",i,pts[0]->iD,pts[1]->iD,pts[2]->iD);
+ ts[i]->getNodes(pts);
+ printf("TR %d : %d %d %d\n",i,pts[0]->iD,pts[1]->iD,pts[2]->iD);
}
}
- if (ts.empty())return false;
+ if (ts.empty())return false;
while (1){
bool found = false;
for (size_t i = 0 ; i < ts.size() ; i++){
@@ -1398,11 +1398,11 @@ static inline bool getOrderedNeighboringVertices (BDS_Point *p,
}
}
}
-
+
if (nbg.size() == ts.size()) break;
if (!found) return false;
}
-
+
if (p->iD == CHECK){
printf("FINALLY : ");
for (size_t i=0;i< nbg.size() ;i++){
@@ -1447,7 +1447,7 @@ static inline void getCentroidUV (const BDS_Point *p, GFace *gf,
factSum += fact;
U += kernel[i].x()*fact;
V += kernel[i].y()*fact;
- LC += lc[i]*fact;
+ LC += lc[i]*fact;
}
U /= factSum;
V /= factSum;
@@ -1488,16 +1488,16 @@ static inline void computeSomeKindOfKernel (const BDS_Point *p,
std::vector<double> &lc,
int check){
-
+
FILE *f = NULL;
if (p->iD == check) f = fopen("kernel.pos","w");
-
+
SPoint2 pp (p->u,p->v);
if (p->iD == check){
fprintf(f,"View \"kernel\"{\n");
fprintf(f,"SP(%g,%g,0){2};\n",p->u,p->v);
}
-
+
double ll = p->lc();
kernel.clear();
lc.clear();
@@ -1512,7 +1512,7 @@ static inline void computeSomeKindOfKernel (const BDS_Point *p,
}
else if (nbg[(i+1)%nbg.size()]->degenerated){
kernel.push_back(SPoint2 (nbg[i]->u,nbg[i]->v));
- kernel.push_back(SPoint2 (nbg[i]->u,nbg[(i+1)%nbg.size()]->v));
+ kernel.push_back(SPoint2 (nbg[i]->u,nbg[(i+1)%nbg.size()]->v));
lc.push_back(nbg[i]->lc());
lc.push_back(nbg[i]->lc());
}
@@ -1552,7 +1552,7 @@ static inline void computeSomeKindOfKernel (const BDS_Point *p,
lc[i] = lc_now;
}
-
+
if (p->iD == check){
// for (size_t i = 0; i<kernel.size();i++){
// fprintf(f,"SL(%g,%g,0,%g,%g,0){3,3};\n",kernel[i].x(),kernel[i].y(),kernel[(i+1)%nbg.size()].x(),kernel[(i+1)%nbg.size()].y());
@@ -1570,8 +1570,8 @@ static inline void getGradientTutteEnergy (const BDS_Point *p,
double dEduv[2]){
// X = X0 + dX/dU (U-U0)
// E = \sum_i (X-X_i)^2 = \sum_i (U-U_i)^T G (U-U_i)
- // dE/dU = G \sum_i (U-U_i)
-
+ // dE/dU = G \sum_i (U-U_i)
+
Pair<SVector3, SVector3> d = gf->firstDer(SPoint2 (p->u,p->v));
double G11 = dot (d.left(),d.left());
double G12 = dot (d.left(),d.right());
@@ -1655,7 +1655,7 @@ static inline bool minimizeTutteEnergyProj (BDS_Point *p,
if (p->iD == check){
printf("%g %g %d\n",p->u,p->v,validityOfCavity (p, nbg));
}
-
+
if (validityOfCavity (p, nbg)){
p->X = gp.x(); p->Y = gp.y(); p->Z = gp.z();
double E_moved = getTutteEnergy (p, nbg, RATIO) ;
@@ -1663,7 +1663,7 @@ static inline bool minimizeTutteEnergyProj (BDS_Point *p,
return true;
}
}
-
+
p->X = oldX; p->Y = oldY; p->Z = oldZ; p->u = oldU; p->v = oldV;
return false;
}
@@ -1683,9 +1683,9 @@ static inline bool minimizeTutteEnergyParam (BDS_Point *p,
double E_moved = getTutteEnergy (p, nbg, RATIO2);
if (p->iD == check)printf("%g vs %g\n",E_unmoved ,E_moved);
-
+
if (E_moved < E_unmoved ) {
- p->u = U; p->v = V;
+ p->u = U; p->v = V;
if (!validityOfCavity (p, nbg)){
p->X = oldX; p->Y = oldY; p->Z = oldZ; p->u = oldU; p->v = oldV;
return false;
@@ -1708,7 +1708,7 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, double threshold)
}
int CHECK = -1;
-
+
std::vector<BDS_Point *> nbg;
std::vector<double> lc;
std::vector<SPoint2> kernel;
@@ -1721,14 +1721,14 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, double threshold)
if (RATIO > threshold)return false;
computeSomeKindOfKernel (p,nbg,kernel,lc, CHECK);
-
- if (! minimizeTutteEnergyParam ( p,E_unmoved, RATIO,nbg,kernel,lc,gf, CHECK)){
- if ( ! minimizeTutteEnergyProj (p,E_unmoved, RATIO,nbg,kernel,lc,gf, CHECK)) {
+
+ if (! minimizeTutteEnergyParam ( p,E_unmoved, RATIO,nbg,kernel,lc,gf, CHECK)){
+ if ( ! minimizeTutteEnergyProj (p,E_unmoved, RATIO,nbg,kernel,lc,gf, CHECK)) {
// __COUNT3++;
return false;
}
else{
- p->config_modified = true;
+ p->config_modified = true;
E_unmoved = getTutteEnergy (p, nbg, RATIO);
minimizeTutteEnergyProj (p,E_unmoved, RATIO,nbg,kernel,lc,gf, CHECK);
// __COUNT2++;
@@ -1740,5 +1740,7 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, double threshold)
}
return true;
-
+
}
+
+
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 4c322e01f597bc7f5398080ec2adca2959690bdb..2cc3a1e8c83b36b6b22f184d71087374f7e00867 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -82,7 +82,7 @@ static inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace
static double computeEdgeLinearLength(BDS_Edge *e, GFace *f)
{
// FIXME !!!
- return
+ return
f->geomType() == GEntity::Plane ?
e->length() :
computeEdgeLinearLength(e->p1, e->p2, f);
@@ -283,7 +283,7 @@ static void getDegeneracy(BDS_Mesh &m, std::vector<BDS_Point*> °)
while(itp != m.points.end()) {
if ((*itp)->degenerated)deg.push_back(*itp);
itp++;
- }
+ }
}
static void setDegeneracy(std::vector<BDS_Point*> °, bool d){
@@ -292,12 +292,12 @@ static void setDegeneracy(std::vector<BDS_Point*> °, bool d){
static int validitiyOfGeodesics(GFace *gf, BDS_Mesh &m)
{
- std::vector< BDS_Edge *> degenerated;
+ std::vector< BDS_Edge *> degenerated;
for (size_t i=0;i< m.edges.size();i++){
BDS_Edge *e = m.edges[i];
if (!e->deleted && e->p1->degenerated + e->p2->degenerated == 1)degenerated.push_back(e);
}
- std::vector< BDS_Edge *> toSplit;
+ std::vector< BDS_Edge *> toSplit;
for (size_t i=0;i< degenerated.size();i++){
BDS_Edge *ed = degenerated[i];
double u3[2] = {ed->p1->degenerated ? ed->p2->u : ed->p1->u,ed->p1->v};
@@ -312,7 +312,7 @@ static int validitiyOfGeodesics(GFace *gf, BDS_Mesh &m)
double rp1 = robustPredicates::orient2d(u1, u2, u3);
double rp2 = robustPredicates::orient2d(u1, u2, u4);
double rq1 = robustPredicates::orient2d(u3, u4, u1);
- double rq2 = robustPredicates::orient2d(u3, u4, u2);
+ double rq2 = robustPredicates::orient2d(u3, u4, u2);
if(rp1*rp2<0 && rq1*rq2<0){
// printf("%d %d -- %d %d\n",ed->p1->iD,ed->p2->iD,e->p1->iD,e->p2->iD);
toSplit.push_back(ed);
@@ -481,7 +481,7 @@ void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP,
count++;
double lone1 = 0.;
bool collapseP1Allowed = true;
-
+
double lone2 = 0.;
bool collapseP2Allowed = false;
if(e->p2->iD > MAXNP) {
@@ -526,7 +526,7 @@ void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q, double thr
// __COUNT2=0;
// __COUNT3=0;
double t1 = Cpu();
-
+
for (int i=0; i<1;i++){
std::set<BDS_Point *, PointLessThan>::iterator itp = m.points.begin();
while(itp != m.points.end()) {
@@ -547,6 +547,48 @@ void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q, double thr
// printf("%d %d %d smooth \n",__COUNT1, __COUNT2, __COUNT3);
}
+static void CHECK_STRANGE(const char *c, BDS_Mesh &m)
+{
+ return;
+#if 0
+ int strange = 0;
+ for (size_t i=0;i<m.triangles.size();++i){
+ BDS_Point *pts[4];
+ m.triangles[i]->getNodes(pts);
+ double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3);
+ if (pts[0]->degenerated+pts[1]->degenerated+pts[2]->degenerated <= 1){
+ if (fabs(surface_triangle_param(pts[0],pts[1],pts[2]))<1.e-8){
+ strange ++;
+ }
+ }
+ }
+ if (strange) printf("strange(%s) = %d\n",c,strange);
+
+ return;
+ for (size_t int i=0;i<m.triangles.size();++i){
+ BDS_Point *pts[4];
+ m.triangles[i]->getNodes(pts);
+ if (pts[0] == pts[1] ||
+ pts[0] == pts[2] ||
+ pts[1] == pts[2]){
+ strange ++;
+ }
+ }
+
+ return;
+ strange = 0;
+ std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
+ while (itp != m.points.end()){
+ if ((*itp)->g && (*itp)->g->classif_degree == 2){
+ std::vector<BDS_Face*> t = (*itp)->getTriangles();
+ if (t.size() == 2)strange++;//printf("vertex %d \n",(*itp)->iD);
+ }
+ ++itp;
+ }
+ printf("strange(%s) = %d\n",c,strange);
+#endif
+}
+
static void computeNodalSizes(
GFace *gf, BDS_Mesh &m,
std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap)
@@ -652,7 +694,7 @@ void modifyInitialMeshToRemoveDegeneracies(
it->first->degenerated = true;
}
}
- for(size_t i = 0; i < degenerated_edges.size(); i++) {
+ for(size_t i = 0; i < degenerated_edges.size(); i++) {
m.collapse_edge_parametric(degenerated_edges[i], degenerated_edges[i]->p1,
true);
recoverMap->erase(degenerated_edges[i]->p1);
@@ -717,9 +759,9 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
setDegeneracy(deg,false);
}
else setDegeneracy(deg,true);
-
+
// printf("degeneracy %d\n",nbSplit);
-
+
// we count the number of local mesh modifs.
int nb_split = 0;
int nb_smooth = 0;
@@ -742,7 +784,7 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
// outputScalarField(m.triangles, "splitsmooth.pos", 1, gf);
-
+
swapEdgePass(gf, m, nb_swap,t_sw);
smoothVertexPass(gf, m, nb_smooth, false,.5,t_sm);
@@ -766,6 +808,8 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
collapseEdgePass(gf, m, .45, MAXNP, nb_collaps,t_col);
smoothVertexPass(gf, m, nb_smooth, false,.5,t_sm);
}
+
+ // CHECK_STRANGE("smmooth", m);
m.cleanup();
@@ -785,7 +829,7 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
}
if((minL > MINE_ && maxL < MAXE_) || IT > (abs(NIT))) break;
-
+
IT++;
Msg::Debug(" iter %3d minL %8.3f/%8.3f maxL %8.3f/%8.3f -- Small/Large/Total (%6d/%6d/%6d): "
"%6d splits, %6d swaps, %6d collapses, %6d moves",
@@ -867,7 +911,7 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
}
setDegeneracy(deg,true);
-
+
double t_total = t_spl + t_sw + t_col + t_sm;
if(!t_total) t_total = 1.e-6;
Msg::Debug(" ---------------------------------------");
diff --git a/contrib/hxt/hxt_mean_values.c b/contrib/hxt/hxt_mean_values.c
index adb2c4761018432884624790c8a097c581380c7b..375b7b226b18bff8faae9491de39bb04c8da828c 100644
--- a/contrib/hxt/hxt_mean_values.c
+++ b/contrib/hxt/hxt_mean_values.c
@@ -34,10 +34,7 @@ HXTStatus hxtMeanValuesCreate(HXTEdges *edges, HXTMeanValues **meanValues)
// be moved into parametrization?
HXTBoundaries *boundaries;
- printf("1 %d %d %d\n",mesh->triangles.node[0],mesh->triangles.node[1],mesh->triangles.node[2]);
- printf("1 %d %d %d\n",mesh->triangles.node[3],mesh->triangles.node[4],mesh->triangles.node[5]);
HXT_CHECK(hxtEdgesSetBoundaries(edges, &boundaries));
- printf("2\n");
map->boundaries = boundaries;
@@ -55,7 +52,7 @@ HXTStatus hxtMeanValuesCreate(HXTEdges *edges, HXTMeanValues **meanValues)
}
}
- HXT_CHECK(hxtMalloc(&map->hole,(nBoundaries+1)*sizeof(int)));
+ HXT_CHECK(hxtMalloc(&map->hole,(nBoundaries-1)*sizeof(int)));
int c = 0;
for(int i=0; i<nBoundaries; i++)
if(i!=map->boundary){
@@ -85,12 +82,11 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
HXT_CHECK(hxtBoundariesGetNumberOfEdgesOfLineLoop(meanValues->boundaries,meanValues->hole[i],&c));
nby += c;
}
-
-#ifdef __FILLERGAP_
+
uint32_t *fillergap;
HXT_CHECK(hxtMalloc(&fillergap, 3*(nTriangles+nby)*sizeof(uint32_t)));
memcpy(fillergap, mesh->triangles.node, 3*nTriangles*sizeof(uint32_t));
-
+
int s=0;
for(int i=0; i<meanValues->nHoles; i++){
uint32_t *cedgs;
@@ -104,27 +100,20 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
}
s += n_;
}
-#endif // __FILLERGAP_
-
/*
for(int i=0; i<nTriangles+nby; i++)
printf("tri %d (nby=%d) \t %u %u %u\n",i,nby,fillergap[3*i+0],fillergap[3*i+1],fillergap[3*i+2]);
*/
-
+
HXTLinearSystem *sys;
-
-#ifdef __FILLERGAP_
+
HXT_CHECK(hxtLinearSystemCreateLU(&sys,nTriangles+nby,3,1,fillergap));
-#else
- HXT_CHECK(hxtLinearSystemCreateLU(&sys,nTriangles,3,1,mesh->triangles.node));
-#endif
-
-
+
uint32_t *flag;
HXT_CHECK(hxtMalloc(&flag,nNodes*sizeof(uint32_t)));
for(uint32_t ii=0; ii<nNodes; ii++)
flag[ii] = 0;
-
+
/* boundary conditions */
double *uv = meanValues->uv;
int n;
@@ -141,28 +130,28 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
uv[2*edges->node[2*edges_ll[i]+0]+1] = sin(angle);
currentLength += hxtEdgesLength(meanValues->initialEdges, edges_ll[i]);
}
-
+
double *U, *V, *Urhs, *Vrhs;
//printf("allocation: %d\n",nNodes+meanValues->nHoles);
HXT_CHECK( hxtMalloc(&U,(nNodes+meanValues->nHoles)*sizeof(double)) );
HXT_CHECK( hxtMalloc(&V,(nNodes+meanValues->nHoles)*sizeof(double)) );
HXT_CHECK( hxtMalloc(&Urhs,(nNodes+meanValues->nHoles)*sizeof(double)) );
HXT_CHECK( hxtMalloc(&Vrhs,(nNodes+meanValues->nHoles)*sizeof(double)) );
-
-
+
+
// init linear system
HXT_CHECK(hxtLinearSystemZeroMatrix(sys));
for(uint32_t i=0; i<(nNodes+meanValues->nHoles); i++){
Urhs[i] = 0.;
Vrhs[i] = 0.;
}
-
+
// setting linear system
for(uint32_t ie=0; ie<edges->numEdges; ie++){
int ik[2] = {-1,-1};
uint64_t *tri = edges->edg2tri + 2*ie;
-
+
// gather local edge index
for(int it=0; it<2; it++){
if(tri[it]==(uint64_t)-1)
@@ -174,21 +163,21 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
}
}
}
-
+
for(int ij=0; ij<2; ij++){
-
+
uint32_t v0 = edges->node[2*ie+ij];
uint32_t v1 = edges->node[2*ie+(1-ij)];
-
+
if (flag[v0]==1){//boundary nodes/conditons
- HXT_CHECK(hxtLinearSystemSetMatrixRowIdentity(sys,v0,0));
+ HXT_CHECK(hxtLinearSystemSetMatrixRowIdentity(sys,v0,0));
HXT_CHECK(hxtLinearSystemSetRhsEntry(sys,Urhs, v0,0, uv[2*v0+0]));
HXT_CHECK(hxtLinearSystemSetRhsEntry(sys,Vrhs, v0,0, uv[2*v0+1]));
}
else {// inner node
- double e[3] = {mesh->vertices.coord[4*v1+0]-mesh->vertices.coord[4*v0+0],mesh->vertices.coord[4*v1+1]-mesh->vertices.coord[4*v0+1],mesh->vertices.coord[4*v1+2]-mesh->vertices.coord[4*v0+2]};
- double ne = sqrt(e[0]*e[0]+e[1]*e[1]+e[2]*e[2]);
-
+ double e[3] = {mesh->vertices.coord[4*v1+0]-mesh->vertices.coord[4*v0+0],mesh->vertices.coord[4*v1+1]-mesh->vertices.coord[4*v0+1],mesh->vertices.coord[4*v1+2]-mesh->vertices.coord[4*v0+2]};
+ double ne = sqrt(e[0]*e[0]+e[1]*e[1]+e[2]*e[2]);
+
uint32_t vLeft = mesh->triangles.node[3*tri[0] + (ik[0]+2)%3];
double a[3] = {mesh->vertices.coord[4*vLeft+0]-mesh->vertices.coord[4*v0+0],mesh->vertices.coord[4*vLeft+1]-mesh->vertices.coord[4*v0+1],mesh->vertices.coord[4*vLeft+2]-mesh->vertices.coord[4*v0+2]};
double na = sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
@@ -210,9 +199,9 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
}
}
}// end for int ie
-
-#ifdef __FILLERGAP_
+
+
// holes are filled somehow
for(int i=0; i<meanValues->nHoles; i++){
@@ -269,8 +258,6 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
HXT_CHECK(hxtFree(&c));
HXT_CHECK(hxtFree(&d));
}
-#endif
-
HXT_CHECK(hxtLinearSystemSolve(sys,Urhs,U));
HXT_CHECK(hxtLinearSystemSolve(sys,Vrhs,V));
for(uint32_t i=0; i<nNodes; i++){
@@ -278,9 +265,7 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
uv[2*i+1] = V[i] ;
}
-#ifdef __FILLERGAP_
HXT_CHECK(hxtFree(&fillergap));
-#endif
HXT_CHECK(hxtFree(&flag));
HXT_CHECK(hxtFree(&U));
HXT_CHECK(hxtFree(&V));
@@ -320,7 +305,7 @@ HXTStatus hxtMeanValueAspectRatio(HXTMeanValues *meanValues, int *aspectRatio)
-HXTStatus hxtMeanValuesGetData(HXTMeanValues *mv, uint64_t **global,uint32_t **gn, double **uv, int *nv, int *ne, int onlyuv){
+HXTStatus hxtMeanValuesGetData(HXTMeanValues *mv, uint64_t **global,uint32_t **gn, double **uv, int *nv, int *ne){
*nv = mv->initialEdges->edg2mesh->vertices.num;
*ne = mv->initialEdges->edg2mesh->triangles.num;
@@ -333,7 +318,7 @@ HXTStatus hxtMeanValuesGetData(HXTMeanValues *mv, uint64_t **global,uint32_t **g
uv_[iv] = mv->uv[iv];
*uv=uv_;
}
- if (onlyuv) return HXT_STATUS_OK;
+ // return HXT_STATUS_OK;
uint64_t *global_;
HXT_CHECK(hxtMalloc(&global_,(*ne)*sizeof(uint64_t)));
diff --git a/contrib/hxt/hxt_mean_values.h b/contrib/hxt/hxt_mean_values.h
index fcbf246ffc23ba277128d423d68c264abba0abb1..c2415b000a3cd69ab7e1b333baf0cc6e42933ea2 100644
--- a/contrib/hxt/hxt_mean_values.h
+++ b/contrib/hxt/hxt_mean_values.h
@@ -11,7 +11,7 @@ HXTStatus hxtMeanValuesCreate(HXTEdges *edges, HXTMeanValues **meanValues);
HXTStatus hxtMeanValuesDelete(HXTMeanValues **meanValues);
HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues);
HXTStatus hxtMeanValueAspectRatio(HXTMeanValues *meanValues, int *aspectRatio);
-HXTStatus hxtMeanValuesGetData(HXTMeanValues *mv, uint64_t **global, uint32_t **gn, double **uv, int *nv, int *ne, int onlyuv);
+HXTStatus hxtMeanValuesGetData(HXTMeanValues *mv, uint64_t **global, uint32_t **gn, double **uv, int *nv, int *ne);
HXTStatus hxtMeanValuesWrite(HXTMeanValues *meanValues, const char* filename);
HXTStatus hxtMeanValuesWriteParamMesh(HXTMeanValues *meanValues, const char* filename);
diff --git a/contrib/hxt/hxt_parametrization.c b/contrib/hxt/hxt_parametrization.c
index 36d4a37920ee608ed6525561d1dc3613db9977f4..9f1ec807f3344c2a7f517041d158215f2bcf8219 100644
--- a/contrib/hxt/hxt_parametrization.c
+++ b/contrib/hxt/hxt_parametrization.c
@@ -10,7 +10,7 @@
int* colors; //gives colors of elements: colors[i] = color of i-th element
int* nNodes; //gives lower index of int* nodes for a color: nNodes[i] = first node index of i-th color
-int* nodes; //gives the global index of initial 'mesh': nodes[i] = if i=nNodes[c]+k, gives the global index of k-th node of c-th color
+int* nodes; //gives the global index of initial 'mesh': nodes[i] = if i=nNodes[c]+k, gives the global index of k-th node of c-th color
double* uv; //give the uv param coordinates of a local node: uv[2*i+p] = u/v component of the global k-node of the c-th color
*/
@@ -29,9 +29,9 @@ static HXTStatus hxtVectorInit(HXTVector **vecptr){
HXT_CHECK(hxtMalloc(vecptr,sizeof(HXTVector)));
-
+
HXTVector *new = *vecptr;
-
+
new->ptr = NULL;
new->last = -1;
new->length = 0;
@@ -43,7 +43,7 @@ static HXTStatus hxtVectorInit(HXTVector **vecptr){
static HXTStatus hxtVectorFree(HXTVector **vecptr){
HXTVector *vector = *vecptr;
-
+
HXT_CHECK(hxtFree(&vector->ptr));
HXT_CHECK(hxtFree(vecptr));
@@ -60,7 +60,7 @@ static HXTStatus hxtVectorPushBack(HXTVector *vector, void *object){
if(i+1 >= length){
if(length==0)
length=1;
- HXT_CHECK(hxtRealloc(&vector->ptr,2*(length)*sizeof(void *)));
+ HXT_CHECK(hxtRealloc(&vector->ptr,2*(length)*sizeof(void *)));
vector->length = 2*length;
}
vector->ptr[i+1] = object;
@@ -87,10 +87,10 @@ static HXTStatus hxtDoubleLongestEdge(HXTEdges *e, int tri, double *l)
{
uint32_t *ed = &e->tri2edg[3*tri];
*l=0;
-
+
for(int i=0; i<3; i++){
uint64_t *tr= &e->edg2tri[2*ed[i]];
-
+
if(tr[1]==(uint64_t)-1)
continue;
@@ -107,7 +107,7 @@ static HXTStatus hxtLongestEdge(HXTEdges *edges, int tri,int *ie,int comp)
*ie=comp;
for(int i=0; i<3; i++){
uint64_t *tr= &edges->edg2tri[2*ed[i]];
-
+
if(tr[1]==(uint64_t)-1)
continue;
@@ -136,7 +136,7 @@ static HXTStatus hxtLongestEdgeBisection(HXTEdges *edges,int nrefinements)
HXTMesh *m = edges->edg2mesh;
-
+
uint64_t *flag=NULL;
HXT_CHECK(hxtMalloc(&flag,m->triangles.num*sizeof(uint64_t)));
for(uint64_t it=0; it<m->triangles.num; it++){
@@ -167,7 +167,7 @@ static HXTStatus hxtLongestEdgeBisection(HXTEdges *edges,int nrefinements)
break;
}
}
- do{
+ do{
HXT_CHECK(hxtLongestEdge(e,ti,&ei,ej));
tj = e->edg2tri[2*ei+0]==ti ? e->edg2tri[2*ei+1] : e->edg2tri[2*ei+0];
HXT_CHECK(hxtLongestEdge(e,tj,&ej,ei));
@@ -272,10 +272,10 @@ static HXTStatus hxtLongestEdgeBisection(HXTEdges *edges,int nrefinements)
if(t_[0]==tj)
t_[0] = ttj;
else
- t_[1] = ttj;
+ t_[1] = ttj;
}
- }//end for
+ }//end for
e->numEdges +=3;
if(e->numEdges > maxEdg){
maxEdg *= 2;
@@ -285,7 +285,7 @@ static HXTStatus hxtLongestEdgeBisection(HXTEdges *edges,int nrefinements)
}
e->node[2*ei+1] = m->vertices.num-1;
e->edg2tri[2*ei+0] = ti;
- e->edg2tri[2*ei+1] = ttj;
+ e->edg2tri[2*ei+1] = ttj;
e->node[2*e0+0] = m->vertices.num-1;
e->node[2*e0+1] = nj;
e->edg2tri[2*e0+0] = tti;
@@ -321,8 +321,8 @@ static HXTStatus hxtLongestEdgeBisection(HXTEdges *edges,int nrefinements)
}//end refinement
-
-
+
+
return HXT_STATUS_OK;
}
@@ -331,29 +331,29 @@ static HXTStatus hxtPartitioning(HXTEdges *edges,int *part, int nPartitions,HXTV
{
HXTMesh *mesh = edges->edg2mesh;
-
+
uint32_t *flagVertices;
HXT_CHECK(hxtMalloc(&flagVertices,mesh->vertices.num*sizeof(uint32_t)));
for(int p = 0; p<nPartitions; p++){
HXTMesh *msh=NULL;
HXT_CHECK(hxtMeshCreate(mesh->ctx,&msh));
-
+
for(uint32_t iv=0; iv<mesh->vertices.num; iv++)
flagVertices[iv] = -1;
-
+
uint32_t cv = 0;
uint64_t ce = 0;
- for(uint64_t ie=0; ie<mesh->triangles.num; ie++){
+ for(uint64_t ie=0; ie<mesh->triangles.num; ie++){
if(part[ie]==p){
for(int jv=0; jv<3; jv++){
uint32_t current = mesh->triangles.node[3*ie+jv];
if(flagVertices[current] == (uint32_t)-1)
- flagVertices[current] = cv++;
+ flagVertices[current] = cv++;
}
ce++;
}
- }
+ }
/*
printf("------------stat \t part %d---------------\n",p);
printf("\t number of vertices:%u\n",cv);
@@ -368,13 +368,13 @@ static HXTStatus hxtPartitioning(HXTEdges *edges,int *part, int nPartitions,HXTV
msh->vertices.coord[4*current+0] = mesh->vertices.coord[4*iv+0];
msh->vertices.coord[4*current+1] = mesh->vertices.coord[4*iv+1];
msh->vertices.coord[4*current+2] = mesh->vertices.coord[4*iv+2];
- msh->vertices.coord[4*current+3] = mesh->vertices.coord[4*iv+3];
+ msh->vertices.coord[4*current+3] = mesh->vertices.coord[4*iv+3];
}
- }
+ }
msh->triangles.num = ce;
//<
- uint64_t *global;
+ uint64_t *global;
HXT_CHECK(hxtMalloc(&global,ce*sizeof(uint64_t)));
//>
HXT_CHECK( hxtAlignedMalloc(&msh->triangles.node,ce*3*sizeof(uint32_t)) );
@@ -418,8 +418,8 @@ static HXTStatus hxtCheckConnectivity(HXTEdges *e,HXTVector *partitions)
uint64_t *idx = NULL;
HXT_CHECK(hxtMalloc(&idx,(m->triangles.num+1)*sizeof(uint64_t)));
-
-
+
+
int count = 0;
flag[0] = count;
queue[0] = 0;
@@ -440,8 +440,8 @@ static HXTStatus hxtCheckConnectivity(HXTEdges *e,HXTVector *partitions)
last++;
idx[iq+1]++;
flag[next] = (int) count;
- }
- }
+ }
+ }
}
count++;
goOn = 0;
@@ -451,7 +451,7 @@ static HXTStatus hxtCheckConnectivity(HXTEdges *e,HXTVector *partitions)
iq=last;
queue[last] = i;
last++;
- flag[i] = count;
+ flag[i] = count;
break;
}
}
@@ -461,8 +461,8 @@ static HXTStatus hxtCheckConnectivity(HXTEdges *e,HXTVector *partitions)
for(uint64_t i=0; i<m->triangles.num; i++){
uint32_t *refVert = &m->triangles.node[3*queue[i]];
for(uint64_t j=idx[i]+flag[queue[i]]; j<idx[i+1]+flag[queue[i]]; j++){
-
- int cond = 0;
+
+ int cond = 0;
uint32_t *checkVert = &m->triangles.node[3*queue[j]];
for(int ii=0; ii<3; ii++){
for(int jj=0; jj<3; jj++){
@@ -474,7 +474,7 @@ static HXTStatus hxtCheckConnectivity(HXTEdges *e,HXTVector *partitions)
if(cond==1)
break;
}//end for ii
- if(cond==0){
+ if(cond==0){
printf("Unconsistent orientation \t %lu (%u %u %u) ; (%u %u)-(%u %u)-(%u %u) [%d] ~ %lu (%u %u %u) ; (%u %u)-(%u %u)-(%u %u) [%d].\n",queue[i],refVert[0],refVert[1],refVert[2],e->node[2*e->tri2edg[3*queue[i]+0]+0],e->node[2*e->tri2edg[3*queue[i]+0]+1],e->node[2*e->tri2edg[3*queue[i]+1]+0],e->node[2*e->tri2edg[3*queue[i]+1]+1],e->node[2*e->tri2edg[3*queue[i]+2]+0],e->node[2*e->tri2edg[3*queue[i]+2]+1],flag[queue[i]],queue[j],checkVert[0],checkVert[1],checkVert[2],e->node[2*e->tri2edg[3*queue[j]+0]+0],e->node[2*e->tri2edg[3*queue[j]+0]+1],e->node[2*e->tri2edg[3*queue[j]+1]+0],e->node[2*e->tri2edg[3*queue[j]+1]+1],e->node[2*e->tri2edg[3*queue[j]+2]+0],e->node[2*e->tri2edg[3*queue[j]+2]+1],flag[queue[j]]);
uint64_t temp = m->triangles.node[3*queue[j]+0];
m->triangles.node[3*queue[j]+0] = m->triangles.node[3*queue[j]+1];
@@ -484,9 +484,9 @@ static HXTStatus hxtCheckConnectivity(HXTEdges *e,HXTVector *partitions)
e->tri2edg[3*queue[j]+2] = te;
printf("Consistent (?) orientation \t %lu (%u %u %u) ; (%u %u)-(%u %u)-(%u %u) [%d] ~ %lu (%u %u %u) ; (%u %u)-(%u %u)-(%u %u) [%d].\n",queue[i],refVert[0],refVert[1],refVert[2],e->node[2*e->tri2edg[3*queue[i]+0]+0],e->node[2*e->tri2edg[3*queue[i]+0]+1],e->node[2*e->tri2edg[3*queue[i]+1]+0],e->node[2*e->tri2edg[3*queue[i]+1]+1],e->node[2*e->tri2edg[3*queue[i]+2]+0],e->node[2*e->tri2edg[3*queue[i]+2]+1],flag[queue[i]],queue[j],checkVert[0],checkVert[1],checkVert[2],e->node[2*e->tri2edg[3*queue[j]+0]+0],e->node[2*e->tri2edg[3*queue[j]+0]+1],e->node[2*e->tri2edg[3*queue[j]+1]+0],e->node[2*e->tri2edg[3*queue[j]+1]+1],e->node[2*e->tri2edg[3*queue[j]+2]+0],e->node[2*e->tri2edg[3*queue[j]+2]+1],flag[queue[j]]);
-
+
}//end if cond==0
- }//end for j
+ }//end for j
}//end for i
if (partitions!=NULL)
@@ -500,10 +500,10 @@ static HXTStatus hxtSplitEdges(HXTEdges *edges,int nPartitions,HXTVector *partit
{
HXTMesh *mesh = edges->edg2mesh;
-
+
uint32_t *tri2edg = edges->tri2edg;
- uint64_t *edg2tri = edges->edg2tri;
-
+ uint64_t *edg2tri = edges->edg2tri;
+
HXTBoundaries *boundaries;
HXT_CHECK(hxtEdgesSetBoundaries(edges, &boundaries));
int nbedg;
@@ -511,42 +511,42 @@ static HXTStatus hxtSplitEdges(HXTEdges *edges,int nPartitions,HXTVector *partit
int nVertex = (int) mesh->triangles.num;
int nEdge = (int) (edges->numEdges) - nbedg;
-
+
int *idx=NULL;
HXT_CHECK(hxtMalloc(&idx,(nVertex+1)*sizeof(int)));
int *nbh=NULL;
HXT_CHECK(hxtMalloc(&nbh, 2*nEdge*sizeof(int)));
int *weights=NULL;
HXT_CHECK(hxtMalloc(&weights,2*nEdge*sizeof(int)));
-
+
idx[0] = 0;
for(int i=0; i<nVertex; i++){// triangle by triangle
uint32_t *current = &tri2edg[3*i];
int temp = 0;
- for(int j=0; j<3; j++){
+ for(int j=0; j<3; j++){
uint64_t *local = &edg2tri[2*current[j]];
if(local[1]==(uint64_t)-1)
continue;
- else{
+ else{
nbh[idx[i]+temp] = i == (int) local[0] ? (int) local[1] : (int) local[0];
weights[idx[i]+temp] = (int) (hxtEdgesLength(edges,current[j])+1);
temp++;
}
}
-
+
idx[i+1] = idx[i] + temp;
}
-
- // only for k-way
+
+ // only for k-way
idx_t options[METIS_NOPTIONS];
METIS_SetDefaultOptions(options);
options[METIS_OPTION_NCUTS] = nPartitions;
//options[METIS_OPTION_MINCONN] = 1;
//options[METIS_OPTION_CCORDER] = 1;
- options[METIS_OPTION_CONTIG] = 1;
+ options[METIS_OPTION_CONTIG] = 1;
//options[METIS_OPTION_DBGLVL] = 1;
-
+
int edgeCut;
int *part;
HXT_CHECK(hxtMalloc(&part,nVertex*sizeof(int)));
@@ -559,7 +559,7 @@ static HXTStatus hxtSplitEdges(HXTEdges *edges,int nPartitions,HXTVector *partit
HXT_CHECK(hxtFree(&idx));
HXT_CHECK(hxtPartitioning(edges,part,nPartitions,partitions));
-
+
return HXT_STATUS_OK;
}
@@ -571,9 +571,9 @@ static HXTStatus hxtCuttingProcess(HXTEdges *edges, int ar, int bool, HXTVector
HXT_CHECK(hxtBoundariesGetNumberOfLineLoops(boundaries,&nll));
int g = (edges->numEdges - edges->edg2mesh->vertices.num - edges->edg2mesh->triangles.num + 2 - nll)/2;
HXT_CHECK(hxtBoundariesGetSeamPoint(boundaries,&seamPoint));
-
+
HXTVector *toSplit;
- HXT_CHECK(hxtVectorInit(&toSplit));
+ HXT_CHECK(hxtVectorInit(&toSplit));
if(seamPoint!=0 || ar==0 || g!=0 || nll == 0)
HXT_CHECK(hxtVectorPushBack(toSplit,edges));
else
@@ -583,7 +583,7 @@ static HXTStatus hxtCuttingProcess(HXTEdges *edges, int ar, int bool, HXTVector
HXT_CHECK(hxtVectorInit(&split));
HXTEdges *beingSplitted = (HXTEdges*) toSplit->ptr[i];
HXT_CHECK(hxtSplitEdges(beingSplitted,2,split));
- if(bool>0){
+ if(bool>0){
HXT_CHECK(hxtMeshDelete(&beingSplitted->edg2mesh));
HXT_CHECK(hxtEdgesDelete(&beingSplitted));
}
@@ -606,7 +606,7 @@ static HXTStatus hxtCuttingProcess(HXTEdges *edges, int ar, int bool, HXTVector
}
HXT_CHECK(hxtVectorFree(&toSplit));
// free boundaries
-
+
return HXT_STATUS_OK;
}
@@ -624,7 +624,7 @@ HXTStatus hxtParametrizationCreate(HXTMesh *mesh, int nrefinements, HXTParametri
// printf("checking connectivity ...\n");
HXT_CHECK(hxtCheckConnectivity(param0->edges,NULL));
// HXT_CHECK(hxtLongestEdgeBisection(param0->edges,nrefinements));
-
+
HXTEdges *edges=NULL;
HXT_CHECK(hxtEdgesCreate(initialMesh, &edges));
uint64_t *global;
@@ -634,43 +634,43 @@ HXTStatus hxtParametrizationCreate(HXTMesh *mesh, int nrefinements, HXTParametri
edges->global = global;
param0->edges->global = global;
- HXTVector *connect;
- HXT_CHECK(hxtVectorInit(&connect));
+ HXTVector *connect;
+ HXT_CHECK(hxtVectorInit(&connect));
HXT_CHECK(hxtCheckConnectivity(edges,connect));
// printf("connectivity is done\n");
-
+
HXTVector *toparam;
HXT_CHECK(hxtVectorInit(&toparam));
for(int i=0; i<connect->last+1; i++){
HXT_CHECK(hxtCuttingProcess(connect->ptr[i],1,0,toparam));
}
- //printf("initial cutting process is done, %d\n", param0->n);
+ printf("initial cutting process is done, %d\n", param0->n);
// return HXT_STATUS_OK;
-
-
+
+
HXTVector *atlas;
HXT_CHECK(hxtVectorInit(&atlas));
for(int i=0; i<toparam->last+1; i++){
HXTMeanValues *param;
HXTEdges *current = (HXTEdges *)(toparam->ptr[i]);
HXT_CHECK(hxtMeanValuesCreate(current,¶m));
- HXT_CHECK(hxtMeanValuesCompute(param));
+ HXT_CHECK(hxtMeanValuesCompute(param));
int ar = 0;
HXT_CHECK(hxtMeanValueAspectRatio(param,&ar));
-
+
if (0 && ar==0)
HXT_CHECK(hxtCuttingProcess(current,ar,toparam->last,toparam));
else
HXT_CHECK(hxtVectorPushBack(atlas,param));
}
hxtVectorFree(&toparam);
-
+
param0->n = atlas->last+1;
HXT_CHECK(hxtMalloc(¶m0->maps,param0->n*sizeof(HXTMeanValues*)));
- for(int i=0; i<param0->n; i++)
+ for(int i=0; i<param0->n; i++)
param0->maps[i] = (HXTMeanValues *) atlas->ptr[i];
-
+
HXT_CHECK(hxtVectorFree(&atlas));
return HXT_STATUS_OK;
}
@@ -688,7 +688,7 @@ HXTStatus hxtParametrizationCompute(HXTParametrization *parametrization, int **_
int *colors=NULL, *nNodes=NULL, *nodes=NULL;
double *uv=NULL;
-
+
*m = parametrization->edges->edg2mesh;
int Ncolors = parametrization->n;
*nc=Ncolors;
@@ -699,33 +699,33 @@ HXTStatus hxtParametrizationCompute(HXTParametrization *parametrization, int **_
for(int c=0; c<Ncolors; c++){
uint64_t *global=NULL;
- int nv, ne;
- HXT_CHECK(hxtMeanValuesGetData(parametrization->maps[c],&global, NULL, NULL, &nv,&ne,0));
+ int nv, ne;
+ HXT_CHECK(hxtMeanValuesGetData(parametrization->maps[c],&global, NULL, NULL, &nv,&ne));
nNodes[c+1] = nNodes[c]+nv;
for(int ie=0; ie<ne; ie++)
colors[global[ie]] = c;
-
+
HXT_CHECK(hxtFree(&global));
-
+
}
-
+
HXT_CHECK(hxtMalloc(&uv,2*nNodes[Ncolors]*sizeof(double)));//
HXT_CHECK(hxtMalloc(&nodes,nNodes[Ncolors]*sizeof(int)));
for(int it=0; it<nNodes[Ncolors]; it++)
nodes[it] = -1;
-
+
for(int c=0; c<Ncolors; c++){
uint64_t *global=NULL;
uint32_t *gn=NULL;
double *uvc=NULL;
- int nv, ne;
- HXT_CHECK(hxtMeanValuesGetData(parametrization->maps[c],&global, &gn, &uvc, &nv,&ne,0));
+ int nv, ne;
+ HXT_CHECK(hxtMeanValuesGetData(parametrization->maps[c],&global, &gn, &uvc, &nv,&ne));
for(int iv=0; iv<2*nv; iv++)
uv[2*nNodes[c]+iv] = uvc[iv];
-
+
for(int ie=0; ie<ne; ie++){
for(int kk=0; kk<3; kk++){
int gvn = (int) parametrization->edges->edg2mesh->triangles.node[3*global[ie]+kk];
@@ -738,7 +738,7 @@ HXTStatus hxtParametrizationCompute(HXTParametrization *parametrization, int **_
HXT_CHECK(hxtFree(&uvc));
}
-
+
*_colors_=colors;
*_nNodes_=nNodes;
*_nodes_=nodes;
@@ -748,7 +748,7 @@ HXTStatus hxtParametrizationCompute(HXTParametrization *parametrization, int **_
}
HXTStatus hxtParametrizationWrite(HXTParametrization *parametrization, const char *filename)
-{
+{
return HXT_STATUS_OK;
}