Forked from
gmsh / gmsh
20321 commits behind the upstream repository.
-
Christophe Geuzaine authoredChristophe Geuzaine authored
3D_Extrude_Old.cpp 15.72 KiB
// $Id: 3D_Extrude_Old.cpp,v 1.18 2001-12-06 08:10:31 geuzaine Exp $
// This is the old extrusion mesh generator -> only available through
// the command line option -extrude (w/o -recombine). This mesh
// generator pre-supposes a definition of surfaces in the XY plane,
// and will extrude everything along the Z axis, taking parameters
// interactively from standard input, e.g.
//
// gmsh test -extrude -recombine < params.ext
//
// The progression ratio defines a geometric progression for the
// definition of the elements heights: a factor of 2 means that the
// next element will be twice as high than the preceding one.
//
// All geometrical entities are automatically numbered:
//
// volumes: 3 * K1 + layer * K2 + surf->num
// New XY surfaces: 2 * K1 + layer * K2 + surf->num
// perp. surfaces: 1 * K1 + layer * K2 + curve->num
// perp. lines: 4 * K1 + layer * K2 + point->Num
// New XY lines: 5 * K1 + layer * K2 + curve->Num
//
#define MAXLAYERS 100
#define K1 100.e6
#define K2 1.e6 // to store MAXLAYERS
#include "Gmsh.h"
#include "Numeric.h"
#include "Geo.h"
#include "CAD.h"
#include "Mesh.h"
#include "Context.h"
#include "Create.h"
extern Context_T CTX ;
extern Mesh *LOCAL, *THEM;
static Tree_T *Tree_Ares, *Tree_Swaps;
Surface *THES;
Volume *THEV;
int TEST_IS_ALL_OK, NbLayer;
int NbElmLayer [MAXLAYERS];
int ZonLayer [MAXLAYERS];
int LineLayer [MAXLAYERS+1];
int SurfLayer [MAXLAYERS+1];
double hLayer [MAXLAYERS];
double parLayer [MAXLAYERS];
typedef struct {
int a,b;
}nxn;
static int compnxn(const void *a, const void *b){
nxn *q,*w;
q = (nxn*)a;
w = (nxn*)b;
if(q->a > w->a)return 1;
if(q->a < w->a)return -1;
if(q->b > w->b)return 1;
if(q->b < w->b)return -1;
return 0;
}
static void InitExtrudeParams (void){
FILE *file;
int i;
printf("Number of layers: ");
scanf("%d",&NbLayer);
if(NbLayer>MAXLAYERS) Msg(FATAL, "Max number of layer (%d) exceeded", MAXLAYERS);
file = fopen("xtrude","w");
if(file) fprintf(file, "%d\n", NbLayer);
for(i=0;i<NbLayer;i++){
printf("Number of elements in layer %d: ",i+1);
scanf("%d",&NbElmLayer[i]);
if(file) fprintf(file, "%d\n", NbElmLayer[i]);
printf("Depth of layer %d: ",i+1);
scanf("%lf",&hLayer[i]);
if(file) fprintf(file, "%g\n", hLayer[i]);
printf("Progresion ratio for layer %d: ",i+1);
scanf("%lf",&parLayer[i]);
if(file) fprintf(file, "%g\n", parLayer[i]);
}
if(file){
fflush(file);
fclose(file);
}
Tree_Ares = Tree_Create(sizeof(nxn),compnxn);
Tree_Swaps = Tree_Create(sizeof(nxn),compnxn);
}
static int are_exist(Vertex *v1,Vertex *v2, Tree_T *t){
nxn n;
n.a = IMAX(v1->Num,v2->Num);
n.b = IMIN(v1->Num,v2->Num);
return Tree_Search(t,&n);
}
static void are_cree(Vertex *v1,Vertex *v2, Tree_T *t){
nxn n;
n.a = IMAX(v1->Num,v2->Num);
n.b = IMIN(v1->Num,v2->Num);
Tree_Replace(t,&n);
}
static void are_del(Vertex *v1,Vertex *v2, Tree_T *t){
nxn n;
n.a = IMAX(v1->Num,v2->Num);
n.b = IMIN(v1->Num,v2->Num);
Tree_Suppress(t,&n);
}
static void Extrude_Simplex_Phase1 (void *data , void *dum){
Simplex **pS , *s;
int i,j,k;
Vertex *v1,*v2,*v3,*v4,*v5,*v6;
pS = (Simplex**)data;
s = *pS;
k = 0;
for(i=0;i<NbLayer;i++){
for(j=0;j<NbElmLayer[i];j++){
List_Read(s->V[0]->Extruded_Points,k,&v1);
List_Read(s->V[1]->Extruded_Points,k,&v2);
List_Read(s->V[2]->Extruded_Points,k,&v3);
List_Read(s->V[0]->Extruded_Points,k+1,&v4);
List_Read(s->V[1]->Extruded_Points,k+1,&v5);
List_Read(s->V[2]->Extruded_Points,k+1,&v6); k++;
if(!are_exist(v1,v5,Tree_Ares))are_cree(v2,v4,Tree_Ares);
if(!are_exist(v5,v3,Tree_Ares))are_cree(v2,v6,Tree_Ares);
if(!are_exist(v4,v3,Tree_Ares))are_cree(v1,v6,Tree_Ares);
}
}
}
static void Extrude_Simplex_Phase3 (void *data , void *dum){
Simplex **pS , *s, *news;
Prism *newp;
Hexahedron *newh;
int i,j,k;
Vertex *v1,*v2,*v3,*v4,*v5,*v6,*v7,*v8;
pS = (Simplex**)data;
s = *pS;
if(s->V[3] && !CTX.mesh.oldxtrude_recombine){
Msg(GERROR, "Non recombined extrusion impossible with quadrangles (use -recombine)");
}
k = 0;
for(i=0;i<=NbLayer;i++){
if(SurfLayer[i]){
List_Read(s->V[0]->Extruded_Points,k,&v1);
List_Read(s->V[1]->Extruded_Points,k,&v2);
List_Read(s->V[2]->Extruded_Points,k,&v3);
news = Create_Simplex(v1,v2,v3,NULL);
if(s->V[3]){
List_Read(s->V[3]->Extruded_Points,k,&v4);
news->V[3] = v4;
}
news->iEnt = SurfLayer[i];
Tree_Add(THEV->Simp_Surf,&news);
}
for(j=0;j<NbElmLayer[i];j++){
k++;
}
}
k = 0;
for(i=0;i<NbLayer;i++){
for(j=0;j<NbElmLayer[i];j++){
if(s->V[3]){
List_Read(s->V[0]->Extruded_Points,k,&v1);
List_Read(s->V[1]->Extruded_Points,k,&v2);
List_Read(s->V[2]->Extruded_Points,k,&v3);
List_Read(s->V[3]->Extruded_Points,k,&v4);
List_Read(s->V[0]->Extruded_Points,k+1,&v5);
List_Read(s->V[1]->Extruded_Points,k+1,&v6);
List_Read(s->V[2]->Extruded_Points,k+1,&v7);
List_Read(s->V[3]->Extruded_Points,k+1,&v8);
}
else{
List_Read(s->V[0]->Extruded_Points,k,&v1);
List_Read(s->V[1]->Extruded_Points,k,&v2);
List_Read(s->V[2]->Extruded_Points,k,&v3);
List_Read(s->V[0]->Extruded_Points,k+1,&v4);
List_Read(s->V[1]->Extruded_Points,k+1,&v5);
List_Read(s->V[2]->Extruded_Points,k+1,&v6);
}
k++;
if(ZonLayer[i]){
if(CTX.mesh.oldxtrude_recombine){
if(s->V[3]){
newh = Create_Hexahedron(v1,v2,v3,v4,v5,v6,v7,v8); newh->iEnt = ZonLayer[i];
Tree_Add(THEV->Hexahedra,&newh);
}
else{
newp = Create_Prism(v1,v2,v3,v4,v5,v6);
newp->iEnt = ZonLayer[i];
Tree_Add(THEV->Prisms,&newp);
}
}
else{
if(are_exist(v4,v2,Tree_Ares) &&
are_exist(v5,v3,Tree_Ares) &&
are_exist(v4,v3,Tree_Ares) ){
news = Create_Simplex(v1,v2,v3,v4);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
news = Create_Simplex(v4,v5,v6,v3);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
news = Create_Simplex(v2,v4,v5,v3);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
}
if(are_exist(v4,v2,Tree_Ares) &&
are_exist(v2,v6,Tree_Ares) &&
are_exist(v4,v3,Tree_Ares) ){
news = Create_Simplex(v1,v2,v3,v4);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
news = Create_Simplex(v4,v5,v6,v2);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
news = Create_Simplex(v4,v2,v6,v3);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
}
if(are_exist(v4,v2,Tree_Ares) &&
are_exist(v2,v6,Tree_Ares) &&
are_exist(v6,v1,Tree_Ares) ){
news = Create_Simplex(v1,v2,v3,v6);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
news = Create_Simplex(v4,v5,v6,v2);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
news = Create_Simplex(v2,v4,v6,v1);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
}
if(are_exist(v5,v1,Tree_Ares) &&
are_exist(v5,v3,Tree_Ares) &&
are_exist(v4,v3,Tree_Ares) ){
news = Create_Simplex(v1,v2,v3,v5);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
news = Create_Simplex(v4,v5,v6,v3);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
news = Create_Simplex(v1,v4,v5,v3);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
}
if(are_exist(v5,v1,Tree_Ares) &&
are_exist(v5,v3,Tree_Ares) &&
are_exist(v6,v1,Tree_Ares) ){
news = Create_Simplex(v1,v2,v3,v5);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
news = Create_Simplex(v4,v5,v6,v1);
news->iEnt = ZonLayer[i]; Tree_Add(THEV->Simplexes,&news);
news = Create_Simplex(v1,v3,v5,v6);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
}
if(are_exist(v5,v1,Tree_Ares) &&
are_exist(v2,v6,Tree_Ares) &&
are_exist(v6,v1,Tree_Ares) ){
news = Create_Simplex(v1,v2,v3,v6);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
news = Create_Simplex(v4,v5,v6,v1);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
news = Create_Simplex(v1,v2,v5,v6);
news->iEnt = ZonLayer[i];
Tree_Add(THEV->Simplexes,&news);
}
}
}
}
}
}
static void Extrude_Simplex_Phase2 (void *data , void *dum){
Simplex **pS , *s;
int i,j,k;
Vertex *v1,*v2,*v3,*v4,*v5,*v6;
pS = (Simplex**)data;
s = *pS;
k = 0;
for(i=0;i<NbLayer;i++){
for(j=0;j<NbElmLayer[i];j++){
List_Read(s->V[0]->Extruded_Points,k,&v1);
List_Read(s->V[1]->Extruded_Points,k,&v2);
List_Read(s->V[2]->Extruded_Points,k,&v3);
List_Read(s->V[0]->Extruded_Points,k+1,&v4);
List_Read(s->V[1]->Extruded_Points,k+1,&v5);
List_Read(s->V[2]->Extruded_Points,k+1,&v6);
k++;
if(are_exist(v4,v2,Tree_Ares) &&
are_exist(v5,v3,Tree_Ares) &&
are_exist(v1,v6,Tree_Ares)){
TEST_IS_ALL_OK++;
if(!are_exist(v4,v2,Tree_Swaps)){
are_del(v4,v2,Tree_Ares);
are_cree(v1,v5,Tree_Ares);
are_cree(v1,v5,Tree_Swaps);
are_cree(v4,v2,Tree_Swaps);
}
else if(!are_exist(v5,v3,Tree_Swaps)){
are_del(v5,v3,Tree_Ares);
are_cree(v2,v6,Tree_Ares);
are_cree(v5,v3,Tree_Swaps);
are_cree(v2,v6,Tree_Swaps);
}
else if(!are_exist(v1,v6,Tree_Swaps)){
are_del(v1,v6,Tree_Ares);
are_cree(v4,v3,Tree_Ares);
are_cree(v1,v6,Tree_Swaps);
are_cree(v4,v3,Tree_Swaps);
}
}
else if(are_exist(v1,v5,Tree_Ares) &&
are_exist(v2,v6,Tree_Ares) &&
are_exist(v4,v3,Tree_Ares)){
TEST_IS_ALL_OK++; if(!are_exist(v1,v5,Tree_Swaps)){
are_del(v1,v5,Tree_Ares);
are_cree(v4,v2,Tree_Ares);
are_cree(v1,v5,Tree_Swaps);
are_cree(v4,v2,Tree_Swaps);
}
else if(!are_exist(v2,v6,Tree_Swaps)){
are_del(v2,v6,Tree_Ares);
are_cree(v5,v3,Tree_Ares);
are_cree(v5,v3,Tree_Swaps);
are_cree(v2,v6,Tree_Swaps);
}
else if(!are_exist(v4,v3,Tree_Swaps)){
are_del(v4,v3,Tree_Ares);
are_cree(v1,v6,Tree_Ares);
are_cree(v1,v6,Tree_Swaps);
are_cree(v4,v3,Tree_Swaps);
}
}
}
}
}
static void Extrude_Vertex (void *data , void *dum){
Vertex **pV , *v, *newv;
int i,j;
double h,a;
pV = (Vertex**)data;
v = *pV;
if(v->Extruded_Points)return;
v->Extruded_Points = List_Create(NbLayer,1,sizeof(Vertex*));
List_Add(v->Extruded_Points,&v);
h = 0.0;
//printf("-extruding vertex %d %p\n", v->Num, v);
for(i=0;i<NbLayer;i++){
// Geometric progression ar^i
// Sum of n (=NbElmLayer[i]) terms = hLayer[i] = a (r^n-1)/(r-1)
if(parLayer[i] == 1.)
a = hLayer[i]/(double)NbElmLayer[i];
else
a = hLayer[i] * (parLayer[i]-1.)/(pow(parLayer[i],NbElmLayer[i])-1.) ;
for(j=0;j<NbElmLayer[i];j++){
//h += hLayer[i]/(double)NbElmLayer[i];
h += a*pow(parLayer[i],j);
newv = Create_Vertex(++THEM->MaxPointNum,v->Pos.X,v->Pos.Y,v->Pos.Z + h, v->lc , v->u );
Tree_Add(THEM->Vertices,&newv);
List_Add(v->Extruded_Points,&newv);
}
}
}
static void Extrude_Surface1 (void *data , void *dum){
Surface **pS , *s;
if(!NbLayer)return;
pS = (Surface**)data;
s = THES = *pS;
Tree_Action(s->Vertices,Extrude_Vertex);
if(!CTX.mesh.oldxtrude_recombine) Tree_Action(s->Simplexes,Extrude_Simplex_Phase1);}
static void Extrude_Surface2 (void *data , void *dum){
Surface **pS , *s;
if(!NbLayer)return;
pS = (Surface**)data;
s = THES = *pS;
Tree_Action(s->Simplexes,Extrude_Simplex_Phase2);
}
static void Extrude_Surface3 (void *data , void *dum){
Surface **pS , *s;
int i;
if(!NbLayer)return;
pS = (Surface**)data;
s = THES = *pS;
/* Numerotation automatique des entites physiques */
Msg(INFO, "Extruding Surface %d", s->Num);
for(i=0;i<NbLayer;i++){
ZonLayer[i] = (int)(3 * K1) + (int) ((i+1) * K2) + s->Num ;
}
SurfLayer[0] = s->Num ;
for(i=0;i<NbLayer;i++){
SurfLayer[i+1] = (int)(2 * K1) + (int)((i+1) * K2) + s->Num ;
}
Tree_Action(s->Simplexes,Extrude_Simplex_Phase3);
}
static void Extrude_Seg(Vertex *V1, Vertex *V2){
int i,j,k;
Vertex *v1,*v2,*v3,*v4;
Simplex *s;
//printf("-curve vertex %d %p %d %p\n", V1->Num, V1, V2->Num, V2);
k = 0;
for(i=0;i<=NbLayer;i++){
if(LineLayer[i]){
List_Read(V1->Extruded_Points,k,&v1);
List_Read(V2->Extruded_Points,k,&v2);
s = Create_Simplex(v1,v2,NULL,NULL);
s->iEnt = LineLayer[i];
Tree_Add(THEV->Simp_Surf,&s);
}
for(j=0;j<NbElmLayer[i];j++){
k++;
}
}
k = 0;
for(i=0;i<NbLayer;i++){
for(j=0;j<NbElmLayer[i];j++){
List_Read(V1->Extruded_Points,k,&v1);
List_Read(V2->Extruded_Points,k,&v2);
List_Read(V1->Extruded_Points,k+1,&v3);
List_Read(V2->Extruded_Points,k+1,&v4);
if(SurfLayer[i]){
if(CTX.mesh.oldxtrude_recombine){
s = Create_Simplex(v1,v2,v4,NULL);
s->V[3] = v3;
s->iEnt = SurfLayer[i];
Tree_Add(THEV->Simp_Surf,&s);
}
else{ if(are_exist(v3,v2,Tree_Ares)){
s = Create_Simplex(v3,v2,v1,NULL);
s->iEnt = SurfLayer[i];
Tree_Add(THEV->Simp_Surf,&s);
s = Create_Simplex(v3,v4,v2,NULL);
s->iEnt = SurfLayer[i];
Tree_Add(THEV->Simp_Surf,&s);
}
else{
s = Create_Simplex(v3,v4,v1,NULL);
s->iEnt = SurfLayer[i];
Tree_Add(THEV->Simp_Surf,&s);
s = Create_Simplex(v1,v4,v2,NULL);
s->iEnt = SurfLayer[i];
Tree_Add(THEV->Simp_Surf,&s);
}
}
}
k++;
}
}
}
static void Extrude_Curve (void *data , void *dum){
Curve **pC , *c;
Vertex *v1,*v2;
int i;
if(!NbLayer)return;
pC = (Curve**)data;
c = *pC;
if (c->Num < 0) return;
/* Numerotation automatique des entites physiques */
Msg(INFO, "Extruding Curve %d", c->Num);
LineLayer[0] = c->Num ;
for(i=0;i<NbLayer;i++){
SurfLayer[i] = (int)(1 * K1) + (int)((i+1) * K2) + c->Num ;
LineLayer[i+1] = (int)(5 * K1) + (int)((i+1) * K2) + c->Num ;
}
for(i=0;i<List_Nbr(c->Vertices)-1;i++){
List_Read(c->Vertices,i,&v1);
List_Read(c->Vertices,i+1,&v2);
Extrude_Seg(v1,v2);
}
}
static void Extrude_Pnt(Vertex *V1){
int i,j,k;
Vertex *v1,*v2;
Simplex *s;
k = 0;
for(i=0;i<NbLayer;i++){
for(j=0;j<NbElmLayer[i];j++){
List_Read(V1->Extruded_Points,k,&v1);
List_Read(V1->Extruded_Points,k+1,&v2);
if(LineLayer[i]){
s = Create_Simplex(v1,v2,NULL,NULL);
s->iEnt = LineLayer[i];
Tree_Add(THEV->Simp_Surf,&s);
}
k++;
}
}
}
static void Extrude_Point (void *data , void *dum){
Vertex **pV, *v, **pV2;
int i;
if(!NbLayer)return;
pV = (Vertex**)data;
v = *pV;
/* Numerotation automatique des entites physiques */
Msg(INFO, "Extruding Vertex %d", v->Num);
for(i=0;i<NbLayer;i++){
LineLayer[i] = (int)(4 * K1) + (int)((i+1) * K2) + v->Num ;
}
if((pV2 = (Vertex**)Tree_PQuery(THEM->Vertices, pV))){
Extrude_Vertex(pV2, NULL);
Extrude_Pnt(*pV2);
}
}
void FreeEP(void *a, void *b){
Vertex *v = *((Vertex**)a);
Free_ExtrudedPoints(v->Extruded_Points);
v->Extruded_Points = NULL;
}
void Extrude_Mesh_Old(Mesh *M){
int j;
Mesh MM;
InitExtrudeParams();
LOCAL = &MM;
THEM = M;
//clean up Extruded_Points stuff (in case another extrusion was
//performed before)
Tree_Action(THEM->Vertices,FreeEP);
Create_BgMesh (WITHPOINTS, .2, LOCAL);
Tree_Left(M->Volumes,&THEV);
Tree_Action(M->Surfaces, Extrude_Surface1);
if(!CTX.mesh.oldxtrude_recombine){
j = 0;
do{
TEST_IS_ALL_OK=0;
Tree_Action ( M->Surfaces , Extrude_Surface2 );
Msg(INFO, "%d swaps", TEST_IS_ALL_OK);
if(TEST_IS_ALL_OK == j)break;
j = TEST_IS_ALL_OK;
}while(TEST_IS_ALL_OK);
}
Tree_Action ( M->Surfaces , Extrude_Surface3 );
Tree_Action ( M->Curves , Extrude_Curve );
Tree_Action ( M->Points , Extrude_Point );
}