diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 997b39fd5672e4717fb129d2167c9c3488588931..34bac55ba99ee99fce1230af4d51bface8b8c814 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.445 2006-08-20 14:12:40 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.446 2006-08-20 17:02:27 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -1112,12 +1112,12 @@ void mesh_options_ok_cb(CALLBACK_ARGS)
opt_mesh_constrained_bgmesh(0, GMSH_SET, WID->mesh_butt[5]->value());
opt_mesh_points(0, GMSH_SET, WID->mesh_butt[6]->value());
opt_mesh_lines(0, GMSH_SET, WID->mesh_butt[7]->value());
- opt_mesh_triangles(0, GMSH_SET, WID->mesh_menu_butt->menu()[0].value());
- opt_mesh_quadrangles(0, GMSH_SET, WID->mesh_menu_butt->menu()[1].value());
- opt_mesh_tetrahedra(0, GMSH_SET, WID->mesh_menu_butt->menu()[2].value());
- opt_mesh_hexahedra(0, GMSH_SET, WID->mesh_menu_butt->menu()[3].value());
- opt_mesh_prisms(0, GMSH_SET, WID->mesh_menu_butt->menu()[4].value());
- opt_mesh_pyramids(0, GMSH_SET, WID->mesh_menu_butt->menu()[5].value());
+ opt_mesh_triangles(0, GMSH_SET, WID->mesh_menu_butt->menu()[0].value() ? 1 : 0);
+ opt_mesh_quadrangles(0, GMSH_SET, WID->mesh_menu_butt->menu()[1].value() ? 1 : 0);
+ opt_mesh_tetrahedra(0, GMSH_SET, WID->mesh_menu_butt->menu()[2].value() ? 1 : 0);
+ opt_mesh_hexahedra(0, GMSH_SET, WID->mesh_menu_butt->menu()[3].value() ? 1 : 0);
+ opt_mesh_prisms(0, GMSH_SET, WID->mesh_menu_butt->menu()[4].value() ? 1 : 0);
+ opt_mesh_pyramids(0, GMSH_SET, WID->mesh_menu_butt->menu()[5].value() ? 1 : 0);
opt_mesh_surfaces_edges(0, GMSH_SET, WID->mesh_butt[8]->value());
opt_mesh_surfaces_faces(0, GMSH_SET, WID->mesh_butt[9]->value());
opt_mesh_volumes_edges(0, GMSH_SET, WID->mesh_butt[10]->value());
diff --git a/Geo/MRep.h b/Geo/MRep.h
index da0a64355deb6d742fe9d95d6e3e8853b1f7ecca..625c8ed7c4d6abd25ee19df77851a01cf1ff24b1 100644
--- a/Geo/MRep.h
+++ b/Geo/MRep.h
@@ -54,7 +54,14 @@ class MRep {
public:
MRep() : va_lines(0), va_triangles(0), va_quads(0), allElementsVisible(true) {}
- virtual ~MRep(){}
+ virtual ~MRep(){ destroy(); }
+
+ // destroy everything
+ void destroy(){
+ resetArrays();
+ edges.clear();
+ allElementsVisible = true;
+ }
// generates the edge representation
virtual void generateEdgeRep() = 0;
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 273162b322a3cba8c82c937dc946a9aba990d214..c43f4a73d9a58de872cf8828ae80243591d6df8a 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.123 2006-08-20 14:12:41 geuzaine Exp $
+# $Id: Makefile,v 1.124 2006-08-20 17:02:28 geuzaine Exp $
#
# Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
#
@@ -472,8 +472,10 @@ SecondOrder.o: SecondOrder.cpp ../Geo/GModel.h ../Geo/GVertex.h \
../Common/Context.h ../DataStr/List.h ../Geo/GFace.h ../Geo/GPoint.h \
../Geo/GEntity.h ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h \
../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/MElement.h \
- ../Geo/SBoundingBox3d.h ../Common/SmoothNormals.h ../Common/Message.h \
- ../Common/OS.h
+ ../Geo/SBoundingBox3d.h ../Common/SmoothNormals.h ../Geo/MRep.h \
+ ../Geo/GEdge.h ../Geo/GFace.h ../Geo/GRegion.h ../Geo/MVertex.h \
+ ../Geo/MEdge.h ../Geo/MElement.h ../Common/VertexArray.h \
+ ../Common/Message.h ../Common/OS.h
# 1 "/Users/geuzaine/.gmsh/Mesh//"
PartitionMesh.o: PartitionMesh.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
diff --git a/Mesh/SecondOrder.cpp b/Mesh/SecondOrder.cpp
index 66e3703fee9a7356c284376a8ae2d319df004144..f737b56fcd6ec5cab69fec2804061851b95956f2 100644
--- a/Mesh/SecondOrder.cpp
+++ b/Mesh/SecondOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: SecondOrder.cpp,v 1.41 2006-08-20 14:12:42 geuzaine Exp $
+// $Id: SecondOrder.cpp,v 1.42 2006-08-20 17:02:28 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -23,139 +23,238 @@
#include "GFace.h"
#include "GEdge.h"
#include "MElement.h"
+#include "MRep.h"
#include "Message.h"
#include "OS.h"
extern GModel *GMODEL;
-// Creation of middle-edge vertices
-
-MVertex *addEdgeVertex(GEdge *e, std::pair<MVertex*, MVertex*> &p, bool linear)
+void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
+ std::map<std::pair<MVertex*,MVertex*>, MVertex* > &edgeVertices,
+ 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(), e);
- }
- else{
- 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);
+ for(int i = 0; i < ele->getNumEdges(); i++){
+ MEdge edge = ele->getEdge(i);
+ std::pair<MVertex*, MVertex*> p(edge.getMinVertex(), edge.getMaxVertex());
+ if(edgeVertices.count(p)){
+ ve.push_back(edgeVertices[p]);
}
else{
- SPoint3 pc = edge.barycenter();
- v = new MVertex(pc.x(), pc.y(), pc.z(), e);
+ MVertex *v, *v0 = edge.getVertex(0), *v1 = edge.getVertex(1);
+ if(linear || ge->geomType() == GEntity::DiscreteCurve){
+ SPoint3 pc = edge.barycenter();
+ v = new MVertex(pc.x(), pc.y(), pc.z(), ge);
+ }
+ else{
+ double u0 = 1e6, u1 = 1e6;
+ Range<double> bounds = ge->parBounds(0);
+ if(ge->getBeginVertex() && ge->getBeginVertex()->mesh_vertices[0] == v0)
+ u0 = bounds.low();
+ else if(ge->getEndVertex() && ge->getEndVertex()->mesh_vertices[0] == v0)
+ u0 = bounds.high();
+ else
+ v0->getParameter(0, u0);
+ if(ge->getBeginVertex() && ge->getBeginVertex()->mesh_vertices[0] == v1)
+ u1 = bounds.low();
+ else if(ge->getEndVertex() && ge->getEndVertex()->mesh_vertices[0] == v1)
+ u1 = bounds.high();
+ else
+ v1->getParameter(0, u1);
+ if(u0 < 1e6 && u1 < 1e6){
+ double uc = 0.5 * (u0 + u1);
+ GPoint pc = ge->point(uc);
+ v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge, uc);
+ }
+ else{
+ // we should normally never end up here
+ SPoint3 pc = edge.barycenter();
+ v = new MVertex(pc.x(), pc.y(), pc.z(), ge);
+ }
+ }
+ edgeVertices[p] = v;
+ ge->mesh_vertices.push_back(v);
+ ve.push_back(v);
}
}
- e->mesh_vertices.push_back(v);
- return v;
-}
+}
-MVertex *addEdgeVertex(GFace *f, std::pair<MVertex* , MVertex*> &p, bool linear)
+void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
+ std::map<std::pair<MVertex*,MVertex*>, MVertex* > &edgeVertices,
+ 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{
- 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);
+ for(int i = 0; i < ele->getNumEdges(); i++){
+ MEdge edge = ele->getEdge(i);
+ std::pair<MVertex*, MVertex*> p(edge.getMinVertex(), edge.getMaxVertex());
+ if(edgeVertices.count(p)){
+ ve.push_back(edgeVertices[p]);
+ }
+ else{
+ MVertex *v, *v0 = edge.getVertex(0), *v1 = edge.getVertex(1);
+ if(linear || gf->geomType() == GEntity::DiscreteSurface){
+ SPoint3 pc = edge.barycenter();
+ v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
+ }
+ else{
+ SPoint2 p0 = gf->parFromPoint(SPoint3(v0->x(), v0->y(), v0->z()));
+ SPoint2 p1 = gf->parFromPoint(SPoint3(v1->x(), v1->y(), v1->z()));
+ double uc = 0.5 * (p0[0] + p1[0]);
+ double vc = 0.5 * (p0[1] + p1[1]);
+ GPoint pc = gf->point(uc, vc);
+ v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, uc, vc);
+ }
+ edgeVertices[p] = v;
+ gf->mesh_vertices.push_back(v);
+ ve.push_back(v);
+ }
}
- f->mesh_vertices.push_back(v);
- return v;
}
-MVertex *addEdgeVertex(GRegion *r, std::pair<MVertex* , MVertex*> &p, bool linear)
+void getEdgeVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &ve,
+ std::map<std::pair<MVertex*,MVertex*>, MVertex* > &edgeVertices,
+ bool linear)
{
- 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;
+ for(int i = 0; i < ele->getNumEdges(); i++){
+ MEdge edge = ele->getEdge(i);
+ std::pair<MVertex*, MVertex*> p(edge.getMinVertex(), edge.getMaxVertex());
+ if(edgeVertices.count(p)){
+ ve.push_back(edgeVertices[p]);
+ }
+ else{
+ SPoint3 pc = edge.barycenter();
+ MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
+ edgeVertices[p] = v;
+ gr->mesh_vertices.push_back(v);
+ ve.push_back(v);
+ }
+ }
}
-template<class T, class U>
-void addEdgeVertices(U *ge, std::vector<T*> &elements, bool linear,
- std::map<std::pair<MVertex*,MVertex*>, MVertex* > &edgeVertices)
+void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
+ std::map<std::vector<MVertex*>, MVertex* > &faceVertices,
+ bool linear)
{
- 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(!edgeVertices.count(p)) edgeVertices[p] = addEdgeVertex(ge, p, linear);
+ vf.clear();
+ for(int i = 0; i < ele->getNumFaces(); i++){
+ MFace face = ele->getFace(i);
+ if(face.getNumVertices() != 4) continue;
+ std::vector<MVertex*> p;
+ face.getOrderedVertices(p);
+ if(faceVertices.count(p)){
+ vf.push_back(faceVertices[p]);
+ }
+ else{
+ MVertex *v;
+ if(linear || gf->geomType() == GEntity::DiscreteSurface){
+ SPoint3 pc = face.barycenter();
+ v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
+ }
+ else{
+ SPoint2 p0 = gf->parFromPoint(SPoint3(p[0]->x(), p[0]->y(), p[0]->z()));
+ SPoint2 p1 = gf->parFromPoint(SPoint3(p[1]->x(), p[1]->y(), p[1]->z()));
+ SPoint2 p2 = gf->parFromPoint(SPoint3(p[2]->x(), p[2]->y(), p[2]->z()));
+ SPoint2 p3 = gf->parFromPoint(SPoint3(p[3]->x(), p[3]->y(), p[3]->z()));
+ double uc = 0.25 * (p0[0] + p1[0] + p2[0] + p3[0]);
+ double vc = 0.25 * (p0[1] + p1[1] + p2[1] + p3[1]);
+ GPoint pc = gf->point(uc, vc);
+ v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, uc, vc);
+ }
+ faceVertices[p] = v;
+ gf->mesh_vertices.push_back(v);
+ vf.push_back(v);
}
}
}
-// Creation of middle-face vertices
-
-MVertex *addFaceVertex(GFace *f, std::vector<MVertex*> &p, bool linear)
+void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &vf,
+ std::map<std::vector<MVertex*>, MVertex* > &faceVertices,
+ bool linear)
{
- //MVertex *v;
- if(linear){
- // just interpolate linear between the 4
+ vf.clear();
+ for(int i = 0; i < ele->getNumFaces(); i++){
+ MFace face = ele->getFace(i);
+ if(face.getNumVertices() != 4) continue;
+ std::vector<MVertex*> p;
+ face.getOrderedVertices(p);
+ if(faceVertices.count(p)){
+ vf.push_back(faceVertices[p]);
+ }
+ else{
+ SPoint3 pc = face.barycenter();
+ MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
+ faceVertices[p] = v;
+ gr->mesh_vertices.push_back(v);
+ vf.push_back(v);
+ }
}
- else{
- //xyz2uv();
- // interpolate
+}
+
+void setSecondOrder(GEdge *ge,
+ std::map<std::pair<MVertex*,MVertex*>, MVertex* > &edgeVertices,
+ bool linear)
+{
+ std::vector<MLine*> l2;
+ for(unsigned int i = 0; i < ge->lines.size(); i++){
+ MLine *l = ge->lines[i];
+ std::vector<MVertex*> ve;
+ getEdgeVertices(ge, l, ve, edgeVertices, linear);
+ l2.push_back(new MLine2(l->getVertex(0), l->getVertex(1), ve[0]));
+ delete l;
}
- //f->mesh_vertices.push_back(v);
- return 0;
+ ge->lines = l2;
+ ge->meshRep->destroy();
}
-MVertex *addFaceVertex(GRegion *r, std::vector<MVertex*> &p, bool linear)
+void setSecondOrder(GFace *gf,
+ std::map<std::pair<MVertex*,MVertex*>, MVertex* > &edgeVertices,
+ std::map<std::vector<MVertex*>, MVertex* > &faceVertices,
+ bool linear)
{
- //MVertex *v;
- // just interpolate linear between the 4
- //r->mesh_vertices.push_back(v);
- return 0;
+ std::vector<MTriangle*> t2;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ MTriangle *t = gf->triangles[i];
+ std::vector<MVertex*> ve;
+ getEdgeVertices(gf, t, ve, edgeVertices, linear);
+ t2.push_back(new MTriangle2(t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ ve[0], ve[1], ve[2]));
+ delete t;
+ }
+ gf->triangles = t2;
+
+ std::vector<MQuadrangle*> q2;
+ for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
+ MQuadrangle *q = gf->quadrangles[i];
+ std::vector<MVertex*> ve, vf;
+ getEdgeVertices(gf, q, ve, edgeVertices, linear);
+ getFaceVertices(gf, q, vf, faceVertices, linear);
+ q2.push_back(new MQuadrangle2(q->getVertex(0), q->getVertex(1), q->getVertex(2),
+ q->getVertex(3), ve[0], ve[1], ve[2], ve[3], vf[0]));
+ delete q;
+ }
+ gf->quadrangles = q2;
+
+ gf->meshRep->destroy();
}
-template<class T, class U>
-void addFaceVertices(U *ge, std::vector<T*> &elements, bool linear,
- std::map<std::vector<MVertex*>, MVertex* > &faceVertices)
+void setSecondOrder(GRegion *gr,
+ std::map<std::pair<MVertex*,MVertex*>, MVertex* > &edgeVertices,
+ std::map<std::vector<MVertex*>, MVertex* > &faceVertices,
+ bool linear)
{
- for(unsigned int i = 0; i < elements.size(); i++){
- for(int j = 0; j < elements[i]->getNumFaces(); j++){
- MFace f = elements[i]->getFace(j);
- if(f.getNumVertices() == 4){
- std::vector<MVertex*> p;
- f.getOrderedVertices(p);
- if(!faceVertices.count(p)) faceVertices[p] = addFaceVertex(ge, p, linear);
- }
- }
+ std::vector<MTetrahedron*> t2;
+ for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
+ MTetrahedron *t = gr->tetrahedra[i];
+ std::vector<MVertex*> ve;
+ getEdgeVertices(gr, t, ve, edgeVertices, linear);
+ t2.push_back(new MTetrahedron2(t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ t->getVertex(3),
+ ve[0], ve[1], ve[2], ve[3], ve[4], ve[5]));
+ delete t;
}
+ gr->tetrahedra = t2;
+
+ // hexa prisms pyr
+
+ gr->meshRep->destroy();
}
// Main routines
@@ -163,11 +262,11 @@ void addFaceVertices(U *ge, std::vector<T*> &elements, bool linear,
void Degre1()
{
// loop on all elements
+ // - if polynomialOrder() == 2
// - 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)
-
}
void Degre2(int dim)
@@ -175,35 +274,18 @@ void Degre2(int dim)
Msg(STATUS1, "Mesh second order...");
double t1 = Cpu();
- //bool linear = true;
- bool linear = false;
+ bool linear = true;
+ //bool linear = false;
std::map<std::pair<MVertex*,MVertex*>, MVertex* > edgeVertices;
std::map<std::vector<MVertex*>, MVertex* > faceVertices;
- // 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, edgeVertices);
- }
- for(GModel::fiter it = GMODEL->firstFace(); it != GMODEL->lastFace(); ++it){
- addEdgeVertices(*it, (*it)->triangles, linear, edgeVertices);
- addEdgeVertices(*it, (*it)->quadrangles, linear, edgeVertices);
- addFaceVertices(*it, (*it)->quadrangles, linear, faceVertices);
- }
- for(GModel::riter it = GMODEL->firstRegion(); it != GMODEL->lastRegion(); ++it){
- addEdgeVertices(*it, (*it)->tetrahedra, linear, edgeVertices);
- addEdgeVertices(*it, (*it)->hexahedra, linear, edgeVertices);
- addEdgeVertices(*it, (*it)->prisms, linear, edgeVertices);
- addEdgeVertices(*it, (*it)->pyramids, linear, edgeVertices);
- addFaceVertices(*it, (*it)->hexahedra, linear, faceVertices);
- addFaceVertices(*it, (*it)->prisms, linear, faceVertices);
- addFaceVertices(*it, (*it)->pyramids, linear, faceVertices);
- }
+ for(GModel::eiter it = GMODEL->firstEdge(); it != GMODEL->lastEdge(); ++it)
+ setSecondOrder(*it, edgeVertices, linear);
+ for(GModel::fiter it = GMODEL->firstFace(); it != GMODEL->lastFace(); ++it)
+ setSecondOrder(*it, edgeVertices, faceVertices, linear);
+ for(GModel::riter it = GMODEL->firstRegion(); it != GMODEL->lastRegion(); ++it)
+ setSecondOrder(*it, edgeVertices, faceVertices, linear);
- printf("%d edges, %d faces\n", edgeVertices.size(), faceVertices.size());
-
- // 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);
}