From 93596920d8cb83c70d9c45e0916138363ddc1c08 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 20 Aug 2006 03:44:15 +0000
Subject: [PATCH] small cleanups + start rewriting 2nd order stuff
---
Geo/MEdge.h | 4 +-
Geo/MElement.h | 18 +-
Geo/MVertex.h | 5 +
Graphics/Mesh.cpp | 4 +-
Mesh/Mesh.h | 2 -
Mesh/SecondOrder.cpp | 600 ++++++++++---------------------------------
6 files changed, 147 insertions(+), 486 deletions(-)
diff --git a/Geo/MEdge.h b/Geo/MEdge.h
index a704c106b2..c5e74bc3ef 100644
--- a/Geo/MEdge.h
+++ b/Geo/MEdge.h
@@ -23,14 +23,12 @@
#include "MVertex.h"
#include "SVector3.h"
-class MElement;
-
class MEdge {
private:
MVertex *_v[2];
public:
- MEdge(MVertex *v0, MVertex *v1, MElement *e)
+ MEdge(MVertex *v0, MVertex *v1)
{
_v[0] = v0; _v[1] = v1;
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index f26241218c..8565992b7e 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -154,7 +154,7 @@ class MLine : public MElement {
inline int getNumVertices(){ return 2; }
inline MVertex *getVertex(int num){ return _v[num]; }
virtual int getNumEdges(){ return 1; }
- virtual MEdge getEdge(int num){ return MEdge(_v[0], _v[1], this); }
+ virtual MEdge getEdge(int num){ return MEdge(_v[0], _v[1]); }
virtual int getNumFaces(){ return 0; }
virtual MFace getFace(int num){ throw; }
int getTypeForMSH(){ return LGN1; }
@@ -184,7 +184,7 @@ class MLine2 : public MLine {
};
int i0 = edges_lin2[num][0];
int i1 = edges_lin2[num][1];
- return MEdge(i0 < 2 ? _v[i0] : _vs[i0 - 2], i1 < 2 ? _v[i1] : _vs[i1 - 2], this);
+ return MEdge(i0 < 2 ? _v[i0] : _vs[i0 - 2], i1 < 2 ? _v[i1] : _vs[i1 - 2]);
}
int getTypeForMSH(){ return LGN2; }
int getTypeForUNV(){ return 24; } // BEAM2
@@ -207,7 +207,7 @@ class MTriangle : public MElement {
virtual int getNumEdges(){ return 3; }
virtual MEdge getEdge(int num)
{
- return MEdge(_v[edges_tetra[num][0]], _v[edges_tetra[num][1]], this);
+ return MEdge(_v[edges_tetra[num][0]], _v[edges_tetra[num][1]]);
}
virtual int getNumFaces(){ return 1; }
virtual MFace getFace(int num)
@@ -244,7 +244,7 @@ class MTriangle2 : public MTriangle {
};
int i0 = edges_tri2[num][0];
int i1 = edges_tri2[num][1];
- return MEdge(i0 < 3 ? _v[i0] : _vs[i0 - 3], i1 < 3 ? _v[i1] : _vs[i1 - 3], this);
+ return MEdge(i0 < 3 ? _v[i0] : _vs[i0 - 3], i1 < 3 ? _v[i1] : _vs[i1 - 3]);
}
int getNumFacesRep(){ return 4; }
MFace getFaceRep(int num)
@@ -283,7 +283,7 @@ class MQuadrangle : public MElement {
virtual int getNumEdges(){ return 4; }
virtual MEdge getEdge(int num)
{
- return MEdge(_v[edges_quad[num][0]], _v[edges_quad[num][1]], this);
+ return MEdge(_v[edges_quad[num][0]], _v[edges_quad[num][1]]);
}
virtual int getNumFaces(){ return 1; }
virtual MFace getFace(int num){ return MFace(_v[0], _v[1], _v[2], _v[3]); }
@@ -330,7 +330,7 @@ class MTetrahedron : public MElement {
virtual int getNumEdges(){ return 6; }
virtual MEdge getEdge(int num)
{
- return MEdge(_v[edges_tetra[num][0]], _v[edges_tetra[num][1]], this);
+ return MEdge(_v[edges_tetra[num][0]], _v[edges_tetra[num][1]]);
}
virtual int getNumFaces(){ return 4; }
virtual MFace getFace(int num)
@@ -417,7 +417,7 @@ class MHexahedron : public MElement {
virtual int getNumEdges(){ return 12; }
virtual MEdge getEdge(int num)
{
- return MEdge(_v[edges_hexa[num][0]], _v[edges_hexa[num][1]], this);
+ return MEdge(_v[edges_hexa[num][0]], _v[edges_hexa[num][1]]);
}
virtual int getNumFaces(){ return 6; }
virtual MFace getFace(int num)
@@ -516,7 +516,7 @@ class MPrism : public MElement {
virtual int getNumEdges(){ return 9; }
virtual MEdge getEdge(int num)
{
- return MEdge(_v[edges_prism[num][0]], _v[edges_prism[num][1]], this);
+ return MEdge(_v[edges_prism[num][0]], _v[edges_prism[num][1]]);
}
virtual int getNumFaces(){ return 5; }
virtual MFace getFace(int num)
@@ -612,7 +612,7 @@ class MPyramid : public MElement {
virtual int getNumEdges(){ return 8; }
virtual MEdge getEdge(int num)
{
- return MEdge(_v[edges_pyramid[num][0]], _v[edges_pyramid[num][1]], this);
+ return MEdge(_v[edges_pyramid[num][0]], _v[edges_pyramid[num][1]]);
}
virtual int getNumFaces(){ return 5; }
virtual MFace getFace(int num)
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 788bbef907..86e6584c40 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -72,6 +72,9 @@ class MVertex{
inline int getNum() const {return _num;}
inline void setNum(int num) { _num = num; }
+ // Get ith parameter
+ bool getParameter(int i, double &par){ return false; }
+
// IO routines
void writeMSH(FILE *fp, bool binary=false, double scalingFactor=1.0);
void writeMSH(FILE *fp, double version, bool binary, int num,
@@ -90,6 +93,7 @@ class MEdgeVertex : public MVertex{
{
}
~MEdgeVertex(){}
+ bool getParameter(int i, double &par){ par = _u; return true; }
};
class MFaceVertex : public MVertex{
@@ -101,6 +105,7 @@ class MFaceVertex : public MVertex{
{
}
~MFaceVertex(){}
+ bool getParameter(int i, double &par){ par = i ? _v : _u; return true; }
};
class MVertexLessThanLexicographic{
diff --git a/Graphics/Mesh.cpp b/Graphics/Mesh.cpp
index dbcc15c5ff..68b7a0d5fd 100644
--- a/Graphics/Mesh.cpp
+++ b/Graphics/Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Mesh.cpp,v 1.180 2006-08-19 18:48:06 geuzaine Exp $
+// $Id: Mesh.cpp,v 1.181 2006-08-20 03:44:15 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -321,7 +321,7 @@ static void drawBarycentricDual(std::vector<T*> &elements)
glVertex3d(p.x(), p.y(), p.z());
for(int k = 0; k < f.getNumVertices(); k++){
MEdge e(f.getVertex(k), (k == f.getNumVertices() - 1) ?
- f.getVertex(0) : f.getVertex(k + 1), ele);
+ f.getVertex(0) : f.getVertex(k + 1));
SPoint3 pe = e.barycenter();
glVertex3d(p.x(), p.y(), p.z());
glVertex3d(pe.x(), pe.y(), pe.z());
diff --git a/Mesh/Mesh.h b/Mesh/Mesh.h
index 1b3743fc75..2af1cbb557 100644
--- a/Mesh/Mesh.h
+++ b/Mesh/Mesh.h
@@ -398,7 +398,5 @@ void ApplyLcFactor();
void Degre1();
void Degre2(int dim);
-void Degre2_Curve(void *a, void *b);
-void Degre2_Surface(void *a, void *b);
#endif
diff --git a/Mesh/SecondOrder.cpp b/Mesh/SecondOrder.cpp
index 33333e7d1c..d29d374d59 100644
--- a/Mesh/SecondOrder.cpp
+++ b/Mesh/SecondOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: SecondOrder.cpp,v 1.39 2006-08-05 10:05:45 geuzaine Exp $
+// $Id: SecondOrder.cpp,v 1.40 2006-08-20 03:44:15 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -19,510 +19,152 @@
//
// Please report all bugs and problems to <gmsh@geuz.org>.
-#include "Gmsh.h"
-#include "Geo.h"
-#include "Mesh.h"
-#include "Utils.h"
-#include "Interpolation.h"
-#include "Numeric.h"
+#include "GModel.h"
+#include "GFace.h"
+#include "GEdge.h"
+#include "MElement.h"
+#include "Message.h"
#include "OS.h"
-extern Mesh *THEM;
+extern GModel *GMODEL;
-extern int edges_tetra[6][2];
-extern int edges_quad[4][2];
-extern int edges_hexa[12][2];
-extern int edges_prism[9][2];
-extern int edges_pyramid[8][2];
+// Creation of middle-edge vertices
-extern int quadfaces_hexa[6][4];
-extern int quadfaces_prism[3][4];
-extern int quadfaces_pyramid[1][4];
-
-static Surface *THES = NULL;
-static Curve *THEC = NULL;
-static EdgesContainer *edges = NULL;
-static QuadFacesContainer *faces = NULL;
-static List_T *VerticesToDelete = NULL;
-
-// using the following onCurve and onSurface functions instead of
-// adding the second order nodes directly in the mesh algorithms is
-// far less efficient... But it's much more general, as we want to be
-// able to couple many algorithms, some of which we don't want
-// to/cannot modify.
-
-Vertex *onCurve(Vertex * v1, Vertex * v2)
+MVertex *addEdgeVertex(GEdge *e, std::pair<MVertex*, MVertex*> &p, bool linear)
{
- Vertex v, *pv;
-
- if(!THEC)
- return NULL;
+ MVertex *v0(p.first), *v1(p.second), *v;
+ MEdge edge(v0, v1);
- int ok1 = 1, ok2 = 1;
- double u1 = 0., u2 = 0.;
-
- if(!THEC->beg || !THEC->end){
- ok1 = ok2 = 0;
+ if(linear){
+ SPoint3 pc = edge.barycenter();
+ v = new MVertex(pc.x(), pc.y(), pc.z(), e);
}
else{
- if(List_Nbr(v1->ListCurves) == 1){
- u1 = v1->u;
- }
- else if(v1->Num == THEC->beg->Num){
- u1 = THEC->ubeg;
- }
- else if(v1->Num == THEC->end->Num){
- u1 = THEC->uend;
+ double u0 = 1e6, u1 = 1e6;
+ Range<double> bounds = e->parBounds(0);
+ if(e->getBeginVertex()->mesh_vertices.size() &&
+ e->getBeginVertex()->mesh_vertices[0] == v0)
+ u0 = bounds.low();
+ else if(e->getEndVertex()->mesh_vertices.size() &&
+ e->getEndVertex()->mesh_vertices[0] == v0)
+ u0 = bounds.high();
+ else
+ v0->getParameter(0, u0);
+ if(e->getBeginVertex()->mesh_vertices.size() &&
+ e->getBeginVertex()->mesh_vertices[0] == v1)
+ u1 = bounds.low();
+ else if(e->getEndVertex()->mesh_vertices.size() &&
+ e->getEndVertex()->mesh_vertices[0] == v1)
+ u1 = bounds.high();
+ else
+ v1->getParameter(0, u1);
+
+ printf("u12 = %g %g\n", u0, u1);
+ if(u0 < 1e6 && u1 < 1e6){
+ double uc = 0.5 * (u0 + u1);
+ GPoint pc = e->point(uc);
+ v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), e, uc);
}
else{
- ok1 = 0;
- }
-
- if(List_Nbr(v2->ListCurves) == 1){
- u2 = v2->u;
- }
- else if(v2->Num == THEC->beg->Num){
- u2 = THEC->ubeg;
- }
- else if(v2->Num == THEC->end->Num){
- u2 = THEC->uend;
- }
- else{
- ok2 = 0;
- }
-
- // Ugly fix for closed curves
- if((THEC->beg->Num == THEC->end->Num) &&
- (u1 == THEC->ubeg || u1 == THEC->uend ||
- u2 == THEC->ubeg || u2 == THEC->uend)){
- ok1 = ok2 = 0;
+ SPoint3 pc = edge.barycenter();
+ v = new MVertex(pc.x(), pc.y(), pc.z(), e);
}
}
+ e->mesh_vertices.push_back(v);
+ return v;
+}
- if(ok1 && ok2){
- v = InterpolateCurve(THEC, 0.5 * (u1 + u2), 0);
+MVertex *addEdgeVertex(GFace *f, std::pair<MVertex* , MVertex*> &p, bool linear)
+{
+ MVertex *v0(p.first), *v1(p.second), *v;
+ MEdge edge(v0, v1);
+
+ if(linear){
+ SPoint3 pc = edge.barycenter();
+ v = new MVertex(pc.x(), pc.y(), pc.z(), f);
}
else{
- v.Pos.X = (v1->Pos.X + v2->Pos.X) * 0.5;
- v.Pos.Y = (v1->Pos.Y + v2->Pos.Y) * 0.5;
- v.Pos.Z = (v1->Pos.Z + v2->Pos.Z) * 0.5;
- v.lc = (v1->lc + v2->lc) * 0.5;
- }
-
- pv = Create_Vertex(++THEM->MaxPointNum, v.Pos.X, v.Pos.Y, v.Pos.Z, v.lc, v.u);
- pv->Degree = 2;
-
- if(!pv->ListCurves) {
- pv->ListCurves = List_Create(1, 1, sizeof(Curve *));
+ SPoint2 p1 = f->parFromPoint(SPoint3(v0->x(), v0->y(), v0->z()));
+ SPoint2 p2 = f->parFromPoint(SPoint3(v1->x(), v1->y(), v1->z()));
+ double uc = 0.5 * (p1[0] + p2[0]);
+ double vc = 0.5 * (p1[1] + p2[1]);
+ GPoint pc = f->point(uc, vc);
+ v = new MFaceVertex(pc.x(), pc.y(), pc.z(), f, uc, vc);
}
- List_Add(pv->ListCurves, &THEC);
- return pv;
-}
-
-Vertex *onSurface(Vertex * v1, Vertex * v2)
-{
- Vertex v, *pv;
- double U, V, U1, U2, V1, V2;
- if(!THES)
- return NULL;
- if(THES->Typ == MSH_SURF_PLAN || THES->Typ == MSH_SURF_DISCRETE)
- return NULL;
-
- XYZtoUV(THES, v1->Pos.X, v1->Pos.Y, v1->Pos.Z, &U1, &V1, 1.0);
- XYZtoUV(THES, v2->Pos.X, v2->Pos.Y, v2->Pos.Z, &U2, &V2, 1.0);
-
- U = 0.5 * (U1 + U2);
- V = 0.5 * (V1 + V2);
- v = InterpolateSurface(THES, U, V, 0, 0);
- pv = Create_Vertex(++THEM->MaxPointNum, v.Pos.X, v.Pos.Y, v.Pos.Z, v.lc, v.u);
- pv->Degree = 2;
-
- return pv;
+ f->mesh_vertices.push_back(v);
+ return v;
}
-Vertex *onSurface(Vertex * v1, Vertex * v2, Vertex * v3, Vertex * v4)
+MVertex *addEdgeVertex(GRegion *r, std::pair<MVertex* , MVertex*> &p, bool linear)
{
- Vertex v, *pv;
- double U, V, U1, U2, U3, U4, V1, V2, V3, V4;
-
- if(!THES)
- return NULL;
- if(THES->Typ == MSH_SURF_PLAN || THES->Typ == MSH_SURF_DISCRETE)
- return NULL;
-
- XYZtoUV(THES, v1->Pos.X, v1->Pos.Y, v1->Pos.Z, &U1, &V1, 1.0);
- XYZtoUV(THES, v2->Pos.X, v2->Pos.Y, v2->Pos.Z, &U2, &V2, 1.0);
- XYZtoUV(THES, v3->Pos.X, v3->Pos.Y, v3->Pos.Z, &U3, &V3, 1.0);
- XYZtoUV(THES, v4->Pos.X, v4->Pos.Y, v4->Pos.Z, &U4, &V4, 1.0);
-
- U = 0.25 * (U1 + U2 + U3 + U4);
- V = 0.25 * (V1 + V2 + V3 + V4);
- v = InterpolateSurface(THES, U, V, 0, 0);
- pv = Create_Vertex(++THEM->MaxPointNum, v.Pos.X, v.Pos.Y, v.Pos.Z, v.lc, v.u);
- pv->Degree = 2;
-
- return pv;
+ MVertex *v0(p.first), *v1(p.second), *v;
+ MEdge edge(v0, v1);
+ SPoint3 pc = edge.barycenter();
+ v = new MVertex(pc.x(), pc.y(), pc.z(), r);
+ r->mesh_vertices.push_back(v);
+ return v;
}
-void putMiddleEdgePoint(void *a, void *b)
+template<class T, class U>
+void addEdgeVertices(U *ge, std::vector<T*> &elements, bool linear,
+ std::map<std::pair<MVertex*,MVertex*>, MVertex* > &edges)
{
- Vertex *v;
- int N;
-
- Edge *ed = (Edge *) a;
-
- if(ed->newv){
- v = ed->newv;
- }
- else if((v = onCurve(ed->V[0], ed->V[1]))){
- }
- else if((v = onSurface(ed->V[0], ed->V[1]))){
- }
- else{
- v = Create_Vertex(++THEM->MaxPointNum,
- 0.5 * (ed->V[0]->Pos.X + ed->V[1]->Pos.X),
- 0.5 * (ed->V[0]->Pos.Y + ed->V[1]->Pos.Y),
- 0.5 * (ed->V[0]->Pos.Z + ed->V[1]->Pos.Z),
- 0.5 * (ed->V[0]->lc + ed->V[1]->lc),
- 0.5 * (ed->V[0]->u + ed->V[1]->u));
- v->Degree = 2;
- }
-
- ed->newv = v;
- Tree_Insert(THEM->Vertices, &v);
-
- for(int i = 0; i < List_Nbr(ed->Simplexes); i++) {
- Simplex *s;
- List_Read(ed->Simplexes, i, &s);
- if(s->V[3]) // tetrahedron
- N = 6;
- else if(s->V[2]) // triangle
- N = 3;
- else // line
- N = 1;
- if(!s->VSUP)
- s->VSUP = (Vertex **) Malloc(N * sizeof(Vertex *));
- for(int j = 0; j < N; j++) {
- if((!compareVertex(&s->V[edges_tetra[j][0]], &ed->V[0]) &&
- !compareVertex(&s->V[edges_tetra[j][1]], &ed->V[1])) ||
- (!compareVertex(&s->V[edges_tetra[j][0]], &ed->V[1]) &&
- !compareVertex(&s->V[edges_tetra[j][1]], &ed->V[0]))) {
- s->VSUP[j] = v;
- }
+ for(unsigned int i = 0; i < elements.size(); i++){
+ for(int j = 0; j < elements[i]->getNumEdges(); j++){
+ MEdge e = elements[i]->getEdge(j);
+ std::pair<MVertex*, MVertex*> p(e.getMinVertex(), e.getMaxVertex());
+ if(!edges.count(p)) edges[p] = addEdgeVertex(ge, p, linear);
}
}
-
- for(int i = 0; i < List_Nbr(ed->Quadrangles); i++) {
- Quadrangle *q;
- List_Read(ed->Quadrangles, i, &q);
- if(!q->VSUP)
- q->VSUP = (Vertex **) Malloc((4 + 1) * sizeof(Vertex *));
- for(int j = 0; j < 4; j++) {
- if((!compareVertex(&q->V[edges_quad[j][0]], &ed->V[0]) &&
- !compareVertex(&q->V[edges_quad[j][1]], &ed->V[1])) ||
- (!compareVertex(&q->V[edges_quad[j][0]], &ed->V[1]) &&
- !compareVertex(&q->V[edges_quad[j][1]], &ed->V[0]))) {
- q->VSUP[j] = v;
- }
- }
- }
-
- for(int i = 0; i < List_Nbr(ed->Hexahedra); i++) {
- Hexahedron *h;
- List_Read(ed->Hexahedra, i, &h);
- if(!h->VSUP)
- h->VSUP = (Vertex **) Malloc((12 + 6 + 1) * sizeof(Vertex *));
- for(int j = 0; j < 12; j++) {
- if((!compareVertex(&h->V[edges_hexa[j][0]], &ed->V[0]) &&
- !compareVertex(&h->V[edges_hexa[j][1]], &ed->V[1])) ||
- (!compareVertex(&h->V[edges_hexa[j][0]], &ed->V[1]) &&
- !compareVertex(&h->V[edges_hexa[j][1]], &ed->V[0]))) {
- h->VSUP[j] = v;
- }
- }
- }
-
- for(int i = 0; i < List_Nbr(ed->Prisms); i++) {
- Prism *p;
- List_Read(ed->Prisms, i, &p);
- if(!p->VSUP)
- p->VSUP = (Vertex **) Malloc((9 + 3) * sizeof(Vertex *));
- for(int j = 0; j < 9; j++) {
- if((!compareVertex(&p->V[edges_prism[j][0]], &ed->V[0]) &&
- !compareVertex(&p->V[edges_prism[j][1]], &ed->V[1])) ||
- (!compareVertex(&p->V[edges_prism[j][0]], &ed->V[1]) &&
- !compareVertex(&p->V[edges_prism[j][1]], &ed->V[0]))) {
- p->VSUP[j] = v;
- }
- }
- }
-
- for(int i = 0; i < List_Nbr(ed->Pyramids); i++) {
- Pyramid *p;
- List_Read(ed->Pyramids, i, &p);
- if(!p->VSUP)
- p->VSUP = (Vertex **) Malloc((8 + 1) * sizeof(Vertex *));
- for(int j = 0; j < 8; j++) {
- if((!compareVertex(&p->V[edges_pyramid[j][0]], &ed->V[0]) &&
- !compareVertex(&p->V[edges_pyramid[j][1]], &ed->V[1])) ||
- (!compareVertex(&p->V[edges_pyramid[j][0]], &ed->V[1]) &&
- !compareVertex(&p->V[edges_pyramid[j][1]], &ed->V[0]))) {
- p->VSUP[j] = v;
- }
- }
- }
-
}
-void putMiddleFacePoint(void *a, void *b)
-{
- QuadFace *fac = (QuadFace*)a;
- Vertex *v;
-
- if(fac->newv){
- v = fac->newv;
- }
- else if((v = onSurface(fac->V[0], fac->V[1], fac->V[2], fac->V[3]))){
- }
- else{
- v = Create_Vertex(++THEM->MaxPointNum,
- 0.25 * (fac->V[0]->Pos.X + fac->V[1]->Pos.X +
- fac->V[2]->Pos.X + fac->V[3]->Pos.X),
- 0.25 * (fac->V[0]->Pos.Y + fac->V[1]->Pos.Y +
- fac->V[2]->Pos.Y + fac->V[3]->Pos.Y),
- 0.25 * (fac->V[0]->Pos.Z + fac->V[1]->Pos.Z +
- fac->V[2]->Pos.Z + fac->V[3]->Pos.Z),
- 0.25 * (fac->V[0]->lc + fac->V[1]->lc +
- fac->V[2]->lc + fac->V[3]->lc),
- 0.25 * (fac->V[0]->u + fac->V[1]->u +
- fac->V[2]->u + fac->V[3]->u));
- v->Degree = 2;
- }
-
- fac->newv = v;
- Tree_Insert(THEM->Vertices, &v);
+// Creation of middle-face vertices
- if(fac->quad && fac->quad->VSUP){
- fac->quad->VSUP[4] = v;
- }
-
- QuadFace F;
-
- for(int i = 0; i < 2; i++) {
- if(fac->hexa[i] && fac->hexa[i]->VSUP){
- for(int j = 0; j < 6; j++) {
- F.V[0] = fac->hexa[i]->V[quadfaces_hexa[j][0]];
- F.V[1] = fac->hexa[i]->V[quadfaces_hexa[j][1]];
- F.V[2] = fac->hexa[i]->V[quadfaces_hexa[j][2]];
- F.V[3] = fac->hexa[i]->V[quadfaces_hexa[j][3]];
- qsort(F.V, 4, sizeof(Vertex *), compareVertex);
- if(!compareQuadFace(fac, &F))
- fac->hexa[i]->VSUP[12 + j] = v;
- }
- }
- }
-
- for(int i = 0; i < 2; i++) {
- if(fac->prism[i] && fac->prism[i]->VSUP){
- for(int j = 0; j < 3; j++) {
- F.V[0] = fac->prism[i]->V[quadfaces_prism[j][0]];
- F.V[1] = fac->prism[i]->V[quadfaces_prism[j][1]];
- F.V[2] = fac->prism[i]->V[quadfaces_prism[j][2]];
- F.V[3] = fac->prism[i]->V[quadfaces_prism[j][3]];
- qsort(F.V, 4, sizeof(Vertex *), compareVertex);
- if(!compareQuadFace(fac, &F))
- fac->prism[i]->VSUP[9 + j] = v;
- }
- }
- }
-
- for(int i = 0; i < 2; i++) {
- if(fac->pyramid[i] && fac->pyramid[i]->VSUP){
- F.V[0] = fac->pyramid[i]->V[quadfaces_pyramid[0][0]];
- F.V[1] = fac->pyramid[i]->V[quadfaces_pyramid[0][1]];
- F.V[2] = fac->pyramid[i]->V[quadfaces_pyramid[0][2]];
- F.V[3] = fac->pyramid[i]->V[quadfaces_pyramid[0][3]];
- qsort(F.V, 4, sizeof(Vertex *), compareVertex);
- if(!compareQuadFace(fac, &F))
- fac->pyramid[i]->VSUP[8] = v;
- }
- }
-}
-
-void putMiddleVolumePoint(void *a, void *b)
+MVertex *addFaceVertex(GFace *f, std::vector<MVertex*> &p, bool linear)
{
- Hexahedron *h = *(Hexahedron**)a;
- double x = 0., y = 0., z = 0., lc = 0., u = 0.;
- for(int i = 0; i < 8; i++){
- x += h->V[i]->Pos.X;
- y += h->V[i]->Pos.Y;
- z += h->V[i]->Pos.Z;
- lc += h->V[i]->lc;
- u += h->V[i]->u;
+ /*
+ MVertex *v;
+ if(linear || e->dim() == 3){
+ // just interpolate linear between the 4
}
- for(int i = 0; i < 18; i++){
- x += h->VSUP[i]->Pos.X;
- y += h->VSUP[i]->Pos.Y;
- z += h->VSUP[i]->Pos.Z;
- lc += h->VSUP[i]->lc;
- u += h->VSUP[i]->u;
+ else if(e->dim() == 2){
+ xyz2uv();
+ xyz2uv();
+ xyz2uv();
+ xyz2uv();
+ // interpolate
}
- x /= 26.;
- y /= 26.;
- z /= 26.;
- lc /= 26.;
- u /= 26.;
- Vertex *v = Create_Vertex(++THEM->MaxPointNum, x, y, z, lc, u);
- v->Degree = 2;
- Tree_Insert(THEM->Vertices, &v);
- h->VSUP[18] = v;
+ e->mesh_vertices.push_back(v);
+ */
+ return 0;
}
-void ResetDegre2_Vertex(void *a, void *b)
+template<class T>
+void addFaceVertices(GEntity *e, std::vector<T*> &elements)
{
- Vertex *v = *(Vertex**)a;
- if(v->Degree == 2)
- List_Add(VerticesToDelete, &v);
}
-void ResetDegre2_Simplex(void *a, void *b)
-{
- Simplex *s = *(Simplex**)a;
- Free(s->VSUP);
- s->VSUP = NULL;
-}
+// Creation of middle-volume vertices
-void ResetDegre2_Quadrangle(void *a, void *b)
+MVertex *addVolumeVertex(GEntity *ge, MElement *elem)
{
- Quadrangle *q = *(Quadrangle**)a;
- Free(q->VSUP);
- q->VSUP = NULL;
+ SPoint3 pc = elem->barycenter();
+ MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), ge);
+ ge->mesh_vertices.push_back(v);
+ return v;
}
-void ResetDegre2_Hexahedron(void *a, void *b)
-{
- Hexahedron *h = *(Hexahedron**)a;
- Free(h->VSUP);
- h->VSUP = NULL;
-}
-
-void ResetDegre2_Prism(void *a, void *b)
-{
- Prism *p = *(Prism**)a;
- Free(p->VSUP);
- p->VSUP = NULL;
-}
-
-void ResetDegre2_Pyramid(void *a, void *b)
-{
- Pyramid *p = *(Pyramid**)a;
- Free(p->VSUP);
- p->VSUP = NULL;
-}
-
-void ResetDegre2_Curve(void *a, void *b)
-{
- Curve *c = *(Curve**)a;
- Tree_Action(c->Simplexes, ResetDegre2_Simplex);
-}
-
-void ResetDegre2_Surface(void *a, void *b)
-{
- Surface *s = *(Surface**)a;
- Tree_Action(s->Simplexes, ResetDegre2_Simplex);
- Tree_Action(s->Quadrangles, ResetDegre2_Quadrangle);
-}
-
-void ResetDegre2_Volume(void *a, void *b)
-{
- Volume *v = *(Volume**)a;
- Tree_Action(v->Simplexes, ResetDegre2_Simplex);
- Tree_Action(v->Hexahedra, ResetDegre2_Hexahedron);
- Tree_Action(v->Prisms, ResetDegre2_Prism);
- Tree_Action(v->Pyramids, ResetDegre2_Pyramid);
-}
+// Main routines
void Degre1()
{
- // transform any SimplexBase into "real" simplices
- Move_SimplexBaseToSimplex(3);
-
- // (re-)initialize the global tree of edges/quadfaces
- if(edges)
- delete edges;
- edges = new EdgesContainer();
-
- if(faces)
- delete faces;
- faces = new QuadFacesContainer();
-
- // reset VSUP in each element
- Tree_Action(THEM->Curves, ResetDegre2_Curve);
- Tree_Action(THEM->Surfaces, ResetDegre2_Surface);
- Tree_Action(THEM->Volumes, ResetDegre2_Volume);
-
- // remove any supp vertex from the database
- if(VerticesToDelete)
- List_Delete(VerticesToDelete);
- VerticesToDelete = List_Create(100, 100, sizeof(Vertex*));
- Tree_Action(THEM->Vertices, ResetDegre2_Vertex);
- for(int i = 0; i < List_Nbr(VerticesToDelete); i++){
- Vertex **v = (Vertex**)List_Pointer(VerticesToDelete, i);
- Tree_Suppress(THEM->Vertices, v);
- Free_Vertex(v, NULL);
- }
-}
-
-void Degre2_Curve(void *a, void *b)
-{
- Curve *c = *(Curve**)a;
- if(c->Num < 0) return;
-
- Msg(STATUS2, "Second order curve %d", c->Num);
-
- edges->AddSimplexTree(c->Simplexes);
- THEC = c;
- THES = NULL;
- Tree_Action(edges->AllEdges, putMiddleEdgePoint);
-}
-
-void Degre2_Surface(void *a, void *b)
-{
- Surface *s = *(Surface**)a;
-
- Msg(STATUS2, "Second order surface %d", s->Num);
-
- edges->AddSimplexTree(s->Simplexes);
- edges->AddQuadrangleTree(s->Quadrangles);
- THEC = NULL;
- THES = s;
- Tree_Action(edges->AllEdges, putMiddleEdgePoint);
-
- faces->AddQuadrangleTree(s->Quadrangles);
- Tree_Action(faces->AllFaces, putMiddleFacePoint);
-}
-
-void Degre2_Volume(void *a, void *b)
-{
- Volume *v = *(Volume**)a;
+ // loop on all elements
+ // - get their edge/face/volume vertices mark them with setVisibility(-1);
+ // - generate a first order element for each element
+ // - swap lists at the end
+ // loop on all vertices and delete all vertices marked (-1)
- Msg(STATUS2, "Second order volume %d", v->Num);
-
- edges->AddSimplexTree(v->Simplexes);
- edges->AddHexahedronTree(v->Hexahedra);
- edges->AddPrismTree(v->Prisms);
- edges->AddPyramidTree(v->Pyramids);
- THEC = NULL;
- THES = NULL;
- Tree_Action(edges->AllEdges, putMiddleEdgePoint);
-
- faces->AddHexahedronTree(v->Hexahedra);
- faces->AddPrismTree(v->Prisms);
- faces->AddPyramidTree(v->Pyramids);
- Tree_Action(faces->AllFaces, putMiddleFacePoint);
-
- Tree_Action(v->Hexahedra, putMiddleVolumePoint);
}
void Degre2(int dim)
@@ -530,15 +172,33 @@ void Degre2(int dim)
Msg(STATUS1, "Mesh second order...");
double t1 = Cpu();
- Degre1();
-
- if(dim >= 1)
- Tree_Action(THEM->Curves, Degre2_Curve);
- if(dim >= 2)
- Tree_Action(THEM->Surfaces, Degre2_Surface);
- if(dim >= 3)
- Tree_Action(THEM->Volumes, Degre2_Volume);
-
+ //bool linear = true;
+ bool linear = false;
+ std::map<std::pair<MVertex*,MVertex*>, MVertex* > edges;
+ std::map<std::vector<MVertex*>, MVertex* > quadFaces;
+
+ // loop over all elements and create unique edges and faces and
+ // their associated new vertices
+ for(GModel::eiter it = GMODEL->firstEdge(); it != GMODEL->lastEdge(); ++it){
+ addEdgeVertices(*it, (*it)->lines, linear, edges);
+ }
+ for(GModel::fiter it = GMODEL->firstFace(); it != GMODEL->lastFace(); ++it){
+ addEdgeVertices(*it, (*it)->triangles, linear, edges);
+ addEdgeVertices(*it, (*it)->quadrangles, linear, edges);
+ //addFaceVertices(*it, (*it)->quadrangles, linear);
+ }
+ for(GModel::riter it = GMODEL->firstRegion(); it != GMODEL->lastRegion(); ++it){
+ addEdgeVertices(*it, (*it)->tetrahedra, linear, edges);
+ addEdgeVertices(*it, (*it)->hexahedra, linear, edges);
+ addEdgeVertices(*it, (*it)->prisms, linear, edges);
+ addEdgeVertices(*it, (*it)->pyramids, linear, edges);
+// addFaceVertices(*it, (*it)->hexahedra, linear);
+// addFaceVertices(*it, (*it)->prisms, linear);
+// addFaceVertices(*it, (*it)->pyramids, linear);
+ }
+
+ // loop on all elements again and create one new element from each
+ // old one, querying the edge/face maps to get the extra vertices
double t2 = Cpu();
Msg(STATUS1, "Mesh second order complete (%g s)", t2 - t1);
}
--
GitLab