diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 30d39d882004cb6518b3513c5273a57a9876143d..21cab6a447b28693b3c323f34db383fa984c13fc 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.18 2006-09-08 17:45:51 geuzaine Exp $
+// $Id: MElement.cpp,v 1.19 2006-09-10 15:34:12 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -121,36 +121,20 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
if(!binary){
fprintf(fp, "%d %d", num ? num : _num, type);
if(version < 2.0)
- fprintf(fp, " %d %d %d", physical, elementary, n);
+ fprintf(fp, " %d %d %d", abs(physical), elementary, n);
else
- fprintf(fp, " 3 %d %d %d", physical, elementary, _partition);
+ fprintf(fp, " 3 %d %d %d", abs(physical), elementary, _partition);
}
else{
- int tags[4] = {num ? num : _num, physical, elementary, _partition};
+ int tags[4] = {num ? num : _num, abs(physical), elementary, _partition};
fwrite(tags, sizeof(int), 4, fp);
}
- int verts[30];
+ if(physical < 0) revert();
- if(physical >= 0){
- for(int i = 0; i < n; i++)
- verts[i] = getVertex(i)->getNum();
- }
- else{
- int nn = n - getNumEdgeVertices() - getNumFaceVertices() - getNumVolumeVertices();
- int j = 0;
- for(int i = 0; i < nn; i++)
- verts[j++] = getVertex(nn - i - 1)->getNum();
- int ne = getNumEdgeVertices();
- for(int i = 0; i < ne; i++)
- verts[j++] = getVertex(nn + ne - i - 1)->getNum();
- int nf = getNumFaceVertices();
- for(int i = 0; i < nf; i++)
- verts[j++] = getVertex(nn + ne + nf - i - 1)->getNum();
- int nv = getNumVolumeVertices();
- for(int i = 0; i < nv; i++)
- verts[j++] = getVertex(n - i - 1)->getNum();
- }
+ int verts[30];
+ for(int i = 0; i < n; i++)
+ verts[i] = getVertex(i)->getNum();
if(!binary){
for(int i = 0; i < n; i++)
@@ -160,6 +144,8 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
else{
fwrite(verts, sizeof(int), n, fp);
}
+
+ if(physical < 0) revert();
}
void MElement::writePOS(FILE *fp, double scalingFactor, int elementary)
@@ -261,7 +247,7 @@ void MElement::writeUNV(FILE *fp, int num, int elementary, int physical)
setVolumePositive();
int n = getNumVertices();
int physical_property = elementary;
- int material_property = physical;
+ int material_property = abs(physical);
int color = 7;
fprintf(fp, "%10d%10d%10d%10d%10d%10d\n",
num ? num : _num, type, physical_property, material_property, color, n);
diff --git a/Geo/fourierModel.cpp b/Geo/fourierModel.cpp
index 27d4289998a94cb3390b474311d5cbc5204d3227..9972b464ce15d5267e7c2bcaa3b648e7109c42c0 100644
--- a/Geo/fourierModel.cpp
+++ b/Geo/fourierModel.cpp
@@ -12,11 +12,11 @@ extern Context_T CTX;
static model *FM = 0;
-class meshFourierFace{
+class meshCartesian{
public:
void operator() (GFace *gf)
{
- int M = 30, N = 30;
+ int M = (int)(30. / CTX.mesh.lc_factor), N = (int)(30. / CTX.mesh.lc_factor);
for(int i = 0; i < M; i++){
for(int j = 0; j < N; j++){
double u = i/(double)(M - 1);
@@ -159,26 +159,33 @@ class createGrout{
public:
void operator() (GFace *gf)
{
- if(gf->tag() != 0) return;
+ if(gf->tag() > 0) return;
printf("processing grout for face %d\n", gf->tag());
- std::vector<int> overlaps;
- FM->GetNeighbor(gf->tag(), overlaps);
+ 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(unsigned int i = 0; i < overlaps.size(); i++){
- GFace *gf2 = gf->model()->faceByTag(overlaps[i]);
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
+ GFace *gf2 = *it;
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]);
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
+ GFace *gf2 = *it;
bool newpart = true;
if(!isConnected(gf, gf2)){
elements.clear();
@@ -237,15 +244,20 @@ public:
}
outside.push_back(outside[0]);
+
+ if(gf->tag() == 0){
// mesh the grout
- fourierFace *tmp = new fourierFace(gf, outside, inside);
+ fourierFace *grout = 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
+ 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 < grout->mesh_vertices.size(); i++)
+ gf->mesh_vertices.push_back(grout->mesh_vertices[i]);
+ delete grout;
+ }
+ // debug
char name[256];
sprintf(name, "aa%d.pos", gf->tag());
FILE *fp = fopen(name, "w");
@@ -323,8 +335,8 @@ fourierModel::fourierModel(const std::string &name)
for(int i = 0; i < FM->GetNumPatches(); i++)
add(new fourierFace(this, i));
- // mesh each face
- std::for_each(firstFace(), lastFace(), meshFourierFace());
+ // mesh each face with quads
+ std::for_each(firstFace(), lastFace(), meshCartesian());
// compute partition of unity
std::for_each(firstFace(), lastFace(), computePartitionOfUnity());
@@ -356,12 +368,14 @@ fourierEdge::fourierEdge(GModel *model, int num, GVertex *v1, GVertex *v2)
fourierFace::fourierFace(GModel *m, int num)
: GFace(m, num)
{
+ for(int i = 0; i < 4; i++){ _v[i] = 0; _e[i] = 0; }
_discrete = 1;
}
fourierFace::fourierFace(GFace *f, std::vector<MVertex*> &loop, std::vector<MVertex*> &hole)
: GFace(f->model(), f->tag())
{
+ for(int i = 0; i < 4; i++){ _v[i] = 0; _e[i] = 0; }
_discrete = 0;
if(!loop.size()){
@@ -369,39 +383,56 @@ fourierFace::fourierFace(GFace *f, std::vector<MVertex*> &loop, std::vector<MVer
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);
+ _v[0] = new fourierVertex(f->model(), loop[0]);
+ _v[1] = new fourierVertex(f->model(), loop[loop.size() - 2]);
+ _e[0] = new fourierEdge(f->model(), 1, _v[0], _v[1]);
+ _e[0]->addFace(this);
+ _e[1] = new fourierEdge(f->model(), 2, _v[1], _v[0]);
+ _e[1]->addFace(this);
for(unsigned int i = 1; i < loop.size() - 2; i++)
- e1->mesh_vertices.push_back(loop[i]);
+ _e[0]->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);
+ l_edges.push_back(_e[0]); l_dirs.push_back(1);
+ l_edges.push_back(_e[1]); 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);
+ _v[2] = new fourierVertex(f->model(), hole[0]);
+ _v[3] = new fourierVertex(f->model(), hole[hole.size() - 2]);
+ _e[2] = new fourierEdge(f->model(), 3, _v[2], _v[3]);
+ _e[2]->addFace(this);
+ _e[3] = new fourierEdge(f->model(), 4, _v[3], _v[2]);
+ _e[3]->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);
+ _e[2]->mesh_vertices.push_back(hole[i]);
- l_edges.push_back(e3); l_dirs.push_back(1);
- l_edges.push_back(e4); l_dirs.push_back(1);
+ l_edges.push_back(_e[2]); l_dirs.push_back(1);
+ l_edges.push_back(_e[3]); l_dirs.push_back(1);
}
}
fourierFace::~fourierFace()
{
+ if(!_discrete){
+ // this face was just used temporarily for meshing, so don't
+ // delete the mesh vertices or the mesh elements: they have been
+ // transferred elsewhere
+ for(int i = 0; i < 4; i++){
+ if(_e[i]){
+ _e[i]->mesh_vertices.clear();
+ delete _e[i];
+ }
+ }
+ for(int i = 0; i < 4; i++){
+ if(_v[i]){
+ _v[i]->mesh_vertices.clear();
+ delete _v[i];
+ }
+ }
+ triangles.clear();
+ quadrangles.clear();
+ mesh_vertices.clear();
+ l_edges.clear();
+ }
}
Range<double> fourierFace::parBounds(int i) const
diff --git a/Geo/fourierModel.h b/Geo/fourierModel.h
index 6111ae8a71d28489f1985d83416636eb236d901f..10b3e072f87cb5732ef452c35a292d04550b7f89 100644
--- a/Geo/fourierModel.h
+++ b/Geo/fourierModel.h
@@ -33,8 +33,6 @@ class fourierVertex : public GVertex {
};
class fourierEdge : public GEdge {
- private:
- GVertex *v0, *v1;
public:
fourierEdge(GModel *m, int num, GVertex *v1, GVertex *v2);
virtual ~fourierEdge() {}
@@ -57,8 +55,11 @@ class fourierEdge : public GEdge {
class fourierFace : public GFace {
private:
- int _num;
+ // flag to know if is the face is already meshed
int _discrete;
+ // vertices and edges purely local to the face (not shared with the model)
+ GVertex *_v[4];
+ GEdge *_e[4];
public:
fourierFace(GModel *m, int num);
fourierFace(GFace *f, std::vector<MVertex*> &loop, std::vector<MVertex*> &hole);
diff --git a/Geo/gmshModel.cpp b/Geo/gmshModel.cpp
index b2288ea93fc6e82544ffa87a6c43f2e591be2209..7c78a0ddce03002afa8d078b6f547ffe37bcb41f 100644
--- a/Geo/gmshModel.cpp
+++ b/Geo/gmshModel.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshModel.cpp,v 1.17 2006-09-07 05:04:38 geuzaine Exp $
+// $Id: gmshModel.cpp,v 1.18 2006-09-10 15:34:12 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -23,6 +23,7 @@
#include "Mesh.h"
#include "Geo.h"
#include "Tools.h"
+#include "Numeric.h"
#include "Message.h"
#include "gmshVertex.h"
#include "gmshFace.h"
@@ -104,14 +105,14 @@ void gmshModel::import()
List_Read(p->Entities, j, &num);
GEntity *ge = 0;
switch(p->Typ){
- case MSH_PHYSICAL_POINT: ge = vertexByTag(num); break;
- case MSH_PHYSICAL_LINE: ge = edgeByTag(num); break;
- case MSH_PHYSICAL_SURFACE: ge = faceByTag(num); break;
- case MSH_PHYSICAL_VOLUME: ge = regionByTag(num); break;
+ case MSH_PHYSICAL_POINT: ge = vertexByTag(abs(num)); break;
+ case MSH_PHYSICAL_LINE: ge = edgeByTag(abs(num)); break;
+ case MSH_PHYSICAL_SURFACE: ge = faceByTag(abs(num)); break;
+ case MSH_PHYSICAL_VOLUME: ge = regionByTag(abs(num)); break;
}
if(ge && std::find(ge->physicals.begin(), ge->physicals.end(), p->Num) ==
ge->physicals.end())
- ge->physicals.push_back(p->Num);
+ ge->physicals.push_back(sign(num) * p->Num);
}
}
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index bdaf872207ea6e52e526fbb6e8c6884005ef66b7..913f8bd012482bad6deb6e4d7de34e922a55a948 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -31,7 +31,6 @@
#define MIN(a,b) (((a)<(b))?(a):(b))
#define MAX(a,b) (((a)<(b))?(b):(a))
#define SQR(a) ((a)*(a))
-#define SIGN(a,b)((b) >= 0.0 ? fabs(a) : -fabs(a))
#define IMIN MIN
#define LMIN MIN
diff --git a/Numeric/gsl_brent.cpp b/Numeric/gsl_brent.cpp
index e1aa2959119114e0f7f9722294e5698935346d2e..15272980608293167a0ac973c5089233ca939ce3 100644
--- a/Numeric/gsl_brent.cpp
+++ b/Numeric/gsl_brent.cpp
@@ -1,4 +1,4 @@
-// $Id: gsl_brent.cpp,v 1.13 2006-01-06 00:34:27 geuzaine Exp $
+// $Id: gsl_brent.cpp,v 1.14 2006-09-10 15:34:12 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -109,6 +109,7 @@ double brent(double ax, double bx, double cx,
#define MYGOLD_ 1.618034
#define MYLIMIT_ 100.0
#define MYTINY_ 1.0e-20
+#define SIGN(a,b)((b) >= 0.0 ? fabs(a) : -fabs(a))
void mnbrak(double *ax, double *bx, double *cx,
double *fa_dummy, double *fb_dummy, double *fc_dummy,