diff --git a/Geo/GModelIO.cpp b/Geo/GModelIO.cpp
index 97c06dfd66e84d2614bc4960cf6a29bb5f6917df..bbb37ef6e7182f30e6211d6b1cd3ed13e47766d5 100644
--- a/Geo/GModelIO.cpp
+++ b/Geo/GModelIO.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO.cpp,v 1.49 2006-09-09 01:10:05 geuzaine Exp $
+// $Id: GModelIO.cpp,v 1.50 2006-09-10 01:49:31 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -1073,7 +1073,6 @@ int GModel::readUNV(const std::string &name)
char buffer[256];
std::map<int, MVertex*> vertexMap;
- std::map<int, std::vector<MVertex*> > points;
std::map<int, std::vector<MElement*> > elements[7];
std::map<int, std::map<int, std::string> > physicals[4];
@@ -1294,8 +1293,6 @@ int GModel::readMESH(const std::string &name)
std::vector<MVertex*> vertexVector;
std::map<int, std::vector<MElement*> > elements[3];
- int nbv, nbe, n[30], cl;
-
while(!feof(fp)) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
if(buffer[0] != '#'){ // skip comments
@@ -1305,22 +1302,26 @@ int GModel::readMESH(const std::string &name)
}
else if(!strcmp(str, "Vertices")){
if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int nbv;
sscanf(buffer, "%d", &nbv);
Msg(INFO, "%d vertices", nbv);
vertexVector.resize(nbv);
for(int i = 0; i < nbv; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int dum;
double x, y, z;
- sscanf(buffer, "%lf %lf %lf %d", &x, &y, &z, &cl);
+ sscanf(buffer, "%lf %lf %lf %d", &x, &y, &z, &dum);
vertexVector[i] = new MVertex(x, y, z);
}
}
else if(!strcmp(str, "Triangles")){
if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int nbe;
sscanf(buffer, "%d", &nbe);
Msg(INFO, "%d triangles", nbe);
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int n[3], cl;
sscanf(buffer, "%d %d %d %d", &n[0], &n[1], &n[2], &cl);
for(int j = 0; j < 3; j++) n[j]--;
std::vector<MVertex*> vertices;
@@ -1330,10 +1331,12 @@ int GModel::readMESH(const std::string &name)
}
else if(!strcmp(str, "Quadrilaterals")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int nbe;
sscanf(buffer, "%d", &nbe);
Msg(INFO, "%d quadrangles", nbe);
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int n[4], cl;
sscanf(buffer, "%d %d %d %d %d", &n[0], &n[1], &n[2], &n[3], &cl);
for(int j = 0; j < 4; j++) n[j]--;
std::vector<MVertex*> vertices;
@@ -1343,10 +1346,12 @@ int GModel::readMESH(const std::string &name)
}
else if(!strcmp(str, "Tetrahedra")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int nbe;
sscanf(buffer, "%d", &nbe);
Msg(INFO, "%d tetrahedra", nbe);
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int n[4], cl;
sscanf(buffer, "%d %d %d %d %d", &n[0], &n[1], &n[2], &n[3], &cl);
for(int j = 0; j < 4; j++) n[j]--;
std::vector<MVertex*> vertices;
@@ -1564,7 +1569,7 @@ static int readElementBDF(FILE *fp, char *buffer, int keySize, int numVertices,
strncpy(tmp, fields[i], cmax); n[i - 2] = atoi(tmp);
}
- // ignore the extra fields when we now how many vertices we need
+ // ignore the extra fields when we know how many vertices we need
int numCheck = (numVertices > 0) ? numVertices : fields.size() - 2;
if(!getVertices(numCheck, n, vertexMap, vertices)) return 0;
return 1;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 3e3886a2ee6ea3740b5a2eeb7627923b5407b436..2745c7b4b39c520259a2352b090cb3dc3f51f223 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -238,12 +238,6 @@ class MTriangle : public MElement {
virtual int getDim(){ return 2; }
virtual int getNumVertices(){ return 3; }
virtual MVertex *getVertex(int num){ return _v[num]; }
- virtual void revert ()
- {
- MVertex *vv = _v[0];
- _v[0] = _v[1];
- _v[1] = vv;
- }
virtual int getNumEdges(){ return 3; }
virtual MEdge getEdge(int num)
{
@@ -258,6 +252,12 @@ class MTriangle : public MElement {
virtual int getTypeForUNV(){ return 91; } // thin shell linear triangle
virtual char *getStringForPOS(){ return "ST"; }
virtual char *getStringForBDF(){ return "CTRIA3"; }
+ virtual void revert()
+ {
+ MVertex *vv = _v[0];
+ _v[0] = _v[1];
+ _v[1] = vv;
+ }
};
class MTriangle6 : public MTriangle {
@@ -348,6 +348,12 @@ class MQuadrangle : public MElement {
virtual int getTypeForUNV(){ return 94; } // thin shell linear quadrilateral
virtual char *getStringForPOS(){ return "SQ"; }
virtual char *getStringForBDF(){ return "CQUAD4"; }
+ virtual void revert()
+ {
+ MVertex *vv = _v[1];
+ _v[1] = _v[3];
+ _v[3] = vv;
+ }
};
class MQuadrangle8 : public MQuadrangle {
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 101430ed72095b10cba05d7285a1f2459ac14a08..ce8da2b242d98d11dd0b68dd58e1fd8f25c5f933 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -73,8 +73,9 @@ class MVertex{
inline int getNum() const {return _num;}
inline void setNum(int num) { _num = num; }
- // Get ith parameter
+ // get/set ith parameter
virtual bool getParameter(int i, double &par){ return false; }
+ virtual bool setParameter(int i, double par){ return false; }
// Get the data associated with this vertex
virtual void *getData(){ return 0; }
@@ -99,6 +100,7 @@ class MEdgeVertex : public MVertex{
}
~MEdgeVertex(){}
virtual bool getParameter(int i, double &par){ par = _u; return true; }
+ virtual bool setParameter(int i, double par){ _u = par; return true; }
};
class MFaceVertex : public MVertex{
@@ -111,6 +113,7 @@ class MFaceVertex : public MVertex{
}
~MFaceVertex(){}
virtual bool getParameter(int i, double &par){ par = (i ? _v : _u); return true; }
+ virtual bool setParameter(int i, double par){ if(!i) _u = par; else _v = par; return true; }
};
template<class T>
diff --git a/Geo/fourierModel.cpp b/Geo/fourierModel.cpp
index e39364e2932456b7b2162f03a6e09b67cbd78936..27d4289998a94cb3390b474311d5cbc5204d3227 100644
--- a/Geo/fourierModel.cpp
+++ b/Geo/fourierModel.cpp
@@ -8,6 +8,7 @@ extern Context_T CTX;
#if defined(HAVE_FOURIER_MODEL)
#include "model.h"
+#include "meshGFace.h"
static model *FM = 0;
@@ -29,10 +30,12 @@ public:
}
for(int i = 0; i < M - 1; i++){
for(int j = 0; j < N - 1; j++){
- gf->quadrangles.push_back(new MQuadrangle(gf->mesh_vertices[i * N + j],
- gf->mesh_vertices[i * N + (j + 1)],
- gf->mesh_vertices[(i + 1) * N + (j + 1)],
- gf->mesh_vertices[(i + 1) * N + j]));
+ MQuadrangle *q = new MQuadrangle(gf->mesh_vertices[i * N + j],
+ gf->mesh_vertices[(i + 1) * N + j],
+ gf->mesh_vertices[(i + 1) * N + (j + 1)],
+ gf->mesh_vertices[i * N + (j + 1)]);
+ if(FM->GetOrientation(gf->tag()) < 0) q->revert();
+ gf->quadrangles.push_back(q);
}
}
}
@@ -94,24 +97,173 @@ public:
void *data = e->getVertex(j)->getData();
if(data){
double pou = *(double*)data;
- if(pou < 0.51) e->setVisibility(0);
+ if(pou < 0.51)
+ e->setVisibility(0);
}
}
}
}
};
+void getOrderedBoundaryVertices(std::vector<MElement*> &elements,
+ std::vector<MVertex*> &boundary)
+{
+ typedef std::pair<MVertex*, MVertex*> vpair;
+ std::map<vpair, vpair> edges;
+ for(unsigned int i = 0; i < elements.size(); i++){
+ if(elements[i]->getVisibility() > 0){
+ for(int j = 0; j < elements[i]->getNumEdges(); j++){
+ MEdge e = elements[i]->getEdge(j);
+ vpair p(e.getMinVertex(), e.getMaxVertex());
+ if(edges.count(p)) edges.erase(p);
+ 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]]);
+}
+
+void addElements(GFace *gf, std::vector<MElement*> &elements)
+{
+ for(unsigned int i = 0; i < gf->triangles.size(); i++)
+ elements.push_back(gf->triangles[i]);
+ for(unsigned int i = 0; i < gf->quadrangles.size(); i++)
+ elements.push_back(gf->quadrangles[i]);
+}
+
+bool isConnected(GFace *gf1, GFace *gf2)
+{
+ 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;
+ }
+ return false;
+}
+
class createGrout{
public:
void operator() (GFace *gf)
{
+ if(gf->tag() != 0) return;
printf("processing grout for face %d\n", gf->tag());
std::vector<int> overlaps;
FM->GetNeighbor(gf->tag(), overlaps);
+
+ std::vector<MElement*> elements;
+ for(unsigned int i = 0; i < overlaps.size(); i++){
+ GFace *gf2 = gf->model()->faceByTag(overlaps[i]);
+ if(isConnected(gf, gf2))
+ addElements(gf2, elements);
+ }
+
+ std::vector<MVertex*> inside;
+ getOrderedBoundaryVertices(elements, inside);
+
+ std::map<int, std::vector<MVertex*> > outsidePart;
+ int part = 0;
for(unsigned int i = 0; i < overlaps.size(); i++){
GFace *gf2 = gf->model()->faceByTag(overlaps[i]);
+ bool newpart = true;
+ if(!isConnected(gf, gf2)){
+ elements.clear();
+ addElements(gf2, elements);
+ 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){
+ part++;
+ newpart = false;
+ }
+ v->setParameter(0, uv[0]);
+ v->setParameter(1, uv[1]);
+ outsidePart[part].push_back(v);
+ }
+ else{
+ newpart = true;
+ }
+ }
+ else{
+ newpart = true;
+ }
+ }
+ }
+ }
+
+ // order exterior parts
+ std::vector<std::pair<double, int> > angle;
+ for(std::map<int, std::vector<MVertex*> >::iterator it = outsidePart.begin();
+ it != outsidePart.end(); it++){
+ SPoint2 barycenter(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);
+ barycenter += SPoint2(u, v);
+ }
+ barycenter *= 1. / (double)it->second.size();
+ double a = atan2(barycenter[1] - 0.5, barycenter[0] - 0.5);
+ angle.push_back(std::make_pair(a, it->first));
+ }
+ std::sort(angle.begin(), angle.end());
+
+ std::vector<MVertex*> outside;
+ for(unsigned int i = 0; i < angle.size(); i++){
+ for(unsigned int j = 0; j < outsidePart[angle[i].second].size(); j++)
+ outside.push_back(outsidePart[angle[i].second][j]);
}
+ outside.push_back(outside[0]);
+
+ // mesh the grout
+ fourierFace *tmp = new fourierFace(gf, outside, inside);
+ meshGFace mesh;
+ mesh(tmp);
+ for(unsigned int i = 0; i < tmp->triangles.size(); i++)
+ gf->triangles.push_back(tmp->triangles[i]);
+
+ // mesh and store elements in gf
+
+ 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 < inside.size(); i++){
+ double x = inside[i]->x(), y = inside[i]->y(), z = inside[i]->z();
+ // plot in parametric space:
+ inside[i]->getParameter(0, x); inside[i]->getParameter(1, y); z = 0;
+ fprintf(fp, "SP(%g,%g,%g){0};\n", x, y, z);
+ }
+ for(unsigned int i = 0; i < outside.size(); i++){
+ double x = outside[i]->x(), y = outside[i]->y(), z = outside[i]->z();
+ // plot in parametric space:
+ outside[i]->getParameter(0, x); outside[i]->getParameter(1, y); z = 0;
+ fprintf(fp, "SP(%g,%g,%g){%d};\n", x, y, z, i);
+ }
+ fprintf(fp, "};\n");
+ fclose(fp);
}
};
@@ -137,6 +289,21 @@ public:
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());
@@ -169,7 +336,7 @@ fourierModel::fourierModel(const std::string &name)
std::for_each(firstFace(), lastFace(), createGrout());
// visualize as a post-pro view
- std::for_each(firstFace(), lastFace(), exportFourierFace());
+ //std::for_each(firstFace(), lastFace(), exportFourierFace());
CTX.mesh.changed = ENT_ALL;
}
@@ -181,8 +348,59 @@ fourierModel::~fourierModel()
FM = 0;
}
+fourierEdge::fourierEdge(GModel *model, int num, GVertex *v1, GVertex *v2)
+ : GEdge(model, num, v1, v2)
+{
+}
+
fourierFace::fourierFace(GModel *m, int num)
: GFace(m, num)
+{
+ _discrete = 1;
+}
+
+fourierFace::fourierFace(GFace *f, std::vector<MVertex*> &loop, std::vector<MVertex*> &hole)
+ : GFace(f->model(), f->tag())
+{
+ _discrete = 0;
+
+ if(!loop.size()){
+ Msg(GERROR, "No vertices in exterior loop");
+ return;
+ }
+
+ fourierVertex *v1 = new fourierVertex(f->model(), loop[0]);
+ fourierVertex *v2 = new fourierVertex(f->model(), loop[loop.size() - 2]);
+
+ fourierEdge *e1 = new fourierEdge(f->model(), 1, v1, v2);
+ e1->addFace(this);
+ for(unsigned int i = 1; i < loop.size() - 2; i++)
+ e1->mesh_vertices.push_back(loop[i]);
+
+ fourierEdge *e2 = new fourierEdge(f->model(), 2, v2, v1);
+ e2->addFace(this);
+
+ l_edges.push_back(e1); l_dirs.push_back(1);
+ l_edges.push_back(e2); l_dirs.push_back(1);
+
+ if(hole.size()){
+ fourierVertex *v3 = new fourierVertex(f->model(), hole[0]);
+ fourierVertex *v4 = new fourierVertex(f->model(), hole[hole.size() - 2]);
+
+ fourierEdge *e3 = new fourierEdge(f->model(), 3, v3, v4);
+ e3->addFace(this);
+ for(unsigned int i = 1; i < hole.size() - 2; i++)
+ e3->mesh_vertices.push_back(hole[i]);
+
+ fourierEdge *e4 = new fourierEdge(f->model(), 4, v4, v3);
+ e3->addFace(this);
+
+ l_edges.push_back(e3); l_dirs.push_back(1);
+ l_edges.push_back(e4); l_dirs.push_back(1);
+ }
+}
+
+fourierFace::~fourierFace()
{
}
@@ -233,7 +451,7 @@ SVector3 fourierFace::normal(const SPoint2 ¶m) const
GEntity::GeomType fourierFace::geomType() const
{
- return GEntity::DiscreteSurface;
+ return _discrete ? GEntity::DiscreteSurface : GEntity::Unknown;
}
SPoint2 fourierFace::parFromPoint(const SPoint3 &p) const
diff --git a/Geo/fourierModel.h b/Geo/fourierModel.h
index 051ee0a6bf0ae50c1844ededd2438063a2c82dca..6111ae8a71d28489f1985d83416636eb236d901f 100644
--- a/Geo/fourierModel.h
+++ b/Geo/fourierModel.h
@@ -11,15 +11,58 @@ class fourierModel : public GModel {
virtual ~fourierModel();
};
+#include "GVertex.h"
+#include "GEdge.h"
#include "GFace.h"
#include "Range.h"
+class fourierVertex : public GVertex {
+ private:
+ MVertex *_v;
+ public:
+ fourierVertex(GModel *m, MVertex *v) : GVertex(m, v->getNum()), _v(v)
+ {
+ mesh_vertices.push_back(v);
+ }
+ virtual ~fourierVertex() {}
+ virtual GPoint point() const { return GPoint(_v->x(), _v->y(), _v->z(), this); }
+ virtual double x() const { return _v->x(); }
+ virtual double y() const { return _v->y(); }
+ virtual double z() const { return _v->z(); }
+ virtual double prescribedMeshSizeAtVertex() const { return 0.1; }
+};
+
+class fourierEdge : public GEdge {
+ private:
+ GVertex *v0, *v1;
+ public:
+ fourierEdge(GModel *m, int num, GVertex *v1, GVertex *v2);
+ virtual ~fourierEdge() {}
+ double period() const { throw ; }
+ virtual bool periodic(int dim=0) const { return false; }
+ virtual Range<double> parBounds(int i) const { throw; }
+ virtual GeomType geomType() const { throw; }
+ virtual bool degenerate(int) const { return false; }
+ virtual bool continuous(int dim) const { return true; }
+ virtual GPoint point(double p) const { throw; }
+ virtual GPoint closestPoint(const SPoint3 & queryPoint) { throw; }
+ virtual int containsPoint(const SPoint3 &pt) const { throw; }
+ virtual int containsParam(double pt) const { throw; }
+ virtual SVector3 firstDer(double par) const { throw; }
+ virtual SPoint2 reparamOnFace(GFace * face, double epar, int dir) const { throw; }
+ virtual double parFromPoint(const SPoint3 &pt) const { throw; }
+ virtual int minimumMeshSegments () const { throw; }
+ virtual int minimumDrawSegments () const { throw; }
+};
+
class fourierFace : public GFace {
private:
int _num;
+ int _discrete;
public:
fourierFace(GModel *m, int num);
- virtual ~fourierFace(){}
+ fourierFace(GFace *f, std::vector<MVertex*> &loop, std::vector<MVertex*> &hole);
+ virtual ~fourierFace();
Range<double> parBounds(int i) const;
virtual int paramDegeneracies(int dir, double *par) { return 0; }