7 #define PI 3.14159265358979323846
144 int IdSubdiviedElement1;
145 int IdSubdiviedElement2;
159 IdNgbhVortex1=LocMesh->
Edges[IdEdge].
V1;
160 IdNgbhVortex2=LocMesh->
Edges[IdEdge].
V2;
163 IdSubdiviedEdge=IdEdge;
168 IdSubdiviedElement1=-1;
169 IdSubdiviedElement2=-1;
171 if (IdSubdiviedElement1==-1){
172 IdSubdiviedElement1=i;
200 IdSubdiviedElement2=i;
237 LocMesh->
Edges[IdSubdiviedEdge].
V1=IdNgbhVortex1;
238 LocMesh->
Edges[IdSubdiviedEdge].
V2=IdNewVortex;
240 LocMesh->
Edges[IdNewEdge1].
V1=IdNewVortex;
241 LocMesh->
Edges[IdNewEdge1].
V2=IdNgbhVortex2;
244 LocMesh->
Edges[IdNewEdge2].
V1=IdNewVortex;
245 LocMesh->
Edges[IdNewEdge2].
V2=IdNgbhVortex4;
248 LocMesh->
Edges[IdNewEdge3].
V1=IdNewVortex;
249 LocMesh->
Edges[IdNewEdge3].
V2=IdNgbhVortex3;
254 LocMesh->
Elements[IdSubdiviedElement1].
E1=IdSubdiviedEdge;
255 LocMesh->
Elements[IdSubdiviedElement1].
E2=IdNewEdge3;
256 LocMesh->
Elements[IdSubdiviedElement1].
E3=IdNgbhEdge1;
258 LocMesh->
Elements[IdSubdiviedElement2].
E1=IdNgbhEdge2;
259 LocMesh->
Elements[IdSubdiviedElement2].
E2=IdNewEdge2;
260 LocMesh->
Elements[IdSubdiviedElement2].
E3=IdSubdiviedEdge;
262 LocMesh->
Elements[IdNewElement1].
E1=IdNewEdge1;
263 LocMesh->
Elements[IdNewElement1].
E2=IdNgbhEdge3;
264 LocMesh->
Elements[IdNewElement1].
E3=IdNewEdge3;
266 LocMesh->
Elements[IdNewElement2].
E1=IdNgbhEdge4;
267 LocMesh->
Elements[IdNewElement2].
E2=IdNewEdge1;
268 LocMesh->
Elements[IdNewElement2].
E3=IdNewEdge2;
274 double Xmax,Xmin,Ymax,Ymin,Zmax,Zmin,Xmean,Ymean,Zmean;
303 LocNorm=sqrt((pow(LocMesh->
Vertexes[i].
X,2.0)/pow(Rx,2.0))+
304 (pow(LocMesh->
Vertexes[i].
Y,2.0)/pow(Ry,2.0))+
305 (pow(LocMesh->
Vertexes[i].
Z,2.0)/pow(Rz,2.0)));
316 int Ve1,Ed1,Ve2,Ed2,Ve3,Ed3;
318 double x1,x2,x3,y1,y2,y3,z1,z2,z3;
319 double u1,u2,u3,v1,v2,v3;
331 temp=Ed2; Ed2=Ed3; Ed3=temp;
357 printf(
"This is not good!!!\n");
365 u1=x2-x1; u2=y2-y1; u3=z2-z1;
366 v1=x3-x1; v2=y3-y1; v3=z3-z1;
375 if (x1*vp1+y1*vp2+z1*vp3<0){
377 temp=Ed3; Ed3=Ed2; Ed2=temp;
400 double SqLengthI,SqLengthJ;
403 order=(
int*)malloc(LocMesh->
NbEdges*
sizeof(
int));
405 for (i=0;i<LocMesh->
NbEdges;i++) order[i]=i;
407 for (i=0;i<LocMesh->
NbEdges-1;i++)
408 for (j=i+1;j<LocMesh->
NbEdges;j++){
416 if (SqLengthI<SqLengthJ){
424 for (i=0;i<OrigNbEdges;i++)
CutEdge(LocMesh,order[i]);
434 double SqLengthI,SqLengthJ;
440 double EdgeV1_X,EdgeV1_Y,EdgeV1_Z,EdgeV2_X,EdgeV2_Y,EdgeV2_Z;
443 order=(
int*)malloc(LocMesh->
NbEdges*
sizeof(
int));
444 for (i=0;i<LocMesh->
NbEdges;i++) order[i]=i;
454 order=(
int*)realloc(order,LocMesh->
NbEdges*
sizeof(
int));
455 for (i=0;i<LocMesh->
NbEdges;i++) order[i]=i;
458 for (i=0;i<LocMesh->
NbEdges-1;i++)
for (j=i+1;j<LocMesh->
NbEdges;j++){
466 if (SqLengthI<SqLengthJ){ k=order[i]; order[i]=order[j]; order[j]=k;}
471 for (i=0;i<OrigNbEdges;i++){
476 if (sqrt(SqLengthI)>MaxEdgeLength){
491 for (i=0;i<LocMesh->
NbEdges;i++){
509 LocMesh->
Edges[i].
SubiVertex[j].
X=EdgeV1_X*(dNbSubdiv-dj)/dNbSubdiv+EdgeV2_X*dj/dNbSubdiv;
510 LocMesh->
Edges[i].
SubiVertex[j].
Y=EdgeV1_Y*(dNbSubdiv-dj)/dNbSubdiv+EdgeV2_Y*dj/dNbSubdiv;
511 LocMesh->
Edges[i].
SubiVertex[j].
Z=EdgeV1_Z*(dNbSubdiv-dj)/dNbSubdiv+EdgeV2_Z*dj/dNbSubdiv;
527 void ExtractSubmesh(
Mesh * LocMesh,
double Xmin,
double Xmax,
double Ymin,
double Ymax,
double Zmin,
double Zmax){
528 int i,j,NbCurvedEdge;
529 int * VertexesToKeep;
531 int * ElementsToKeep;
540 VertexesToKeep=(
int*)malloc(LocMesh->
NbVertexes*
sizeof(
int));
541 for (i=0;i<LocMesh->
NbVertexes;i++) VertexesToKeep[i]=0;
543 EdgesToKeep=(
int*)malloc(LocMesh->
NbEdges*
sizeof(
int));
544 for (i=0;i<LocMesh->
NbEdges;i++) EdgesToKeep[i]=0;
546 ElementsToKeep=(
int*)malloc(LocMesh->
NbElements*
sizeof(
int));
547 for (i=0;i<LocMesh->
NbElements;i++) ElementsToKeep[i]=0;
553 (LocMesh->
Vertexes[i].
Z>Zmin-epsilon)&&(LocMesh->
Vertexes[i].
Z<Zmax+epsilon)) VertexesToKeep[i]=1;
556 for (i=0;i<LocMesh->
NbEdges;i++)
if ((VertexesToKeep[LocMesh->
Edges[i].
V1]==1)&&(VertexesToKeep[LocMesh->
Edges[i].
V2]==1))
568 for (i=0;i<LocMesh->
NbVertexes;i++)
if (VertexesToKeep[i]==1){
569 VertexesToKeep[i]=NewValue;
574 for (i=0;i<LocMesh->
NbEdges;i++)
if (EdgesToKeep[i]==1){
575 EdgesToKeep[i]=NewValue;
580 for (i=0;i<LocMesh->
NbElements;i++)
if (ElementsToKeep[i]==1){
581 ElementsToKeep[i]=NewValue;
596 for (i=0;i<TempMesh.
NbEdges;i++){
620 for (i=0;i<TempMesh.
NbVertexes;i++)
if (VertexesToKeep[i]>0){
630 for (i=0;i<TempMesh.
NbEdges;i++)
if (EdgesToKeep[i]>0){
644 for (i=0;i<TempMesh.
NbElements;i++)
if (ElementsToKeep[i]>0){
658 FILE * XmlMeshGeomFile;
659 int i,j,NbCurvedEdge;
660 int * BoundaryEdgesList;
667 XmlMeshGeomFile = fopen(FileName,
"w");
670 fprintf(XmlMeshGeomFile,
"<?xml version=\"1.0\" encoding=\"utf-8\"?>\n");
671 fprintf(XmlMeshGeomFile,
"\n");
672 fprintf(XmlMeshGeomFile,
"<NEKTAR>\n");
673 fprintf(XmlMeshGeomFile,
"<!-- Embed a 2-dimensional object in a 2-dimensional space -->\n");
674 fprintf(XmlMeshGeomFile,
"<!-- DIM <= SPACE -->\n");
675 fprintf(XmlMeshGeomFile,
"<!-- This provides a method of optimizing code for a 1-D curve embedded in 3-space. -->\n");
676 fprintf(XmlMeshGeomFile,
"<GEOMETRY DIM=\"2\" SPACE=\"3\">\n");
677 fprintf(XmlMeshGeomFile,
"\n");
680 fprintf(XmlMeshGeomFile,
" <VERTEX>\n");
681 fprintf(XmlMeshGeomFile,
" <!-- Always must have four values per entry. -->\n");
684 fprintf(XmlMeshGeomFile,
" </VERTEX>\n");
685 fprintf(XmlMeshGeomFile,
" \n");
688 fprintf(XmlMeshGeomFile,
" <EDGE>\n");
689 fprintf(XmlMeshGeomFile,
" <!--Edges are vertex pairs -->\n");
690 for (i=0;i<LocMesh->
NbEdges;i++)
691 fprintf(XmlMeshGeomFile,
" <E ID=\"%d\"> %d %d </E>\n",i,LocMesh->
Edges[i].
V1,LocMesh->
Edges[i].
V2);
692 fprintf(XmlMeshGeomFile,
" </EDGE>\n");
693 fprintf(XmlMeshGeomFile,
" \n");
696 fprintf(XmlMeshGeomFile,
" <ELEMENT>\n");
699 fprintf(XmlMeshGeomFile,
" </ELEMENT>\n");
700 fprintf(XmlMeshGeomFile,
" \n");
705 fprintf(XmlMeshGeomFile,
" <CURVED>\n");
708 fprintf(XmlMeshGeomFile,
" <E ID=\"%d\" EDGEID=\"%d\" TYPE=\"PolyEvenlySpaced\" NUMPOINTS=\"%d\">\n",NbCurvedEdge,i,LocMesh->
Edges[i].
NbSubdiv+1);
712 fprintf(XmlMeshGeomFile,
" </E>\n\n");
715 fprintf(XmlMeshGeomFile,
" </CURVED>\n");
719 BoundaryEdgesList=(
int *)malloc(LocMesh->
NbEdges*
sizeof(
int));
721 for (i=0;i<LocMesh->
NbEdges;i++) BoundaryEdgesList[i]=-1;
724 if (BoundaryEdgesList[LocMesh->
Elements[i].
E1]==1) printf(
"Problem between the elements and the edges\n");
725 if (BoundaryEdgesList[LocMesh->
Elements[i].
E2]==1) printf(
"Problem between the elements and the edges\n");
726 if (BoundaryEdgesList[LocMesh->
Elements[i].
E3]==1) printf(
"Problem between the elements and the edges\n");
727 if (BoundaryEdgesList[LocMesh->
Elements[i].
E1]==0) BoundaryEdgesList[LocMesh->
Elements[i].
E1]=1;
728 if (BoundaryEdgesList[LocMesh->
Elements[i].
E2]==0) BoundaryEdgesList[LocMesh->
Elements[i].
E2]=1;
729 if (BoundaryEdgesList[LocMesh->
Elements[i].
E3]==0) BoundaryEdgesList[LocMesh->
Elements[i].
E3]=1;
730 if (BoundaryEdgesList[LocMesh->
Elements[i].
E1]==-1) BoundaryEdgesList[LocMesh->
Elements[i].
E1]=0;
731 if (BoundaryEdgesList[LocMesh->
Elements[i].
E2]==-1) BoundaryEdgesList[LocMesh->
Elements[i].
E2]=0;
732 if (BoundaryEdgesList[LocMesh->
Elements[i].
E3]==-1) BoundaryEdgesList[LocMesh->
Elements[i].
E3]=0;
736 fprintf(XmlMeshGeomFile,
"<!-- V - vertex, E - edge, F - face, L - element -->\n");
737 fprintf(XmlMeshGeomFile,
" <COMPOSITE>\n");
738 fprintf(XmlMeshGeomFile,
" <C ID=\"0\"> T[0-%d]\n",LocMesh->
NbElements-1);
739 fprintf(XmlMeshGeomFile,
" </C>\n");
740 fprintf(XmlMeshGeomFile,
" <C ID=\"1\"> E[");
742 for (i=0;i<LocMesh->
NbEdges;i++)
if (BoundaryEdgesList[i]==0){
743 if (FirstBoundary==1) FirstBoundary=0;
744 else fprintf(XmlMeshGeomFile,
",");
745 fprintf(XmlMeshGeomFile,
"%d",i);
747 fprintf(XmlMeshGeomFile,
"]\n");
748 fprintf(XmlMeshGeomFile,
" </C>\n");
749 fprintf(XmlMeshGeomFile,
" </COMPOSITE>\n");
750 fprintf(XmlMeshGeomFile,
" \n");
751 fprintf(XmlMeshGeomFile,
" <DOMAIN> C[0] </DOMAIN>\n");
752 fprintf(XmlMeshGeomFile,
" \n");
755 fprintf(XmlMeshGeomFile,
"</GEOMETRY>\n");
758 fclose(XmlMeshGeomFile);
772 printf(
"Usage: makeEllipsoidNektar <options>\n");
773 printf(
"Where <options> are one or more of the following:\n");
774 printf(
"\t<-NameFile n> Name of the output file containg the ellipsoid (default=\"toto.xml\")\n");
775 printf(
"\t<-MaxLength n> Maximum length allowed for an edge (default=0.5)\n");
776 printf(
"\t<-NbSubdiv n> Number of subdivisions of each edge to make it curved (default=1)\n");
777 printf(
"\t<-RadiusX n> Radius of the ellipsoid in the X direction (default=1.)\n");
778 printf(
"\t<-RadiusY n> Radius of the ellipsoid in the Y direction (default=1.)\n");
779 printf(
"\t<-RadiusZ n> Radius of the ellipsoid in the Z direction (default=1.)\n");
780 printf(
"\t<-MinX n> Bounding box - min X (default=-1000.)\n");
781 printf(
"\t<-MaxX n> Bounding box - max X (default=1000.)\n");
782 printf(
"\t<-MinY n> Bounding box - min Y (default=-1000.)\n");
783 printf(
"\t<-MaxY n> Bounding box - max Y (default=1000.)\n");
784 printf(
"\t<-MinZ n> Bounding box - min Z (default=-1000.)\n");
785 printf(
"\t<-MaxZ n> Bounding box - max Z (default=1000.)\n");
789 int main(
int argc,
char **argv)
792 char *output_name =
"toto.xml";
810 if ((ok == 0) && (strcmp(argv[1],
"help") == 0)) {
813 if ((ok == 0) && (strcmp(argv[1],
"-NameFile") == 0)) {
816 output_name = argv[1];
821 if ((ok == 0) && (strcmp(argv[1],
"-RadiusX") == 0)) {
824 RadiusX = atof(argv[1]);
825 if (RadiusX<0) RadiusX=1.;
830 if ((ok == 0) && (strcmp(argv[1],
"-RadiusY") == 0)) {
833 RadiusY = atof(argv[1]);
834 if (RadiusY<0) RadiusY=1.;
839 if ((ok == 0) && (strcmp(argv[1],
"-RadiusZ") == 0)) {
842 RadiusZ = atof(argv[1]);
843 if (RadiusZ<0) RadiusZ=1.;
848 if ((ok == 0) && (strcmp(argv[1],
"-MinX") == 0)) {
851 MinX = atof(argv[1]);
856 if ((ok == 0) && (strcmp(argv[1],
"-MaxX") == 0)) {
859 MaxX = atof(argv[1]);
864 if ((ok == 0) && (strcmp(argv[1],
"-MinY") == 0)) {
867 MinY = atof(argv[1]);
872 if ((ok == 0) && (strcmp(argv[1],
"-MaxY") == 0)) {
875 MaxY = atof(argv[1]);
880 if ((ok == 0) && (strcmp(argv[1],
"-MinZ") == 0)) {
883 MinZ = atof(argv[1]);
888 if ((ok == 0) && (strcmp(argv[1],
"-MaxZ") == 0)) {
891 MaxZ = atof(argv[1]);
896 if ((ok == 0) && (strcmp(argv[1],
"-MaxLength") == 0)) {
899 Prec = atof(argv[1]);
900 if ((Prec<(RadiusX+RadiusY+RadiusZ)*0.001/3.)) Prec=(RadiusX+RadiusY+RadiusZ)*0.001/3;
905 if ((ok == 0) && (strcmp(argv[1],
"-NbSubdiv") == 0)) {
908 NbSubdiv = atoi(argv[1]);
909 if ((NbSubdiv<1)) NbSubdiv=1;
915 printf(
"Unknown option: \n");
920 printf(
"Output file: %s / Radius: %lf / Radius Y: %lf / Radius Z: %lf / Precision: %lf / Subdivisions: %d\n",output_name,RadiusX,RadiusY,RadiusZ,Prec,NbSubdiv);
for(iterator iter=lhs.begin();iter!=lhs.end();++iter)
void ProjectMeshEllipsoidAndRefine(Mesh *LocMesh, double Rx, double Ry, double Rz, double MaxEdgeLength, int NbSubdiv)
void WriteMesh(Mesh *LocMesh, char FileName[256])
int main(int argc, char **argv)
Represents a vertex in the mesh.
void MakeClockwiseOrder(Mesh *LocMesh)
void ProjectMeshEllipsoid(Mesh *LocMesh, double Rx, double Ry, double Rz)
void ExtractSubmesh(Mesh *LocMesh, double Xmin, double Xmax, double Ymin, double Ymax, double Zmin, double Zmax)
void CutEdge(Mesh *LocMesh, int IdEdge)
Mesh CreateBasicCubicMesh()
void RefineMesh(Mesh *LocMesh)