Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Classes | Macros | Functions
makeEllipsoidNektar.c File Reference
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <float.h>
Include dependency graph for makeEllipsoidNektar.c:

Go to the source code of this file.

Classes

struct  Vertex
 Represents a vertex in the mesh. More...
struct  Edge
struct  Element
struct  Mesh

Macros

#define PI   3.14159265358979323846

Functions

Mesh CreateBasicCubicMesh ()
void CutEdge (Mesh *LocMesh, int IdEdge)
void ProjectMeshEllipsoid (Mesh *LocMesh, double Rx, double Ry, double Rz)
void MakeClockwiseOrder (Mesh *LocMesh)
void RefineMesh (Mesh *LocMesh)
void ProjectMeshEllipsoidAndRefine (Mesh *LocMesh, double Rx, double Ry, double Rz, double MaxEdgeLength, int NbSubdiv)
void ExtractSubmesh (Mesh *LocMesh, double Xmin, double Xmax, double Ymin, double Ymax, double Zmin, double Zmax)
void WriteMesh (Mesh *LocMesh, char FileName[256])
void usage ()
int main (int argc, char **argv)

Macro Definition Documentation

#define PI   3.14159265358979323846

Definition at line 7 of file makeEllipsoidNektar.c.

Function Documentation

Mesh CreateBasicCubicMesh ( )

Definition at line 63 of file makeEllipsoidNektar.c.

References Element::E1, Element::E2, Element::E3, Mesh::Edges, Mesh::Elements, Mesh::NbEdges, Mesh::NbElements, Edge::NbSubdiv, Mesh::NbVertexes, Edge::V1, Edge::V2, Mesh::Vertexes, Vertex::X, Vertex::Y, and Vertex::Z.

Referenced by main().

{
Mesh LocMesh;
int i;
//allocations...
LocMesh.NbVertexes=8;
LocMesh.Vertexes=(Vertex *)malloc(LocMesh.NbVertexes*sizeof(Vertex));
LocMesh.NbEdges=18;
LocMesh.Edges=(Edge *)malloc(LocMesh.NbEdges*sizeof(Edge));
LocMesh.NbElements=12;
LocMesh.Elements=(Element *)malloc(LocMesh.NbElements*sizeof(Element));
//fill the coordinates of each point
LocMesh.Vertexes[0].X=0; LocMesh.Vertexes[0].Y=0; LocMesh.Vertexes[0].Z=0;
LocMesh.Vertexes[1].X=1; LocMesh.Vertexes[1].Y=0; LocMesh.Vertexes[1].Z=0;
LocMesh.Vertexes[2].X=1; LocMesh.Vertexes[2].Y=1; LocMesh.Vertexes[2].Z=0;
LocMesh.Vertexes[3].X=0; LocMesh.Vertexes[3].Y=1; LocMesh.Vertexes[3].Z=0;
LocMesh.Vertexes[4].X=0; LocMesh.Vertexes[4].Y=0; LocMesh.Vertexes[4].Z=1;
LocMesh.Vertexes[5].X=1; LocMesh.Vertexes[5].Y=0; LocMesh.Vertexes[5].Z=1;
LocMesh.Vertexes[6].X=1; LocMesh.Vertexes[6].Y=1; LocMesh.Vertexes[6].Z=1;
LocMesh.Vertexes[7].X=0; LocMesh.Vertexes[7].Y=1; LocMesh.Vertexes[7].Z=1;
//fill each edge
LocMesh.Edges[0].V1=0; LocMesh.Edges[0].V2=1;
LocMesh.Edges[1].V1=1; LocMesh.Edges[1].V2=2;
LocMesh.Edges[2].V1=2; LocMesh.Edges[2].V2=3;
LocMesh.Edges[3].V1=3; LocMesh.Edges[3].V2=0;
LocMesh.Edges[4].V1=4; LocMesh.Edges[4].V2=5;
LocMesh.Edges[5].V1=5; LocMesh.Edges[5].V2=6;
LocMesh.Edges[6].V1=6; LocMesh.Edges[6].V2=7;
LocMesh.Edges[7].V1=7; LocMesh.Edges[7].V2=4;
LocMesh.Edges[8].V1=0; LocMesh.Edges[8].V2=4;
LocMesh.Edges[9].V1=1; LocMesh.Edges[9].V2=5;
LocMesh.Edges[10].V1=2; LocMesh.Edges[10].V2=6;
LocMesh.Edges[11].V1=3; LocMesh.Edges[11].V2=7;
LocMesh.Edges[12].V1=1; LocMesh.Edges[12].V2=4;
LocMesh.Edges[13].V1=1; LocMesh.Edges[13].V2=3;
LocMesh.Edges[14].V1=1; LocMesh.Edges[14].V2=6;
LocMesh.Edges[15].V1=7; LocMesh.Edges[15].V2=0;
LocMesh.Edges[16].V1=7; LocMesh.Edges[16].V2=2;
LocMesh.Edges[17].V1=7; LocMesh.Edges[17].V2=5;
for (i=0;i<LocMesh.NbEdges;i++) LocMesh.Edges[i].NbSubdiv=1;
//fill each Element
LocMesh.Elements[0].E1=0; LocMesh.Elements[0].E2=8; LocMesh.Elements[0].E3=12;
LocMesh.Elements[1].E1=9; LocMesh.Elements[1].E2=4; LocMesh.Elements[1].E3=12;
LocMesh.Elements[2].E1=0; LocMesh.Elements[2].E2=3; LocMesh.Elements[2].E3=13;
LocMesh.Elements[3].E1=1; LocMesh.Elements[3].E2=2; LocMesh.Elements[3].E3=13;
LocMesh.Elements[4].E1=1; LocMesh.Elements[4].E2=10; LocMesh.Elements[4].E3=14;
LocMesh.Elements[5].E1=9; LocMesh.Elements[5].E2=5; LocMesh.Elements[5].E3=14;
LocMesh.Elements[6].E1=3; LocMesh.Elements[6].E2=11; LocMesh.Elements[6].E3=15;
LocMesh.Elements[7].E1=8; LocMesh.Elements[7].E2=7; LocMesh.Elements[7].E3=15;
LocMesh.Elements[8].E1=2; LocMesh.Elements[8].E2=11; LocMesh.Elements[8].E3=16;
LocMesh.Elements[9].E1=10; LocMesh.Elements[9].E2=6; LocMesh.Elements[9].E3=16;
LocMesh.Elements[10].E1=4; LocMesh.Elements[10].E2=7; LocMesh.Elements[10].E3=17;
LocMesh.Elements[11].E1=5; LocMesh.Elements[11].E2=6; LocMesh.Elements[11].E3=17;
return LocMesh;
}
void CutEdge ( Mesh LocMesh,
int  IdEdge 
)

Definition at line 129 of file makeEllipsoidNektar.c.

References Element::E1, Element::E2, Element::E3, Mesh::Edges, Mesh::Elements, Mesh::NbEdges, Mesh::NbElements, Edge::NbSubdiv, Mesh::NbVertexes, Edge::V1, Edge::V2, Mesh::Vertexes, Vertex::X, Vertex::Y, and Vertex::Z.

Referenced by ProjectMeshEllipsoidAndRefine(), ProjectMeshSphereAndRefine(), and RefineMesh().

{
int i;
int IdNgbhVortex1;
int IdNgbhVortex2;
int IdNgbhVortex3;
int IdNgbhVortex4;
int IdNewVortex;
int IdSubdiviedEdge;
int IdNewEdge1;
int IdNewEdge2;
int IdNewEdge3;
int IdNgbhEdge1;
int IdNgbhEdge2;
int IdNgbhEdge3;
int IdNgbhEdge4;
int IdSubdiviedElement1;
int IdSubdiviedElement2;
int IdNewElement2;
int IdNewElement1;
//reallocations in LocMesh
LocMesh->NbVertexes++;
LocMesh->NbEdges=LocMesh->NbEdges+3;
LocMesh->NbElements=LocMesh->NbElements+2;
LocMesh->Vertexes=(Vertex *)realloc(LocMesh->Vertexes,LocMesh->NbVertexes*sizeof(Vertex));
LocMesh->Edges=(Edge *)realloc(LocMesh->Edges,LocMesh->NbEdges*sizeof(Edge));
LocMesh->Elements=(Element *)realloc(LocMesh->Elements,LocMesh->NbElements*sizeof(Element));
//identifiers to take into account
IdNgbhVortex1=LocMesh->Edges[IdEdge].V1;
IdNgbhVortex2=LocMesh->Edges[IdEdge].V2;
IdNewVortex=LocMesh->NbVertexes-1;
IdSubdiviedEdge=IdEdge;
IdNewEdge1=LocMesh->NbEdges-3;
IdNewEdge2=LocMesh->NbEdges-2;
IdNewEdge3=LocMesh->NbEdges-1;
IdSubdiviedElement1=-1;
IdSubdiviedElement2=-1;
for (i=0;i<LocMesh->NbElements-2;i++) if ((LocMesh->Elements[i].E1==IdEdge)||(LocMesh->Elements[i].E2==IdEdge)||(LocMesh->Elements[i].E3==IdEdge)){
if (IdSubdiviedElement1==-1){
IdSubdiviedElement1=i;
IdNewElement1=LocMesh->NbElements-2;
if ((LocMesh->Elements[i].E1!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E1].V1==IdNgbhVortex1)||(LocMesh->Edges[LocMesh->Elements[i].E1].V2==IdNgbhVortex1))){
IdNgbhEdge1=LocMesh->Elements[i].E1;
if (LocMesh->Edges[LocMesh->Elements[i].E1].V1==IdNgbhVortex1) IdNgbhVortex3=LocMesh->Edges[LocMesh->Elements[i].E1].V2;
else IdNgbhVortex3=LocMesh->Edges[LocMesh->Elements[i].E1].V1;
}
if ((LocMesh->Elements[i].E2!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E2].V1==IdNgbhVortex1)||(LocMesh->Edges[LocMesh->Elements[i].E2].V2==IdNgbhVortex1))){
IdNgbhEdge1=LocMesh->Elements[i].E2;
if (LocMesh->Edges[LocMesh->Elements[i].E2].V1==IdNgbhVortex1) IdNgbhVortex3=LocMesh->Edges[LocMesh->Elements[i].E2].V2;
else IdNgbhVortex3=LocMesh->Edges[LocMesh->Elements[i].E2].V1;
}
if ((LocMesh->Elements[i].E3!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E3].V1==IdNgbhVortex1)||(LocMesh->Edges[LocMesh->Elements[i].E3].V2==IdNgbhVortex1))){
IdNgbhEdge1=LocMesh->Elements[i].E3;
if (LocMesh->Edges[LocMesh->Elements[i].E3].V1==IdNgbhVortex1) IdNgbhVortex3=LocMesh->Edges[LocMesh->Elements[i].E3].V2;
else IdNgbhVortex3=LocMesh->Edges[LocMesh->Elements[i].E3].V1;
}
if ((LocMesh->Elements[i].E1!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E1].V1==IdNgbhVortex2)||(LocMesh->Edges[LocMesh->Elements[i].E1].V2==IdNgbhVortex2))){
IdNgbhEdge3=LocMesh->Elements[i].E1;
}
if ((LocMesh->Elements[i].E2!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E2].V1==IdNgbhVortex2)||(LocMesh->Edges[LocMesh->Elements[i].E2].V2==IdNgbhVortex2))){
IdNgbhEdge3=LocMesh->Elements[i].E2;
}
if ((LocMesh->Elements[i].E3!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E3].V1==IdNgbhVortex2)||(LocMesh->Edges[LocMesh->Elements[i].E3].V2==IdNgbhVortex2))){
IdNgbhEdge3=LocMesh->Elements[i].E3;
}
}
else{
IdSubdiviedElement2=i;
IdNewElement2=LocMesh->NbElements-1;
if ((LocMesh->Elements[i].E1!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E1].V1==IdNgbhVortex1)||(LocMesh->Edges[LocMesh->Elements[i].E1].V2==IdNgbhVortex1))){
IdNgbhEdge2=LocMesh->Elements[i].E1;
if (LocMesh->Edges[LocMesh->Elements[i].E1].V1==IdNgbhVortex1) IdNgbhVortex4=LocMesh->Edges[LocMesh->Elements[i].E1].V2;
else IdNgbhVortex4=LocMesh->Edges[LocMesh->Elements[i].E1].V1;
}
if ((LocMesh->Elements[i].E2!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E2].V1==IdNgbhVortex1)||(LocMesh->Edges[LocMesh->Elements[i].E2].V2==IdNgbhVortex1))){
IdNgbhEdge2=LocMesh->Elements[i].E2;
if (LocMesh->Edges[LocMesh->Elements[i].E2].V1==IdNgbhVortex1) IdNgbhVortex4=LocMesh->Edges[LocMesh->Elements[i].E2].V2;
else IdNgbhVortex4=LocMesh->Edges[LocMesh->Elements[i].E2].V1;
}
if ((LocMesh->Elements[i].E3!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E3].V1==IdNgbhVortex1)||(LocMesh->Edges[LocMesh->Elements[i].E3].V2==IdNgbhVortex1))){
IdNgbhEdge2=LocMesh->Elements[i].E3;
if (LocMesh->Edges[LocMesh->Elements[i].E3].V1==IdNgbhVortex1) IdNgbhVortex4=LocMesh->Edges[LocMesh->Elements[i].E3].V2;
else IdNgbhVortex4=LocMesh->Edges[LocMesh->Elements[i].E3].V1;
}
if ((LocMesh->Elements[i].E1!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E1].V1==IdNgbhVortex2)||(LocMesh->Edges[LocMesh->Elements[i].E1].V2==IdNgbhVortex2))){
IdNgbhEdge4=LocMesh->Elements[i].E1;
}
if ((LocMesh->Elements[i].E2!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E2].V1==IdNgbhVortex2)||(LocMesh->Edges[LocMesh->Elements[i].E2].V2==IdNgbhVortex2))){
IdNgbhEdge4=LocMesh->Elements[i].E2;
}
if ((LocMesh->Elements[i].E3!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E3].V1==IdNgbhVortex2)||(LocMesh->Edges[LocMesh->Elements[i].E3].V2==IdNgbhVortex2))){
IdNgbhEdge4=LocMesh->Elements[i].E3;
}
}
}
IdNewElement1=LocMesh->NbElements-2;
IdNewElement2=LocMesh->NbElements-1;
//Coordinates of the new vertex
LocMesh->Vertexes[IdNewVortex].X=(LocMesh->Vertexes[IdNgbhVortex1].X+LocMesh->Vertexes[IdNgbhVortex2].X)/2.;
LocMesh->Vertexes[IdNewVortex].Y=(LocMesh->Vertexes[IdNgbhVortex1].Y+LocMesh->Vertexes[IdNgbhVortex2].Y)/2.;
LocMesh->Vertexes[IdNewVortex].Z=(LocMesh->Vertexes[IdNgbhVortex1].Z+LocMesh->Vertexes[IdNgbhVortex2].Z)/2.;
//update the vertex identifiers of the edge ends
LocMesh->Edges[IdSubdiviedEdge].V1=IdNgbhVortex1;
LocMesh->Edges[IdSubdiviedEdge].V2=IdNewVortex;
LocMesh->Edges[IdNewEdge1].V1=IdNewVortex;
LocMesh->Edges[IdNewEdge1].V2=IdNgbhVortex2;
LocMesh->Edges[IdNewEdge1].NbSubdiv=1;
LocMesh->Edges[IdNewEdge2].V1=IdNewVortex;
LocMesh->Edges[IdNewEdge2].V2=IdNgbhVortex4;
LocMesh->Edges[IdNewEdge2].NbSubdiv=1;
LocMesh->Edges[IdNewEdge3].V1=IdNewVortex;
LocMesh->Edges[IdNewEdge3].V2=IdNgbhVortex3;
LocMesh->Edges[IdNewEdge3].NbSubdiv=1;
//update the edges identifiers of the elements
LocMesh->Elements[IdSubdiviedElement1].E1=IdSubdiviedEdge;
LocMesh->Elements[IdSubdiviedElement1].E2=IdNewEdge3;
LocMesh->Elements[IdSubdiviedElement1].E3=IdNgbhEdge1;
LocMesh->Elements[IdSubdiviedElement2].E1=IdNgbhEdge2;
LocMesh->Elements[IdSubdiviedElement2].E2=IdNewEdge2;
LocMesh->Elements[IdSubdiviedElement2].E3=IdSubdiviedEdge;
LocMesh->Elements[IdNewElement1].E1=IdNewEdge1;
LocMesh->Elements[IdNewElement1].E2=IdNgbhEdge3;
LocMesh->Elements[IdNewElement1].E3=IdNewEdge3;
LocMesh->Elements[IdNewElement2].E1=IdNgbhEdge4;
LocMesh->Elements[IdNewElement2].E2=IdNewEdge1;
LocMesh->Elements[IdNewElement2].E3=IdNewEdge2;
}
void ExtractSubmesh ( Mesh LocMesh,
double  Xmin,
double  Xmax,
double  Ymin,
double  Ymax,
double  Zmin,
double  Zmax 
)

Definition at line 527 of file makeEllipsoidNektar.c.

References Element::E1, Element::E2, Element::E3, Mesh::Edges, Mesh::Elements, Mesh::NbEdges, Mesh::NbElements, Edge::NbSubdiv, Mesh::NbVertexes, Edge::SubiVertex, Edge::V1, Edge::V2, Mesh::Vertexes, Vertex::X, Vertex::Y, and Vertex::Z.

Referenced by main().

{
int i,j,NbCurvedEdge;
int * VertexesToKeep;
int * EdgesToKeep;
int * ElementsToKeep;
int NewValue;
int NewI;
double epsilon;
Mesh TempMesh;
//1) allocations and initializations
epsilon=0.001;
VertexesToKeep=(int*)malloc(LocMesh->NbVertexes*sizeof(int));
for (i=0;i<LocMesh->NbVertexes;i++) VertexesToKeep[i]=0;
EdgesToKeep=(int*)malloc(LocMesh->NbEdges*sizeof(int));
for (i=0;i<LocMesh->NbEdges;i++) EdgesToKeep[i]=0;
ElementsToKeep=(int*)malloc(LocMesh->NbElements*sizeof(int));
for (i=0;i<LocMesh->NbElements;i++) ElementsToKeep[i]=0;
//2) find the vertexes to keep
for (i=0;i<LocMesh->NbVertexes;i++)
if ((LocMesh->Vertexes[i].X>Xmin-epsilon)&&(LocMesh->Vertexes[i].X<Xmax+epsilon)&&
(LocMesh->Vertexes[i].Y>Ymin-epsilon)&&(LocMesh->Vertexes[i].Y<Ymax+epsilon)&&
(LocMesh->Vertexes[i].Z>Zmin-epsilon)&&(LocMesh->Vertexes[i].Z<Zmax+epsilon)) VertexesToKeep[i]=1;
//3) find the edges to keep
for (i=0;i<LocMesh->NbEdges;i++) if ((VertexesToKeep[LocMesh->Edges[i].V1]==1)&&(VertexesToKeep[LocMesh->Edges[i].V2]==1))
EdgesToKeep[i]=1;
//4) find the elements to keep
for (i=0;i<LocMesh->NbElements;i++) if ((EdgesToKeep[LocMesh->Elements[i].E1]==1)&&(EdgesToKeep[LocMesh->Elements[i].E2]==1)&&(EdgesToKeep[LocMesh->Elements[i].E3]==1))
ElementsToKeep[i]=1;
//5) extraction of the submesh
//5.1) new IDs (ID+1 actually)
NewValue=1;
for (i=0;i<LocMesh->NbVertexes;i++) if (VertexesToKeep[i]==1){
VertexesToKeep[i]=NewValue;
NewValue++;
}
NewValue=1;
for (i=0;i<LocMesh->NbEdges;i++) if (EdgesToKeep[i]==1){
EdgesToKeep[i]=NewValue;
NewValue++;
}
NewValue=1;
for (i=0;i<LocMesh->NbElements;i++) if (ElementsToKeep[i]==1){
ElementsToKeep[i]=NewValue;
NewValue++;
}
//5.2) copy the mesh in a temporary mesh
TempMesh.NbVertexes=LocMesh->NbVertexes;
TempMesh.Vertexes=(Vertex *)malloc(TempMesh.NbVertexes*sizeof(Vertex));
for (i=0;i<TempMesh.NbVertexes;i++){
TempMesh.Vertexes[i].X=LocMesh->Vertexes[i].X;
TempMesh.Vertexes[i].Y=LocMesh->Vertexes[i].Y;
TempMesh.Vertexes[i].Z=LocMesh->Vertexes[i].Z;
}
TempMesh.NbEdges=LocMesh->NbEdges;
TempMesh.Edges=(Edge *)malloc(TempMesh.NbEdges*sizeof(Edge));
for (i=0;i<TempMesh.NbEdges;i++){
TempMesh.Edges[i].V1=LocMesh->Edges[i].V1;
TempMesh.Edges[i].V2=LocMesh->Edges[i].V2;
}
for (i=0;i<LocMesh->NbEdges;i++) if (LocMesh->Edges[i].NbSubdiv>1){
TempMesh.Edges[i].NbSubdiv=LocMesh->Edges[i].NbSubdiv;
TempMesh.Edges[i].SubiVertex=(Vertex *)malloc((TempMesh.Edges[i].NbSubdiv+1)*sizeof(Vertex));
for (j=0;j<TempMesh.Edges[i].NbSubdiv+1;j++){
TempMesh.Edges[i].SubiVertex[j].X=LocMesh->Edges[i].SubiVertex[j].X;
TempMesh.Edges[i].SubiVertex[j].Y=LocMesh->Edges[i].SubiVertex[j].Y;
TempMesh.Edges[i].SubiVertex[j].Z=LocMesh->Edges[i].SubiVertex[j].Z;
}
}
TempMesh.NbElements=LocMesh->NbElements;
TempMesh.Elements=(Element *)malloc(TempMesh.NbElements*sizeof(Element));
for (i=0;i<TempMesh.NbElements;i++){
TempMesh.Elements[i].E1=LocMesh->Elements[i].E1;
TempMesh.Elements[i].E2=LocMesh->Elements[i].E2;
TempMesh.Elements[i].E3=LocMesh->Elements[i].E3;
}
//5.3) change the mesh
NewI=0;
for (i=0;i<TempMesh.NbVertexes;i++) if (VertexesToKeep[i]>0){
LocMesh->Vertexes[NewI].X=TempMesh.Vertexes[i].X;
LocMesh->Vertexes[NewI].Y=TempMesh.Vertexes[i].Y;
LocMesh->Vertexes[NewI].Z=TempMesh.Vertexes[i].Z;
//printf("Vertex %d - %f %f %f\n",NewI,LocMesh->Vertexes[NewI].X,LocMesh->Vertexes[NewI].Y,LocMesh->Vertexes[NewI].Z);
NewI++;
}
LocMesh->NbVertexes=NewI;
NewI=0;
for (i=0;i<TempMesh.NbEdges;i++) if (EdgesToKeep[i]>0){
LocMesh->Edges[NewI].V1=VertexesToKeep[TempMesh.Edges[i].V1]-1;
LocMesh->Edges[NewI].V2=VertexesToKeep[TempMesh.Edges[i].V2]-1;
if (TempMesh.Edges[i].NbSubdiv>1) for (j=0;j<TempMesh.Edges[i].NbSubdiv+1;j++){
LocMesh->Edges[NewI].SubiVertex[j].X=TempMesh.Edges[i].SubiVertex[j].X;
LocMesh->Edges[NewI].SubiVertex[j].Y=TempMesh.Edges[i].SubiVertex[j].Y;
LocMesh->Edges[NewI].SubiVertex[j].Z=TempMesh.Edges[i].SubiVertex[j].Z;
} //we suppose that all edges have the same number of subdivisions
//printf("Edge %d - %d %d\n",NewI,LocMesh->Edges[NewI].V1,LocMesh->Edges[NewI].V2);
NewI++;
}
LocMesh->NbEdges=NewI;
NewI=0;
for (i=0;i<TempMesh.NbElements;i++) if (ElementsToKeep[i]>0){
LocMesh->Elements[NewI].E1=EdgesToKeep[TempMesh.Elements[i].E1]-1;
LocMesh->Elements[NewI].E2=EdgesToKeep[TempMesh.Elements[i].E2]-1;
LocMesh->Elements[NewI].E3=EdgesToKeep[TempMesh.Elements[i].E3]-1;
//printf("Element %d - %d %d %d \n",NewI,LocMesh->Elements[NewI].E1,LocMesh->Elements[NewI].E2,LocMesh->Elements[NewI].E3);
NewI++;
}
LocMesh->NbElements=NewI;
}
int main ( int  argc,
char **  argv 
)

Definition at line 789 of file makeEllipsoidNektar.c.

References CreateBasicCubicMesh(), ExtractSubmesh(), ProjectMeshEllipsoidAndRefine(), usage(), and WriteMesh().

{
int ok,i;
char *output_name = "toto.xml";
double RadiusX=1.;
double RadiusY=1.;
double RadiusZ=1.;
double MinX=-1000;
double MaxX=1000;
double MinY=-1000;
double MaxY=1000;
double MinZ=-1000;
double MaxZ=1000;
double Prec=0.5;
int NbSubdiv=1;
Mesh LocMesh;
//1) Parse parameters
while (argc > 1) {
ok = 0;
if ((ok == 0) && (strcmp(argv[1], "help") == 0)) {
usage();
}
if ((ok == 0) && (strcmp(argv[1], "-NameFile") == 0)) {
argc--;
argv++;
output_name = argv[1];
argc--;
argv++;
ok = 1;
}
if ((ok == 0) && (strcmp(argv[1], "-RadiusX") == 0)) {
argc--;
argv++;
RadiusX = atof(argv[1]);
if (RadiusX<0) RadiusX=1.;
argc--;
argv++;
ok = 1;
}
if ((ok == 0) && (strcmp(argv[1], "-RadiusY") == 0)) {
argc--;
argv++;
RadiusY = atof(argv[1]);
if (RadiusY<0) RadiusY=1.;
argc--;
argv++;
ok = 1;
}
if ((ok == 0) && (strcmp(argv[1], "-RadiusZ") == 0)) {
argc--;
argv++;
RadiusZ = atof(argv[1]);
if (RadiusZ<0) RadiusZ=1.;
argc--;
argv++;
ok = 1;
}
if ((ok == 0) && (strcmp(argv[1], "-MinX") == 0)) {
argc--;
argv++;
MinX = atof(argv[1]);
argc--;
argv++;
ok = 1;
}
if ((ok == 0) && (strcmp(argv[1], "-MaxX") == 0)) {
argc--;
argv++;
MaxX = atof(argv[1]);
argc--;
argv++;
ok = 1;
}
if ((ok == 0) && (strcmp(argv[1], "-MinY") == 0)) {
argc--;
argv++;
MinY = atof(argv[1]);
argc--;
argv++;
ok = 1;
}
if ((ok == 0) && (strcmp(argv[1], "-MaxY") == 0)) {
argc--;
argv++;
MaxY = atof(argv[1]);
argc--;
argv++;
ok = 1;
}
if ((ok == 0) && (strcmp(argv[1], "-MinZ") == 0)) {
argc--;
argv++;
MinZ = atof(argv[1]);
argc--;
argv++;
ok = 1;
}
if ((ok == 0) && (strcmp(argv[1], "-MaxZ") == 0)) {
argc--;
argv++;
MaxZ = atof(argv[1]);
argc--;
argv++;
ok = 1;
}
if ((ok == 0) && (strcmp(argv[1], "-MaxLength") == 0)) {
argc--;
argv++;
Prec = atof(argv[1]);
if ((Prec<(RadiusX+RadiusY+RadiusZ)*0.001/3.)) Prec=(RadiusX+RadiusY+RadiusZ)*0.001/3;
argc--;
argv++;
ok = 1;
}
if ((ok == 0) && (strcmp(argv[1], "-NbSubdiv") == 0)) {
argc--;
argv++;
NbSubdiv = atoi(argv[1]);
if ((NbSubdiv<1)) NbSubdiv=1;
argc--;
argv++;
ok = 1;
}
if (ok == 0) {
printf("Unknown option: \n");
usage();
}
}
printf("Output file: %s / Radius: %lf / Radius Y: %lf / Radius Z: %lf / Precision: %lf / Subdivisions: %d\n",output_name,RadiusX,RadiusY,RadiusZ,Prec,NbSubdiv);
//2) run functions
ProjectMeshEllipsoidAndRefine(&LocMesh,RadiusX,RadiusY,RadiusZ,Prec,NbSubdiv);
ExtractSubmesh(&LocMesh,MinX,MaxX,MinY,MaxY,MinZ,MaxZ);
WriteMesh(&LocMesh,output_name);
return 0;
}
void MakeClockwiseOrder ( Mesh LocMesh)

Definition at line 314 of file makeEllipsoidNektar.c.

References Element::E1, Element::E2, Element::E3, Mesh::Edges, Mesh::Elements, Mesh::NbElements, Edge::NbSubdiv, Edge::SubiVertex, Edge::V1, Edge::V2, Mesh::Vertexes, Vertex::X, Vertex::Y, and Vertex::Z.

Referenced by WriteMesh().

{
int i,j,jswap;
int Ve1,Ed1,Ve2,Ed2,Ve3,Ed3; //Order: Vertex 1 -> Edge1 -> Vertex 2 -> Edge 2 -> Vertex 3 -> Edge 3 -> Vertex 1 ...
int temp;
double x1,x2,x3,y1,y2,y3,z1,z2,z3; //coordinates
double u1,u2,u3,v1,v2,v3; //vectors
double vp1,vp2,vp3; //cross product
double tempDbl;
//test all elements of the mesh
for (i=0;i<LocMesh->NbElements;i++){
//order the vertexes and edges
Ed1=LocMesh->Elements[i].E1;
Ed2=LocMesh->Elements[i].E2;
Ed3=LocMesh->Elements[i].E3;
if ((LocMesh->Edges[Ed1].V2!=LocMesh->Edges[Ed2].V1)&&(LocMesh->Edges[Ed1].V2!=LocMesh->Edges[Ed2].V2)){ //swap edges
temp=Ed2; Ed2=Ed3; Ed3=temp;
}
if (LocMesh->Edges[Ed1].V2!=LocMesh->Edges[Ed2].V1){ //swap vertexes
temp=LocMesh->Edges[Ed2].V1; LocMesh->Edges[Ed2].V1=LocMesh->Edges[Ed2].V2; LocMesh->Edges[Ed2].V2=temp;
if (LocMesh->Edges[Ed2].NbSubdiv>1)
for (j=0;j<=LocMesh->Edges[Ed2].NbSubdiv/2;j++){
jswap=LocMesh->Edges[Ed2].NbSubdiv-j;
tempDbl=LocMesh->Edges[Ed2].SubiVertex[j].X; LocMesh->Edges[Ed2].SubiVertex[j].X=LocMesh->Edges[Ed2].SubiVertex[jswap].X; LocMesh->Edges[Ed2].SubiVertex[jswap].X=tempDbl;
tempDbl=LocMesh->Edges[Ed2].SubiVertex[j].Y; LocMesh->Edges[Ed2].SubiVertex[j].Y=LocMesh->Edges[Ed2].SubiVertex[jswap].Y; LocMesh->Edges[Ed2].SubiVertex[jswap].Y=tempDbl;
tempDbl=LocMesh->Edges[Ed2].SubiVertex[j].Z; LocMesh->Edges[Ed2].SubiVertex[j].Z=LocMesh->Edges[Ed2].SubiVertex[jswap].Z; LocMesh->Edges[Ed2].SubiVertex[jswap].Z=tempDbl;
}
}
if (LocMesh->Edges[Ed2].V2!=LocMesh->Edges[Ed3].V1){ //swap vertexes
temp=LocMesh->Edges[Ed3].V1; LocMesh->Edges[Ed3].V1=LocMesh->Edges[Ed3].V2; LocMesh->Edges[Ed3].V2=temp;
if (LocMesh->Edges[Ed3].NbSubdiv>1)
for (j=0;j<=LocMesh->Edges[Ed3].NbSubdiv/2;j++){
jswap=LocMesh->Edges[Ed3].NbSubdiv-j;
tempDbl=LocMesh->Edges[Ed3].SubiVertex[j].X; LocMesh->Edges[Ed3].SubiVertex[j].X=LocMesh->Edges[Ed3].SubiVertex[jswap].X; LocMesh->Edges[Ed3].SubiVertex[jswap].X=tempDbl;
tempDbl=LocMesh->Edges[Ed3].SubiVertex[j].Y; LocMesh->Edges[Ed3].SubiVertex[j].Y=LocMesh->Edges[Ed3].SubiVertex[jswap].Y; LocMesh->Edges[Ed3].SubiVertex[jswap].Y=tempDbl;
tempDbl=LocMesh->Edges[Ed3].SubiVertex[j].Z; LocMesh->Edges[Ed3].SubiVertex[j].Z=LocMesh->Edges[Ed3].SubiVertex[jswap].Z; LocMesh->Edges[Ed3].SubiVertex[jswap].Z=tempDbl;
}
}
if (LocMesh->Edges[Ed3].V2!=LocMesh->Edges[Ed1].V1){
printf("This is not good!!!\n");
}
//invert the order if necessary
x1=LocMesh->Vertexes[LocMesh->Edges[Ed1].V1].X; y1=LocMesh->Vertexes[LocMesh->Edges[Ed1].V1].Y; z1=LocMesh->Vertexes[LocMesh->Edges[Ed1].V1].Z;
x2=LocMesh->Vertexes[LocMesh->Edges[Ed1].V2].X; y2=LocMesh->Vertexes[LocMesh->Edges[Ed1].V2].Y; z2=LocMesh->Vertexes[LocMesh->Edges[Ed1].V2].Z;
x3=LocMesh->Vertexes[LocMesh->Edges[Ed3].V1].X; y3=LocMesh->Vertexes[LocMesh->Edges[Ed3].V1].Y; z3=LocMesh->Vertexes[LocMesh->Edges[Ed3].V1].Z;
u1=x2-x1; u2=y2-y1; u3=z2-z1;
v1=x3-x1; v2=y3-y1; v3=z3-z1;
vp1=u2*v3-u3*v2;
vp2=u3*v1-u1*v3;
vp3=u1*v2-u2*v1;
//printf("%3.2lf %3.2lf %3.2lf | %3.2lf %3.2lf %3.2lf | %3.2lf %3.2lf %3.2lf || %3.2lf %3.2lf %3.2lf || %3.2lf %3.2lf %3.2lf | %lf\n",x1,y1,z1,x2,y2,z2,x3,y3,z3,u1,u2,u3,v1,v2,v3,x1*vp1+y1*vp2+z1*vp3);
//printf("%3.2lf %3.2lf %3.2lf | %3.2lf %3.2lf %3.2lf | %3.2lf\n",x1,y1,z1,vp1,vp2,vp3,x1*vp1+y1*vp2+z1*vp3);
if (x1*vp1+y1*vp2+z1*vp3<0){//we consider that the origin is within the spheric shape so the vector (x1,y1,z1) points out of the shape
//swap edges
temp=Ed3; Ed3=Ed2; Ed2=temp;
//swap vertexes
temp=LocMesh->Edges[Ed1].V1; LocMesh->Edges[Ed1].V1=LocMesh->Edges[Ed1].V2; LocMesh->Edges[Ed1].V2=temp;
if (LocMesh->Edges[Ed1].NbSubdiv>1)
for (j=0;j<=LocMesh->Edges[Ed1].NbSubdiv/2;j++){
jswap=LocMesh->Edges[Ed1].NbSubdiv-j;
tempDbl=LocMesh->Edges[Ed1].SubiVertex[j].X; LocMesh->Edges[Ed1].SubiVertex[j].X=LocMesh->Edges[Ed1].SubiVertex[jswap].X; LocMesh->Edges[Ed1].SubiVertex[jswap].X=tempDbl;
tempDbl=LocMesh->Edges[Ed1].SubiVertex[j].Y; LocMesh->Edges[Ed1].SubiVertex[j].Y=LocMesh->Edges[Ed1].SubiVertex[jswap].Y; LocMesh->Edges[Ed1].SubiVertex[jswap].Y=tempDbl;
tempDbl=LocMesh->Edges[Ed1].SubiVertex[j].Z; LocMesh->Edges[Ed1].SubiVertex[j].Z=LocMesh->Edges[Ed1].SubiVertex[jswap].Z; LocMesh->Edges[Ed1].SubiVertex[jswap].Z=tempDbl;
}
}
//new order of the elements
LocMesh->Elements[i].E1=Ed1;
LocMesh->Elements[i].E2=Ed2;
LocMesh->Elements[i].E3=Ed3;
}
}
void ProjectMeshEllipsoid ( Mesh LocMesh,
double  Rx,
double  Ry,
double  Rz 
)

Definition at line 272 of file makeEllipsoidNektar.c.

References Mesh::NbVertexes, Mesh::Vertexes, Vertex::X, Vertex::Y, and Vertex::Z.

Referenced by ProjectMeshEllipsoidAndRefine().

{
int i;
double Xmax,Xmin,Ymax,Ymin,Zmax,Zmin,Xmean,Ymean,Zmean;
double LocNorm;
//translate the mesh centre to the origin
Xmax=LocMesh->Vertexes[0].X; Xmin=LocMesh->Vertexes[0].X;
Ymax=LocMesh->Vertexes[0].Y; Ymin=LocMesh->Vertexes[0].Y;
Zmax=LocMesh->Vertexes[0].Z; Zmin=LocMesh->Vertexes[0].Z;
for (i=1;i<LocMesh->NbVertexes;i++){
if (LocMesh->Vertexes[i].X>Xmax) Xmax=LocMesh->Vertexes[i].X;
if (LocMesh->Vertexes[i].X<Xmin) Xmin=LocMesh->Vertexes[i].X;
if (LocMesh->Vertexes[i].Y>Ymax) Ymax=LocMesh->Vertexes[i].Y;
if (LocMesh->Vertexes[i].Y<Ymin) Ymin=LocMesh->Vertexes[i].Y;
if (LocMesh->Vertexes[i].Z>Zmax) Zmax=LocMesh->Vertexes[i].Z;
if (LocMesh->Vertexes[i].Z<Zmin) Zmin=LocMesh->Vertexes[i].Z;
}
Xmean=(Xmax+Xmin)/2;
Ymean=(Ymax+Ymin)/2;
Zmean=(Zmax+Zmin)/2;
for (i=0;i<LocMesh->NbVertexes;i++){
LocMesh->Vertexes[i].X-=Xmean;
LocMesh->Vertexes[i].Y-=Ymean;
LocMesh->Vertexes[i].Z-=Zmean;
}
//project each vertex on a ellipsoid of radius 'Rx,Ry,Rz'
for (i=0;i<LocMesh->NbVertexes;i++){
LocNorm=sqrt((pow(LocMesh->Vertexes[i].X,2.0)/pow(Rx,2.0))+
(pow(LocMesh->Vertexes[i].Y,2.0)/pow(Ry,2.0))+
(pow(LocMesh->Vertexes[i].Z,2.0)/pow(Rz,2.0)));
LocMesh->Vertexes[i].X=LocMesh->Vertexes[i].X/LocNorm;
LocMesh->Vertexes[i].Y=LocMesh->Vertexes[i].Y/LocNorm;
LocMesh->Vertexes[i].Z=LocMesh->Vertexes[i].Z/LocNorm;
}
}
void ProjectMeshEllipsoidAndRefine ( Mesh LocMesh,
double  Rx,
double  Ry,
double  Rz,
double  MaxEdgeLength,
int  NbSubdiv 
)

Definition at line 431 of file makeEllipsoidNektar.c.

References CutEdge(), Mesh::Edges, Nektar::for(), Mesh::NbEdges, Edge::NbSubdiv, ProjectMeshEllipsoid(), Edge::SubiVertex, Edge::V1, Edge::V2, Mesh::Vertexes, Vertex::X, Vertex::Y, and Vertex::Z.

Referenced by main().

{
int i,j,k;
int * order;
double SqLengthI,SqLengthJ;
int OrigNbEdges;
int OK;
double LocNorm;
double dj,dNbSubdiv;
int jm1;
double EdgeV1_X,EdgeV1_Y,EdgeV1_Z,EdgeV2_X,EdgeV2_Y,EdgeV2_Z;
//initialize
order=(int*)malloc(LocMesh->NbEdges*sizeof(int));
for (i=0;i<LocMesh->NbEdges;i++) order[i]=i;
OK=1;
//mesh projection on an ellipsoid
ProjectMeshEllipsoid(LocMesh,Rx,Ry,Rz);
//Big loop to subdivide the mesh until each edge has a length smaller than 'MaxEdgeLength'
while (OK==1){
OK=0;
order=(int*)realloc(order,LocMesh->NbEdges*sizeof(int));
for (i=0;i<LocMesh->NbEdges;i++) order[i]=i;
//order the edges as a function of their length
for (i=0;i<LocMesh->NbEdges-1;i++) for (j=i+1;j<LocMesh->NbEdges;j++){
SqLengthI=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].X,2.0);
SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Y,2.0);
SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Z,2.0);
SqLengthJ=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].X,2.0);
SqLengthJ+=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].Y,2.0);
SqLengthJ+=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].Z,2.0);
if (SqLengthI<SqLengthJ){ k=order[i]; order[i]=order[j]; order[j]=k;}
}
//subdivide the edges if their length is larger than 'MaxEdgeLength'
OrigNbEdges=LocMesh->NbEdges;
for (i=0;i<OrigNbEdges;i++){
SqLengthI=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].X,2.0);
SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Y,2.0);
SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Z,2.0);
if (sqrt(SqLengthI)>MaxEdgeLength){
//subdvision
OK=1;
CutEdge(LocMesh,order[i]); //rem: 'LocMesh' is changed here. In particular LocMesh->NbEdges is larger.
}
}
//projection of the new edges on an ellipsoid
ProjectMeshEllipsoid(LocMesh,Rx,Ry,Rz);
}
//subdivide the edges end project them on the ellipsoid if required
if (NbSubdiv>1){
for (i=0;i<LocMesh->NbEdges;i++){
//init
LocMesh->Edges[i].NbSubdiv=NbSubdiv;
LocMesh->Edges[i].SubiVertex=(Vertex *)malloc((LocMesh->Edges[i].NbSubdiv+1)*sizeof(Vertex));
//coordinates of the subdivised edge
EdgeV1_X=LocMesh->Vertexes[LocMesh->Edges[i].V1].X;
EdgeV1_Y=LocMesh->Vertexes[LocMesh->Edges[i].V1].Y;
EdgeV1_Z=LocMesh->Vertexes[LocMesh->Edges[i].V1].Z;
EdgeV2_X=LocMesh->Vertexes[LocMesh->Edges[i].V2].X;
EdgeV2_Y=LocMesh->Vertexes[LocMesh->Edges[i].V2].Y;
EdgeV2_Z=LocMesh->Vertexes[LocMesh->Edges[i].V2].Z;
dNbSubdiv=(double)LocMesh->Edges[i].NbSubdiv;
for (j=0;j<LocMesh->Edges[i].NbSubdiv+1;j++){
dj=(double)j;
//interpolation
LocMesh->Edges[i].SubiVertex[j].X=EdgeV1_X*(dNbSubdiv-dj)/dNbSubdiv+EdgeV2_X*dj/dNbSubdiv;
LocMesh->Edges[i].SubiVertex[j].Y=EdgeV1_Y*(dNbSubdiv-dj)/dNbSubdiv+EdgeV2_Y*dj/dNbSubdiv;
LocMesh->Edges[i].SubiVertex[j].Z=EdgeV1_Z*(dNbSubdiv-dj)/dNbSubdiv+EdgeV2_Z*dj/dNbSubdiv;
//projection on an ellipsoid
LocNorm=sqrt((pow(LocMesh->Edges[i].SubiVertex[j].X,2.0)/pow(Rx,2.0))+
(pow(LocMesh->Edges[i].SubiVertex[j].Y,2.0)/pow(Ry,2.0))+
(pow(LocMesh->Edges[i].SubiVertex[j].Z,2.0)/pow(Rz,2.0)));
LocMesh->Edges[i].SubiVertex[j].X=LocMesh->Edges[i].SubiVertex[j].X/LocNorm;
LocMesh->Edges[i].SubiVertex[j].Y=LocMesh->Edges[i].SubiVertex[j].Y/LocNorm;
LocMesh->Edges[i].SubiVertex[j].Z=LocMesh->Edges[i].SubiVertex[j].Z/LocNorm;
}
}
}
}
void RefineMesh ( Mesh LocMesh)

Definition at line 397 of file makeEllipsoidNektar.c.

References CutEdge(), Mesh::Edges, Mesh::NbEdges, Edge::V1, Edge::V2, Mesh::Vertexes, Vertex::X, Vertex::Y, and Vertex::Z.

{
int i,j,k;
int * order;
double SqLengthI,SqLengthJ;
int OrigNbEdges;
order=(int*)malloc(LocMesh->NbEdges*sizeof(int));
for (i=0;i<LocMesh->NbEdges;i++) order[i]=i;
for (i=0;i<LocMesh->NbEdges-1;i++)
for (j=i+1;j<LocMesh->NbEdges;j++){
SqLengthI=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].X,2.0);
SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Y,2.0);
SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Z,2.0);
SqLengthJ=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].X,2.0);
SqLengthJ+=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].Y,2.0);
SqLengthJ+=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].Z,2.0);
if (SqLengthI<SqLengthJ){
k=order[i];
order[i]=order[j];
order[j]=k;
}
}
OrigNbEdges=LocMesh->NbEdges;
for (i=0;i<OrigNbEdges;i++) CutEdge(LocMesh,order[i]);
}
void usage ( )

Definition at line 770 of file makeEllipsoidNektar.c.

Referenced by main().

{
printf("Usage: makeEllipsoidNektar <options>\n");
printf("Where <options> are one or more of the following:\n");
printf("\t<-NameFile n> Name of the output file containg the ellipsoid (default=\"toto.xml\")\n");
printf("\t<-MaxLength n> Maximum length allowed for an edge (default=0.5)\n");
printf("\t<-NbSubdiv n> Number of subdivisions of each edge to make it curved (default=1)\n");
printf("\t<-RadiusX n> Radius of the ellipsoid in the X direction (default=1.)\n");
printf("\t<-RadiusY n> Radius of the ellipsoid in the Y direction (default=1.)\n");
printf("\t<-RadiusZ n> Radius of the ellipsoid in the Z direction (default=1.)\n");
printf("\t<-MinX n> Bounding box - min X (default=-1000.)\n");
printf("\t<-MaxX n> Bounding box - max X (default=1000.)\n");
printf("\t<-MinY n> Bounding box - min Y (default=-1000.)\n");
printf("\t<-MaxY n> Bounding box - max Y (default=1000.)\n");
printf("\t<-MinZ n> Bounding box - min Z (default=-1000.)\n");
printf("\t<-MaxZ n> Bounding box - max Z (default=1000.)\n");
exit(1);
}
void WriteMesh ( Mesh LocMesh,
char  FileName[256] 
)

Definition at line 657 of file makeEllipsoidNektar.c.

References Element::E1, Element::E2, Element::E3, Mesh::Edges, Mesh::Elements, MakeClockwiseOrder(), Mesh::NbEdges, Mesh::NbElements, Edge::NbSubdiv, Mesh::NbVertexes, Edge::SubiVertex, Edge::V1, Edge::V2, Mesh::Vertexes, Vertex::X, Vertex::Y, and Vertex::Z.

Referenced by main().

{
FILE * XmlMeshGeomFile;
int i,j,NbCurvedEdge;
int * BoundaryEdgesList;
int FirstBoundary;
//to have a proper order of the edges within the elements
//open file
XmlMeshGeomFile = fopen(FileName,"w");
//write header
fprintf(XmlMeshGeomFile,"<?xml version=\"1.0\" encoding=\"utf-8\"?>\n");
fprintf(XmlMeshGeomFile,"\n");
fprintf(XmlMeshGeomFile,"<NEKTAR>\n");
fprintf(XmlMeshGeomFile,"<!-- Embed a 2-dimensional object in a 2-dimensional space -->\n");
fprintf(XmlMeshGeomFile,"<!-- DIM <= SPACE -->\n");
fprintf(XmlMeshGeomFile,"<!-- This provides a method of optimizing code for a 1-D curve embedded in 3-space. -->\n");
fprintf(XmlMeshGeomFile,"<GEOMETRY DIM=\"2\" SPACE=\"3\">\n");
fprintf(XmlMeshGeomFile,"\n");
//write vertex
fprintf(XmlMeshGeomFile," <VERTEX>\n");
fprintf(XmlMeshGeomFile," <!-- Always must have four values per entry. -->\n");
for (i=0;i<LocMesh->NbVertexes;i++)
fprintf(XmlMeshGeomFile," <V ID=\"%d\"> %f %f %f </V>\n",i,LocMesh->Vertexes[i].X,LocMesh->Vertexes[i].Y,LocMesh->Vertexes[i].Z);
fprintf(XmlMeshGeomFile," </VERTEX>\n");
fprintf(XmlMeshGeomFile," \n");
//write edges
fprintf(XmlMeshGeomFile," <EDGE>\n");
fprintf(XmlMeshGeomFile," <!--Edges are vertex pairs -->\n");
for (i=0;i<LocMesh->NbEdges;i++)
fprintf(XmlMeshGeomFile," <E ID=\"%d\"> %d %d </E>\n",i,LocMesh->Edges[i].V1,LocMesh->Edges[i].V2);
fprintf(XmlMeshGeomFile," </EDGE>\n");
fprintf(XmlMeshGeomFile," \n");
//write Element
fprintf(XmlMeshGeomFile," <ELEMENT>\n");
for (i=0;i<LocMesh->NbElements;i++)
fprintf(XmlMeshGeomFile," <T ID=\"%d\"> %d %d %d </T>\n",i,LocMesh->Elements[i].E1,LocMesh->Elements[i].E2,LocMesh->Elements[i].E3);
fprintf(XmlMeshGeomFile," </ELEMENT>\n");
fprintf(XmlMeshGeomFile," \n");
//write curved elements
fprintf(XmlMeshGeomFile," <CURVED>\n");
NbCurvedEdge=0;
for (i=0;i<LocMesh->NbEdges;i++) if (LocMesh->Edges[i].NbSubdiv>1){
fprintf(XmlMeshGeomFile," <E ID=\"%d\" EDGEID=\"%d\" TYPE=\"PolyEvenlySpaced\" NUMPOINTS=\"%d\">\n",NbCurvedEdge,i,LocMesh->Edges[i].NbSubdiv+1);
for (j=0;j<LocMesh->Edges[i].NbSubdiv+1;j++){
fprintf(XmlMeshGeomFile," %lf %lf %lf\n",LocMesh->Edges[i].SubiVertex[j].X,LocMesh->Edges[i].SubiVertex[j].Y,LocMesh->Edges[i].SubiVertex[j].Z);
}
fprintf(XmlMeshGeomFile," </E>\n\n");
NbCurvedEdge++;
}
fprintf(XmlMeshGeomFile," </CURVED>\n");
//find the boundaries (points at 0 in BoundaryEdgesList)
BoundaryEdgesList=(int *)malloc(LocMesh->NbEdges*sizeof(int));
for (i=0;i<LocMesh->NbEdges;i++) BoundaryEdgesList[i]=-1;
for (i=0;i<LocMesh->NbElements;i++){
if (BoundaryEdgesList[LocMesh->Elements[i].E1]==1) printf("Problem between the elements and the edges\n");
if (BoundaryEdgesList[LocMesh->Elements[i].E2]==1) printf("Problem between the elements and the edges\n");
if (BoundaryEdgesList[LocMesh->Elements[i].E3]==1) printf("Problem between the elements and the edges\n");
if (BoundaryEdgesList[LocMesh->Elements[i].E1]==0) BoundaryEdgesList[LocMesh->Elements[i].E1]=1;
if (BoundaryEdgesList[LocMesh->Elements[i].E2]==0) BoundaryEdgesList[LocMesh->Elements[i].E2]=1;
if (BoundaryEdgesList[LocMesh->Elements[i].E3]==0) BoundaryEdgesList[LocMesh->Elements[i].E3]=1;
if (BoundaryEdgesList[LocMesh->Elements[i].E1]==-1) BoundaryEdgesList[LocMesh->Elements[i].E1]=0;
if (BoundaryEdgesList[LocMesh->Elements[i].E2]==-1) BoundaryEdgesList[LocMesh->Elements[i].E2]=0;
if (BoundaryEdgesList[LocMesh->Elements[i].E3]==-1) BoundaryEdgesList[LocMesh->Elements[i].E3]=0;
}
//write the end (could be more evolved)
fprintf(XmlMeshGeomFile,"<!-- V - vertex, E - edge, F - face, L - element -->\n");
fprintf(XmlMeshGeomFile," <COMPOSITE>\n");
fprintf(XmlMeshGeomFile," <C ID=\"0\"> T[0-%d]\n",LocMesh->NbElements-1);
fprintf(XmlMeshGeomFile," </C>\n");
fprintf(XmlMeshGeomFile," <C ID=\"1\"> E[");
FirstBoundary=1;
for (i=0;i<LocMesh->NbEdges;i++) if (BoundaryEdgesList[i]==0){
if (FirstBoundary==1) FirstBoundary=0;
else fprintf(XmlMeshGeomFile,",");
fprintf(XmlMeshGeomFile,"%d",i);
}
fprintf(XmlMeshGeomFile,"]\n");
fprintf(XmlMeshGeomFile," </C>\n");
fprintf(XmlMeshGeomFile," </COMPOSITE>\n");
fprintf(XmlMeshGeomFile," \n");
fprintf(XmlMeshGeomFile," <DOMAIN> C[0] </DOMAIN>\n");
fprintf(XmlMeshGeomFile," \n");
//EOF
fprintf(XmlMeshGeomFile,"</GEOMETRY>\n");
//close file
fclose(XmlMeshGeomFile);
}