diff --git a/Geo/fourierModel.cpp b/Geo/fourierModel.cpp
index 0a4c958a2d306161c450aca862c35667731a4790..420a363abb13be6e1867ad31004f54daf6902b79 100644
--- a/Geo/fourierModel.cpp
+++ b/Geo/fourierModel.cpp
@@ -3,6 +3,23 @@
#include "Context.h"
#include "Views.h"
+/*
+ Investigate the following approach:
+
+ compute and store the pou of each overlapping patch in the nodes of
+ all the patches
+
+ for each pair of overlapping patches, find the line pou1=pou2 by
+ interpolation on the overlapping grids
+
+ compute the intersections of these lines
+
+ this should define a non-overlapping partitioning of the grid, which
+ could be used as the boundary constrain for the unstructured algo
+
+ Would that work??
+*/
+
extern Context_T CTX;
#if defined(HAVE_FOURIER_MODEL)
@@ -12,8 +29,55 @@ extern Context_T CTX;
static model *FM = 0;
+void debugVertices(std::vector<MVertex*> &vertices, std::string file,
+ bool parametric, int num=0)
+{
+ char name[256];
+ sprintf(name, "%s_%d.pos", file.c_str(), num);
+ FILE *fp = fopen(name, "w");
+ fprintf(fp, "View \"debug\"{\n");
+ for(unsigned int i = 0; i < vertices.size(); i++){
+ double x, y, z;
+ if(parametric){
+ vertices[i]->getParameter(0, x);
+ vertices[i]->getParameter(1, y);
+ z = 0;
+ }
+ else{
+ x = vertices[i]->x();
+ y = vertices[i]->y();
+ z = vertices[i]->z();
+ }
+ fprintf(fp, "SP(%g,%g,%g){%d};\n", x, y, z, i);
+ }
+ fprintf(fp, "};\n");
+ fclose(fp);
+}
+
+void debugElements(std::vector<MElement*> &elements, std::string file,
+ bool parametric, int num=0)
+{
+ char name[256];
+ sprintf(name, "%s_%d.pos", file.c_str(), num);
+ FILE *fp = fopen(name, "w");
+ fprintf(fp, "View \"debug\"{\n");
+ for(unsigned int i = 0; i < elements.size(); i++)
+ elements[i]->writePOS(fp);
+ fprintf(fp, "};\n");
+ fclose(fp);
+}
+
class meshCartesian{
public:
+ typedef MDataFaceVertex<double> DVertex;
+ static std::set<DVertex*, MVertexLessThanLexicographic> vPosition;
+ static std::set<MVertex*> vDelete;
+ static void deleteUnusedVertices()
+ {
+ for(std::set<MVertex*>::iterator it = vDelete.begin(); it != vDelete.end(); it++)
+ delete *it;
+ vDelete.clear();
+ }
void operator() (GFace *gf)
{
int M = (int)(30. / CTX.mesh.lc_factor), N = (int)(30. / CTX.mesh.lc_factor);
@@ -24,8 +88,17 @@ public:
GPoint p = gf->point(u, v);
double pou;
FM->GetPou(gf->tag(), u, v, pou);
- gf->mesh_vertices.push_back
- (new MDataFaceVertex<double>(p.x(), p.y(), p.z(), gf, u, v, pou));
+ DVertex *w = new DVertex(p.x(), p.y(), p.z(), gf, u, v, pou);
+ // eliminate duplicate vertices on hard edges
+ std::set<DVertex*, MVertexLessThanLexicographic>::iterator it = vPosition.find(w);
+ if(it != vPosition.end()){
+ delete w;
+ gf->mesh_vertices.push_back(*it);
+ }
+ else{
+ vPosition.insert(w);
+ gf->mesh_vertices.push_back(w);
+ }
}
}
for(int i = 0; i < M - 1; i++){
@@ -41,6 +114,9 @@ public:
}
};
+std::set<meshCartesian::DVertex*, MVertexLessThanLexicographic> meshCartesian::vPosition;
+std::set<MVertex*> meshCartesian::vDelete;
+
class computePartitionOfUnity{
public:
void operator() (GFace *gf)
@@ -57,7 +133,7 @@ public:
SPoint2 uv2 = gf2->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
if(gf2->containsParam(uv2)){
GPoint xyz2 = gf2->point(uv2);
- const double tol = 1.e-2; // need to adapt this w.r.t size of model
+ const double tol = 1.e-2; // Warning: tolerance
if(fabs(xyz2.x() - v->x()) < tol &&
fabs(xyz2.y() - v->y()) < tol &&
fabs(xyz2.z() - v->z()) < tol)
@@ -87,10 +163,11 @@ public:
}
};
-class createGrooves{
+class createGroove{
public:
void operator() (GFace *gf)
{
+ // mark elements in the groove with '-1' visibility tag
for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
MElement *e = gf->quadrangles[i];
for(int j = 0; j < e->getNumVertices(); j++){
@@ -98,23 +175,39 @@ public:
if(data){
double pou = *(double*)data;
if(pou < 0.51)
- e->setVisibility(0);
+ e->setVisibility(-1);
}
}
}
- std::vector<MQuadrangle*> remaining;
+
+ // remove elements in the groove
+ std::vector<MQuadrangle*> newq;
for(unsigned int i = 0; i < gf->quadrangles.size(); i++)
- if(gf->quadrangles[i]->getVisibility())
- remaining.push_back(gf->quadrangles[i]);
- else
+ if(gf->quadrangles[i]->getVisibility() < 0)
delete gf->quadrangles[i];
- gf->quadrangles.clear();
- gf->quadrangles = remaining;
+ else
+ newq.push_back(gf->quadrangles[i]);
+ gf->quadrangles = newq;
+
+ // remove vertices in the groove (we cannot delete them right away
+ // since some can be shared between faces on hard edges)
+ std::vector<MVertex*> newv;
+ for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++)
+ gf->mesh_vertices[i]->setVisibility(-1);
+ for(unsigned int i = 0; i < gf->quadrangles.size(); i++)
+ for(int j = 0; j < gf->quadrangles[i]->getNumVertices(); j++)
+ gf->quadrangles[i]->getVertex(j)->setVisibility(1);
+ for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++)
+ if(gf->mesh_vertices[i]->getVisibility() < 0)
+ meshCartesian::vDelete.insert(gf->mesh_vertices[i]);
+ else
+ newv.push_back(gf->mesh_vertices[i]);
+ gf->mesh_vertices = newv;
}
};
-void getOrderedBoundaryVertices(std::vector<MElement*> &elements,
- std::vector<MVertex*> &boundary)
+void getOrderedBoundaryLoops(std::vector<MElement*> &elements,
+ std::vector<std::vector<MVertex*> > &loops)
{
typedef std::pair<MVertex*, MVertex*> vpair;
std::map<vpair, vpair> edges;
@@ -126,12 +219,27 @@ void getOrderedBoundaryVertices(std::vector<MElement*> &elements,
else edges[p] = vpair(e.getVertex(0), e.getVertex(1));
}
}
+
std::map<MVertex*, MVertex*> connect;
for(std::map<vpair, vpair>::iterator it = edges.begin(); it != edges.end(); it++)
connect[it->second.first] = it->second.second;
- boundary.push_back(connect.begin()->first);
- for(unsigned int i = 0; i < connect.size(); i++)
- boundary.push_back(connect[boundary[boundary.size() - 1]]);
+
+ loops.resize(1);
+ while(connect.size()){
+ if(loops[loops.size() - 1].empty())
+ loops[loops.size() - 1].push_back(connect.begin()->first);
+ MVertex *prev = loops[loops.size() - 1][loops[loops.size() - 1].size() - 1];
+ MVertex *next = connect[prev];
+ connect.erase(prev);
+ loops[loops.size() - 1].push_back(next);
+ if(next == loops[loops.size() - 1][0])
+ loops.resize(loops.size() + 1);
+ }
+
+ if(loops[loops.size() - 1].empty())
+ loops.pop_back();
+
+ Msg(INFO, "Found %d loop(s) in boundary", loops.size());
}
void addElements(GFace *gf, std::vector<MElement*> &elements)
@@ -142,55 +250,78 @@ void addElements(GFace *gf, std::vector<MElement*> &elements)
elements.push_back(gf->quadrangles[i]);
}
-bool isConnected(GFace *gf1, GFace *gf2)
+void classifyFaces(GFace *gf, std::set<GFace*> &connected, std::set<GFace*> &other)
{
- std::set<MVertex*> set;
- std::vector<MElement*> elements1, elements2;
- std::vector<MVertex*> boundary1, boundary2;
- addElements(gf1, elements1);
- addElements(gf2, elements2);
- getOrderedBoundaryVertices(elements1, boundary1);
- getOrderedBoundaryVertices(elements2, boundary2);
-
- for(unsigned int i = 0; i < boundary1.size(); i++)
- set.insert(boundary1[i]);
- for(unsigned int i = 0; i < boundary2.size(); i++){
- std::set<MVertex*>::iterator it = set.find(boundary2[i]);
- if(it != set.end()) return true;
+ connected.insert(gf);
+ std::set<MVertex*> verts;
+ for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++)
+ verts.insert(gf->mesh_vertices[i]);
+
+ for(int i = 0; i < gf->model()->numFace() - 1; i++){ // worst case
+ for(GModel::fiter it = gf->model()->firstFace(); it != gf->model()->lastFace(); it++){
+ if(connected.find(*it) == connected.end()){
+ for(unsigned int j = 0; j < (*it)->mesh_vertices.size(); j++){
+ if(std::find(verts.begin(), verts.end(), (*it)->mesh_vertices[j]) != verts.end()){
+ connected.insert(*it);
+ for(unsigned int k = 0; k < (*it)->mesh_vertices.size(); k++)
+ verts.insert((*it)->mesh_vertices[k]);
+ break;
+ }
+ }
+ }
+ }
}
- return false;
+
+ for(GModel::fiter it = gf->model()->firstFace(); it != gf->model()->lastFace(); it++)
+ if(connected.find(*it) == connected.end())
+ other.insert(*it);
+
+ Msg(INFO, "Found %d face(s) connected to face %d", (int)connected.size(), gf->tag());
+ for(std::set<GFace*>::iterator it = connected.begin(); it != connected.end(); it++)
+ Msg(INFO, " %d", (*it)->tag());
}
-void getIntersectingBoundaryParts(GFace *gf, std::vector<MElement*> &elements, int &start,
- std::map<int, std::vector<MVertex*> > &parts)
+void getIntersectingBoundaryParts(GFace *gf, std::vector<MElement*> &elements,
+ std::vector<std::vector<std::vector<MVertex*> > > &parts)
{
- bool newpart = true;
- std::vector<MVertex*> vertices;
- getOrderedBoundaryVertices(elements, vertices);
- for(unsigned int i = 0; i < vertices.size() - 1; i++){
- MVertex *v = vertices[i];
- SPoint2 uv = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
- if(gf->containsParam(uv)){
- GPoint xyz = gf->point(uv);
- const double tol = 1.e-2; // need to adapt this w.r.t size of model
- if(fabs(xyz.x() - v->x()) < tol &&
- fabs(xyz.y() - v->y()) < tol &&
- fabs(xyz.z() - v->z()) < tol){
- if(newpart){
- start++;
- newpart = false;
+ std::vector<std::vector<MVertex*> > loops;
+ getOrderedBoundaryLoops(elements, loops);
+ parts.resize(loops.size());
+
+ if(gf->tag() == 0){
+ debugElements(elements, "elements", false);
+ for(unsigned int i = 0; i < loops.size(); i++)
+ debugVertices(loops[i], "boundary", false, i);
+ }
+
+ for(unsigned int i = 0; i < loops.size(); i++){
+ // last vertex in loop is equal to the firt vertex
+ bool newpart = true;
+ for(unsigned int j = 0; j < loops[i].size() - 1; j++){
+ MVertex *v = loops[i][j];
+ SPoint2 uv = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+ if(gf->containsParam(uv)){
+ GPoint xyz = gf->point(uv);
+ const double tol = 1.e-2; // Warning: tolerance
+ if(fabs(xyz.x() - v->x()) < tol &&
+ fabs(xyz.y() - v->y()) < tol &&
+ fabs(xyz.z() - v->z()) < tol){
+ if(newpart){
+ parts[i].resize(parts[i].size() + 1);
+ newpart = false;
+ }
+ v->setParameter(0, uv[0]);
+ v->setParameter(1, uv[1]);
+ parts[i][parts[i].size() - 1].push_back(v);
+ }
+ else{
+ newpart = true;
}
- v->setParameter(0, uv[0]);
- v->setParameter(1, uv[1]);
- parts[start].push_back(v);
}
else{
newpart = true;
}
}
- else{
- newpart = true;
- }
}
}
@@ -199,169 +330,122 @@ public:
void operator() (GFace *gf)
{
if(gf->tag() > 1) return;
- printf("processing grout for face %d\n", gf->tag());
-
- GModel *m = gf->model();
-
- // need to find the disconnected parts in inside too!!! then
- // we'll pair them with disconnected parts in outside (by looking
- // at the value of the POU)
- // Distinguish 2 cases:
- // - patch with no hard edges and and no connections: we have a hole
- // THIS IS WHAT WE ASSUME NOW
- // - patch with hard edges or connections: we don't have a hole, but
- // we might have non-connected parts
- std::vector<MElement*> elements;
-
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
- if(isConnected(gf, *it)){
- printf("face %d: mesh connected to %d\n", gf->tag(), (*it)->tag());
- addElements(*it, elements);
- }
- }
- int numinside = 0;
- std::map<int, std::vector<MVertex*> > inside;
- getIntersectingBoundaryParts(gf, elements, numinside, inside);
-
- printf("numinside = %d\n", numinside);
-
- int numoutside = 0;
- std::map<int, std::vector<MVertex*> > outside;
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
- if(!isConnected(gf, *it)){
- elements.clear();
- addElements(*it, elements);
- getIntersectingBoundaryParts(gf, elements, numoutside, outside);
- }
- }
- printf("numoutside = %d\n", numoutside);
+ Msg(INFO, "Processing grout for face %d", gf->tag());
+
+ std::set<GFace*> connected, other;
+ classifyFaces(gf, connected, other);
+
+ std::vector<MElement*> connectedElements, otherElements;
+ for(std::set<GFace*>::iterator it = connected.begin(); it != connected.end(); it++)
+ addElements(*it, connectedElements);
+ for(std::set<GFace*>::iterator it = other.begin(); it != other.end(); it++)
+ addElements(*it, otherElements);
+
+ std::vector<std::vector<std::vector<MVertex*> > > inside;
+ getIntersectingBoundaryParts(gf, connectedElements, inside);
+ Msg(INFO, "%d inside loop(s)", (int)inside.size());
+ for(unsigned int i = 0; i < inside.size(); i++)
+ Msg(INFO, " inside loop %d intersect has %d part(s)", i, (int)inside[i].size());
+
+ std::vector<std::vector<std::vector<MVertex*> > > outside;
+ getIntersectingBoundaryParts(gf, otherElements, outside);
+ Msg(INFO, "%d outside loop(s)", (int)outside.size());
+ for(unsigned int i = 0; i < outside.size(); i++)
+ Msg(INFO, " outside loop %d intersect has %d part(s)", i, (int)outside[i].size());
std::vector<MVertex*> hole, loop;
- if(numinside == 1){
+ if(inside.size() == 1 && inside[0].size() == 1){
// create hole
SPoint2 ic(0., 0.);
- for(unsigned int i = 0; i < inside[1].size(); i++){
- hole.push_back(inside[1][i]);
- double u, v;
- inside[1][i]->getParameter(0, u);
- inside[1][i]->getParameter(1, v);
- ic += SPoint2(u, v);
- }
- ic *= 1. / (double) inside[1].size();
- hole.push_back(hole[0]);
-
- // order exterior parts
- std::vector<std::pair<double, int> > angle;
- std::map<int, std::vector<MVertex*> >::iterator it = outside.begin();
- for(; it != outside.end(); it++){
- SPoint2 oc(0., 0.);
- for(unsigned int i = 0; i < it->second.size(); i++){
+ {
+ for(unsigned int i = 0; i < inside[0][0].size(); i++){
+ hole.push_back(inside[0][0][i]);
double u, v;
- it->second[i]->getParameter(0, u);
- it->second[i]->getParameter(1, v);
- oc += SPoint2(u, v);
+ inside[0][0][i]->getParameter(0, u);
+ inside[0][0][i]->getParameter(1, v);
+ ic += SPoint2(u, v);
}
- oc *= 1. / (double)it->second.size();
- double a = atan2(oc[1] - ic[1], oc[0] - ic[0]);
- angle.push_back(std::make_pair(a, it->first));
+ ic *= 1. / (double)inside[0][0].size();
+ hole.push_back(hole[0]);
}
- std::sort(angle.begin(), angle.end());
-
- // create exterior loop
- for(unsigned int i = 0; i < angle.size(); i++){
- for(unsigned int j = 0; j < outside[angle[i].second].size(); j++)
- loop.push_back(outside[angle[i].second][j]);
+ // order exterior parts and create exterior loop
+ {
+ std::vector<std::pair<double, int> > angle;
+ std::map<int, std::vector<MVertex*> > outside_flat;
+ int part = 0;
+ for(unsigned int i = 0; i < outside.size(); i++){
+ for(unsigned int j = 0; j < outside[i].size(); j++){
+ for(unsigned int k = 0; k < outside[i][j].size(); k++){
+ outside_flat[part].push_back(outside[i][j][k]);
+ }
+ part++;
+ }
+ }
+ std::map<int, std::vector<MVertex*> >::iterator it = outside_flat.begin();
+ for(; it != outside_flat.end(); it++){
+ SPoint2 oc(0., 0.);
+ for(unsigned int i = 0; i < it->second.size(); i++){
+ double u, v;
+ it->second[i]->getParameter(0, u);
+ it->second[i]->getParameter(1, v);
+ oc += SPoint2(u, v);
+ }
+ oc *= 1. / (double)it->second.size();
+ double a = atan2(oc[1] - ic[1], oc[0] - ic[0]);
+ angle.push_back(std::make_pair(a, it->first));
+ }
+ std::sort(angle.begin(), angle.end());
+ for(unsigned int i = 0; i < angle.size(); i++){
+ for(unsigned int j = 0; j < outside_flat[angle[i].second].size(); j++)
+ loop.push_back(outside_flat[angle[i].second][j]);
+ }
+ loop.push_back(hole[0]);
}
- loop.push_back(hole[0]);
-
// mesh the grout
fourierFace *grout = new fourierFace(gf, loop, hole);
meshGFace mesh;
mesh(grout);
for(unsigned int i = 0; i < grout->triangles.size(); i++)
gf->triangles.push_back(grout->triangles[i]);
+ for(unsigned int i = 0; i < loop.size(); i++)
+ gf->mesh_vertices.push_back(loop[i]);
+ for(unsigned int i = 0; i < hole.size(); i++)
+ gf->mesh_vertices.push_back(hole[i]);
for(unsigned int i = 0; i < grout->mesh_vertices.size(); i++)
gf->mesh_vertices.push_back(grout->mesh_vertices[i]);
delete grout;
}
else{
-
Msg(WARNING, "Faces with no holes not implemented yet!");
- }
+ // num individual meshes = num parts with onWhat() == gf!
- // debug
- char name[256];
- sprintf(name, "aa%d.pos", gf->tag());
- FILE *fp = fopen(name, "w");
- fprintf(fp, "View \"aa\"{\n");
- for(unsigned int i = 0; i < hole.size(); i++){
- double x = hole[i]->x(), y = hole[i]->y(), z = hole[i]->z();
- // plot in parametric space:
- hole[i]->getParameter(0, x); hole[i]->getParameter(1, y); z = 0;
- fprintf(fp, "SP(%g,%g,%g){0};\n", x, y, z);
- }
- for(unsigned int i = 0; i < loop.size(); i++){
- double x = loop[i]->x(), y = loop[i]->y(), z = loop[i]->z();
- // plot in parametric space:
- loop[i]->getParameter(0, x); loop[i]->getParameter(1, y); z = 0;
- fprintf(fp, "SP(%g,%g,%g){%d};\n", x, y, z, i);
- }
- fprintf(fp, "};\n");
- fclose(fp);
- }
-};
+ // for each one
+ // - find other parts that are "close" (using POUs?)
+ // - sort parts w.r.t barycenter of each group
-class exportFourierFace{
-public:
- void operator() (GFace *gf)
- {
- Post_View *view = BeginView(1);
- for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
- MElement *e = gf->quadrangles[i];
- double x[4], y[4], z[4], val[4];
- for(int j = 0; j < 4; j++){
- MVertex *v = e->getVertex(j);
- x[j] = v->x();
- y[j] = v->y();
- z[j] = v->z();
- val[j] = *(double*)v->getData();
+ for(unsigned int i = 0; i < inside.size(); i++){
+ for(unsigned int j = 0; j < inside[i].size(); j++){
+ for(unsigned int k = 0; k < inside[i][j].size(); k++){
+ loop.push_back(inside[i][j][k]);
+ }
+ }
}
- for(int j = 0; j < 4; j++) List_Add(view->SQ, &x[j]);
- for(int j = 0; j < 4; j++) List_Add(view->SQ, &y[j]);
- for(int j = 0; j < 4; j++) List_Add(view->SQ, &z[j]);
- for(int j = 0; j < 4; j++) List_Add(view->SQ, &val[j]);
- view->NbSQ++;
}
- std::vector<MElement*> elements;
- std::vector<MVertex*> vertices;
- addElements(gf, elements);
- getOrderedBoundaryVertices(elements, vertices);
- for(unsigned int i = 0; i < vertices.size() - 1; i++){
- double x[2] = {vertices[i]->x(), vertices[i + 1]->x()};
- double y[2] = {vertices[i]->y(), vertices[i + 1]->y()};
- double z[2] = {vertices[i]->z(), vertices[i + 1]->z()};
- double val[2] = {10, 10};
- for(int j = 0; j < 2; j++) List_Add(view->SL, &x[j]);
- for(int j = 0; j < 2; j++) List_Add(view->SL, &y[j]);
- for(int j = 0; j < 2; j++) List_Add(view->SL, &z[j]);
- for(int j = 0; j < 2; j++) List_Add(view->SL, &val[j]);
- view->NbSL++;
- }
- char name[1024], filename[1024];
- sprintf(name, "patch%d", gf->tag());
- sprintf(filename, "patch%d", gf->tag());
- EndView(view, 1, filename, name);
+
+ debugVertices(hole, "hole", false, gf->tag());
+ debugVertices(loop, "loop", false, gf->tag());
}
};
-
fourierModel::fourierModel(const std::string &name)
: GModel(name)
{
FM = new model(name);
+
+ CTX.terminal = 1;
Msg(INFO, "Fourier model created: %d patches", FM->GetNumPatches());
@@ -369,6 +453,9 @@ fourierModel::fourierModel(const std::string &name)
for(int i = 0; i < FM->GetNumPatches(); i++)
add(new fourierFace(this, i));
+ meshCartesian::vPosition.clear();
+ MVertexLessThanLexicographic::tolerance = 1.e-12; // Warning: tolerance
+
// mesh each face with quads
std::for_each(firstFace(), lastFace(), meshCartesian());
@@ -376,14 +463,13 @@ fourierModel::fourierModel(const std::string &name)
std::for_each(firstFace(), lastFace(), computePartitionOfUnity());
// create grooves
- std::for_each(firstFace(), lastFace(), createGrooves());
+ meshCartesian::vDelete.clear();
+ std::for_each(firstFace(), lastFace(), createGroove());
+ meshCartesian::deleteUnusedVertices();
// create grout
std::for_each(firstFace(), lastFace(), createGrout());
- // visualize as a post-pro view
- //std::for_each(firstFace(), lastFace(), exportFourierFace());
-
CTX.mesh.changed = ENT_ALL;
}