Skip to content
Snippets Groups Projects
Commit 81f72610 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files 

second order
parent 1a62b35f
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Fltk/Callbacks.cpp
+ 7
− 7
View file @ 81f72610
// $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());
......
Geo/MRep.h
+ 8
− 1
View file @ 81f72610
......@@ -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;
......
Mesh/Makefile
+ 5
− 3
View file @ 81f72610
# $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 \
......
Mesh/SecondOrder.cpp
+ 209
− 127
View file @ 81f72610
// $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){
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 || ge->geomType() == GEntity::DiscreteCurve){
SPoint3 pc = edge.barycenter();
v = new MVertex(pc.x(), pc.y(), pc.z(), e);
v = new MVertex(pc.x(), pc.y(), pc.z(), ge);
}
else{
double u0 = 1e6, u1 = 1e6;
Range<double> bounds = e->parBounds(0);
if(e->getBeginVertex()->mesh_vertices.size() &&
e->getBeginVertex()->mesh_vertices[0] == v0)
Range<double> bounds = ge->parBounds(0);
if(ge->getBeginVertex() && ge->getBeginVertex()->mesh_vertices[0] == v0)
u0 = bounds.low();
else if(e->getEndVertex()->mesh_vertices.size() &&
e->getEndVertex()->mesh_vertices[0] == v0)
else if(ge->getEndVertex() && ge->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)
if(ge->getBeginVertex() && ge->getBeginVertex()->mesh_vertices[0] == v1)
u1 = bounds.low();
else if(e->getEndVertex()->mesh_vertices.size() &&
e->getEndVertex()->mesh_vertices[0] == v1)
else if(ge->getEndVertex() && ge->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);
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(), e);
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){
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(), f);
v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
}
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);
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);
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(), r);
r->mesh_vertices.push_back(v);
return v;
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{
//xyz2uv();
// interpolate
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);
}
}
//f->mesh_vertices.push_back(v);
return 0;
}
MVertex *addFaceVertex(GRegion *r, std::vector<MVertex*> &p, bool linear)
void setSecondOrder(GEdge *ge,
std::map<std::pair<MVertex*,MVertex*>, MVertex* > &edgeVertices,
bool linear)
{
//MVertex *v;
// just interpolate linear between the 4
//r->mesh_vertices.push_back(v);
return 0;
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;
}
ge->lines = l2;
ge->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(GFace *gf,
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<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();
}
void setSecondOrder(GRegion *gr,
std::map<std::pair<MVertex*,MVertex*>, MVertex* > &edgeVertices,
std::map<std::vector<MVertex*>, MVertex* > &faceVertices,
bool 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);
}
printf("%d edges, %d faces\n", edgeVertices.size(), faceVertices.size());
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);
// 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);
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment