diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp index a2568a7458f95125882d95da75727e97751f2266..488015fac9b1c9077cb72d1b4a0858909f886aef 100644 --- a/Mesh/Field.cpp +++ b/Mesh/Field.cpp @@ -1,4 +1,4 @@ -// $Id: Field.cpp,v 1.36 2008-05-04 08:31:15 geuzaine Exp $ +// $Id: Field.cpp,v 1.37 2008-05-22 12:26:59 remacle Exp $ // // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle // @@ -453,7 +453,7 @@ public: double operator() (double x, double y, double z) { return (*GModel::current()->getFields()->get(field_id)) - (atan2(x, y), asin(z / sqrt(x * x + y * y + z * z)), 0); + (atan2(y, x), asin(z / sqrt(x * x + y * y + z * z)), 0); } FieldDialogBox *&dialog_box() { diff --git a/Plugin/GSHHS.cpp b/Plugin/GSHHS.cpp index d70042d48aee18c0addf6c4052fb65fd00f0132c..04fc69c04a8c2f69ed3ed94de365f901ff460177 100644 --- a/Plugin/GSHHS.cpp +++ b/Plugin/GSHHS.cpp @@ -203,7 +203,7 @@ class GeoEarthImport buff << "IFI = newf;\n"; if(write_polar_sphere){ buff << "Point ( IP + " << ip++ << " ) = {0, 0, 0 , lc};\n"; - buff << "Point ( IP + " << ip++ << " ) = {0, 0, 6.371e6 , lc};\n"; + buff << "Point ( IP + " << ip++ <<" ) = {0, 0, 6.371e6 , lc};\n"; buff << "PolarSphere ( IS + " << is++ << " ) = {IP , IP+1};\n"; } *file << buff.str(); @@ -273,7 +273,7 @@ class GeoEarthImport } }; -/* $Id: GSHHS.cpp,v 1.9 2008-05-04 08:31:23 geuzaine Exp $ +/* $Id: GSHHS.cpp,v 1.10 2008-05-22 12:26:59 remacle Exp $ * * PROGRAM: gshhs.c * AUTHOR: Paul Wessel (pwessel@hawaii.edu) @@ -378,6 +378,8 @@ void import_gshhs(FILE * fp, GeoEarthImport & geo_import) p.x = swabi4((unsigned int)p.x); p.y = swabi4((unsigned int)p.y); } + double lat=M_PI / 180*((h.greenwich && p.x > max_east) ? p.x * 1.0e-6 - 360.0 : p.x * 1.0e-6); + double lon= (p.y * 1.0e-6) * M_PI / 180; geo_import.add_point(SPoint3( M_PI / 180*((h.greenwich && p.x > max_east) ? p.x * 1.0e-6 - 360.0 : p.x * 1.0e-6), diff --git a/Plugin/TransformLatLon.cpp b/Plugin/TransformLatLon.cpp index 49bd57444a9770614ab8e99af8bca633644ab55a..ee63c382dc56a7292056f9ef659dbdf6f9cb79b2 100644 --- a/Plugin/TransformLatLon.cpp +++ b/Plugin/TransformLatLon.cpp @@ -1,4 +1,4 @@ -// $Id: TransformLatLon.cpp,v 1.5 2008-03-20 11:44:14 geuzaine Exp $ +// $Id: TransformLatLon.cpp,v 1.6 2008-05-22 12:26:59 remacle Exp $ // // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle // @@ -22,8 +22,8 @@ #include "TransformLatLon.h" StringXNumber TransformLatLonOptions_Number[] = { - {GMSH_FULLRC,(char*) "iView", NULL, -1.}, - {GMSH_FULLRC,(char*) "Longitude0", NULL, 0.} + {GMSH_FULLRC, (char *)"iView", NULL, -1.}, + {GMSH_FULLRC, (char *)"Longitude0", NULL, 0.} }; extern "C" @@ -45,7 +45,7 @@ void GMSH_TransformLatLonPlugin::getName(char *name) const } void GMSH_TransformLatLonPlugin::getInfos(char *author, char *copyright, - char *help_text) const + char *help_text) const { strcpy(author, "C. Geuzaine"); strcpy(copyright, "DGR (www.multiphysics.com)"); @@ -59,8 +59,7 @@ void GMSH_TransformLatLonPlugin::getInfos(char *author, char *copyright, "If `SwapOrientation' is set, the orientation of the\n" "elements is reversed. If `iView' < 0, the plugin\n" "is run on the current view.\n" - "\n" - "Plugin(TransformLatLon) is executed in-place.\n"); + "\n" "Plugin(TransformLatLon) is executed in-place.\n"); } int GMSH_TransformLatLonPlugin::getNbOptions() const @@ -79,112 +78,172 @@ void GMSH_TransformLatLonPlugin::catchErrorMessage(char *errorMessage) const } // TransformLatLonation - -static void transform(double v[3], double *x, double *y, double *z,double rotate) +static void vector_transform(double *vin, double *vout, double x, + double y, double z) { - double r = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); - *y=asin(v[2]/r)*180/M_PI; - *x=atan2(v[1],v[0])*180/M_PI+rotate; - while(*x>180)*x-=360; - *z=0; + double theta = y * M_PI / 180; + double phi = x * M_PI / 180; + vout[0] = -sin(phi) * vin[0] + cos(phi) * vin[1]; + vout[1] = + -sin(theta) * (cos(phi) * vin[0] + sin(phi) * vin[1]) + + cos(theta) * vin[2]; + vout[2] = + cos(theta) * (cos(phi) * vin[0] + sin(phi) * vin[1]) + + sin(theta) * vin[2]; } -static void transform_list(List_T *list, int nbList, List_T *list2, int *nbList2, - int nbVert, int nbComp,double rotate) +static void transform_list(List_T * list, int nbList, List_T * list2, + int *nbList2, int nbVert, int nbComp, + double rotate) { - if(!nbList) return; + if(!nbList) + return; int nb = List_Nbr(list) / nbList; - double *x2=new double[nbVert]; - double *y2=new double[nbVert]; - double *z2=new double[nbVert]; + double *ll_out=new double[nbVert*3]; + double *lon=ll_out; + double *lat=ll_out+nbVert; + double *r=ll_out+nbVert*2; + double *din = new double[nbComp]; + double *dout = new double[nbComp]; for(int i = 0; i < List_Nbr(list); i += nb) { double *x = (double *)List_Pointer_Fast(list, i); double *y = (double *)List_Pointer_Fast(list, i + nbVert); double *z = (double *)List_Pointer_Fast(list, i + 2 * nbVert); for(int j = 0; j < nbVert; j++) { - double v[3] = { x[j], y[j], z[j] }; - transform(v, &x2[j], &y2[j], &z2[j],rotate); + r[j] = sqrt(x[j]*x[j] + y[j]*y[j]+ z[j]*z[j]); + lat[j] = asin(z[j] / r[j]) * 180 / M_PI; + lon[j] = atan2(y[j], x[j]) * 180 / M_PI + rotate; + while(lon[j] > 180) + lon[j] -= 360; + lon[j]+=rotate; + } + // check the periodicity + bool cross = false; + for(int j = 1; j < nbVert; j++) + cross |= (fabs(lon[j] - lon[0]) > 180); + if(cross) + for(int j = 0; j < nbVert; j++) + if(lon[j] < 0) + lon[j] += 360; + for(int j = 0; j < nbVert*3; j++) + List_Add(list2, &ll_out[j]); + for(int k = 0; k < nbVert; k++) { + for(int l = 0; l < nbComp; l++) + List_Read(list, i + 3 * nbVert + nbComp * k + l, &din[l]); + if(nbComp == 1) + dout[0]=din[0]; + else if(nbComp == 3) { + vector_transform(din, dout, lon[k], lat[k], r[k]); + }else if(nbComp == 9){ + + }else{ + + } + for(int l = 0; l < nbComp; l++) + List_Add(list2, &dout[l]); } - // check the periodicity - bool cross=false; - for(int j=1;j<nbVert;j++) cross|=(fabs(x2[j]-x2[0])>180); - if(cross) for(int j=0;j<nbVert;j++) if(x2[j]<0)x2[j]+=360; - for(int j = 0; j < nbVert; j++)List_Add(list2,&x2[j]); - for(int j = 0; j < nbVert; j++)List_Add(list2,&y2[j]); - for(int j = 0; j < nbVert; j++)List_Add(list2,&z2[j]); - for(int k = 0; k < nbVert; k++){ - double d; - for(int l = 0; l < nbComp; l++){ - List_Read(list, i + 3 * nbVert + nbComp * k + l, &d); - List_Add(list2, &d); - } - } (*nbList2)++; - if(cross){ - for(int j=0;j<nbVert;j++) x2[j]-=360; - for(int j = 0; j < nbVert; j++)List_Add(list2,&x2[j]); - for(int j = 0; j < nbVert; j++)List_Add(list2,&y2[j]); - for(int j = 0; j < nbVert; j++)List_Add(list2,&z2[j]); - for(int k = 0; k < nbVert; k++){ - double d; - for(int l = 0; l < nbComp; l++){ - List_Read(list, i + 3 * nbVert + nbComp * k + l, &d); - List_Add(list2, &d); - } - } - (*nbList2)++; - } + if(cross) { + for(int j = 0; j < nbVert; j++) + lon[j] -= 360; + for(int j = 0; j < nbVert*3; j++) + List_Add(list2, &ll_out[j]); + for(int k = 0; k < nbVert; k++) { + for(int l = 0; l < nbComp; l++) + List_Read(list, i + 3 * nbVert + nbComp * k + l, &din[l]); + if(nbComp == 1){ + dout[0]=din[0]; + } + else if(nbComp == 3) { + vector_transform(din, dout, lon[k], lat[k], r[k]); + }else if(nbComp == 9){ + + }else{ + + } + for(int l = 0; l < nbComp; l++) + List_Add(list2, &dout[l]); + } + (*nbList2)++; + } } - delete []x2; - delete []y2; - delete []z2; + delete[]dout; + delete[]din; + delete[]ll_out; } -PView *GMSH_TransformLatLonPlugin::execute(PView *v) +PView *GMSH_TransformLatLonPlugin::execute(PView * v) { int iView = (int)TransformLatLonOptions_Number[0].def; - double rotate = (double)TransformLatLonOptions_Number[1].def; + double rotate = (double)TransformLatLonOptions_Number[1].def; PView *v1 = getView(iView, v); - if(!v1) return v; - PView *v2 = new PView(true); + if(!v1) + return v; + PView *v2 = new PView(true); PViewDataList *data2 = getDataList(v2); - if(!data2) return v; + if(!data2) + return v; PViewDataList *data1 = getDataList(v1); - if(!data1) return v; + if(!data1) + return v; // points - transform_list(data1->SP, data1->NbSP, data2->SP, &data2->NbSP, 1, 1, rotate); - transform_list(data1->VP, data1->NbVP, data2->VP, &data2->NbVP, 1, 3, rotate); - transform_list(data1->TP, data1->NbTP, data2->TP, &data2->NbTP, 1, 9, rotate); + transform_list(data1->SP, data1->NbSP, data2->SP, &data2->NbSP, 1, 1, + rotate); + transform_list(data1->VP, data1->NbVP, data2->VP, &data2->NbVP, 1, 3, + rotate); + transform_list(data1->TP, data1->NbTP, data2->TP, &data2->NbTP, 1, 9, + rotate); // lines - transform_list(data1->SL, data1->NbSL, data2->SL, &data2->NbSL, 2, 1, rotate); - transform_list(data1->VL, data1->NbVL, data2->VL, &data2->NbVL, 2, 3, rotate); - transform_list(data1->TL, data1->NbTL, data2->TL, &data2->NbTL, 2, 9, rotate); + transform_list(data1->SL, data1->NbSL, data2->SL, &data2->NbSL, 2, 1, + rotate); + transform_list(data1->VL, data1->NbVL, data2->VL, &data2->NbVL, 2, 3, + rotate); + transform_list(data1->TL, data1->NbTL, data2->TL, &data2->NbTL, 2, 9, + rotate); // triangles - transform_list(data1->ST, data1->NbST, data2->ST, &data2->NbST, 3, 1, rotate); - transform_list(data1->VT, data1->NbVT, data2->VT, &data2->NbVT, 3, 3, rotate); - transform_list(data1->TT, data1->NbTT, data2->TT, &data2->NbTT, 3, 9, rotate); + transform_list(data1->ST, data1->NbST, data2->ST, &data2->NbST, 3, 1, + rotate); + transform_list(data1->VT, data1->NbVT, data2->VT, &data2->NbVT, 3, 3, + rotate); + transform_list(data1->TT, data1->NbTT, data2->TT, &data2->NbTT, 3, 9, + rotate); // quadrangles - transform_list(data1->SQ, data1->NbSQ, data2->SQ, &data2->NbSQ, 4, 1, rotate); - transform_list(data1->VQ, data1->NbVQ, data2->VQ, &data2->NbVQ, 4, 3, rotate); - transform_list(data1->TQ, data1->NbTQ, data2->TQ, &data2->NbTQ, 4, 9, rotate); + transform_list(data1->SQ, data1->NbSQ, data2->SQ, &data2->NbSQ, 4, 1, + rotate); + transform_list(data1->VQ, data1->NbVQ, data2->VQ, &data2->NbVQ, 4, 3, + rotate); + transform_list(data1->TQ, data1->NbTQ, data2->TQ, &data2->NbTQ, 4, 9, + rotate); // tets - transform_list(data1->SS, data1->NbSS, data2->SS, &data2->NbSS, 4, 1, rotate); - transform_list(data1->VS, data1->NbVS, data2->VS, &data2->NbVS, 4, 3, rotate); - transform_list(data1->TS, data1->NbTS, data2->TS, &data2->NbTS, 4, 9, rotate); + transform_list(data1->SS, data1->NbSS, data2->SS, &data2->NbSS, 4, 1, + rotate); + transform_list(data1->VS, data1->NbVS, data2->VS, &data2->NbVS, 4, 3, + rotate); + transform_list(data1->TS, data1->NbTS, data2->TS, &data2->NbTS, 4, 9, + rotate); // hexas - transform_list(data1->SH, data1->NbSH, data2->SH, &data2->NbSH, 8, 1, rotate); - transform_list(data1->VH, data1->NbVH, data2->VH, &data2->NbVH, 8, 3, rotate); - transform_list(data1->TH, data1->NbTH, data2->TH, &data2->NbTH, 8, 9, rotate); + transform_list(data1->SH, data1->NbSH, data2->SH, &data2->NbSH, 8, 1, + rotate); + transform_list(data1->VH, data1->NbVH, data2->VH, &data2->NbVH, 8, 3, + rotate); + transform_list(data1->TH, data1->NbTH, data2->TH, &data2->NbTH, 8, 9, + rotate); // prisms - transform_list(data1->SI, data1->NbSI, data2->SI, &data2->NbSI, 6, 1, rotate); - transform_list(data1->VI, data1->NbVI, data2->VI, &data2->NbVI, 6, 3, rotate); - transform_list(data1->TI, data1->NbTI, data2->TI, &data2->NbTI, 6, 9, rotate); + transform_list(data1->SI, data1->NbSI, data2->SI, &data2->NbSI, 6, 1, + rotate); + transform_list(data1->VI, data1->NbVI, data2->VI, &data2->NbVI, 6, 3, + rotate); + transform_list(data1->TI, data1->NbTI, data2->TI, &data2->NbTI, 6, 9, + rotate); // pyramids - transform_list(data1->SY, data1->NbSY, data2->SY, &data2->NbSY, 5, 1, rotate); - transform_list(data1->VY, data1->NbVY, data2->VY, &data2->NbVY, 5, 3, rotate); - transform_list(data1->TY, data1->NbTY, data2->TY, &data2->NbTY, 5, 9, rotate); + transform_list(data1->SY, data1->NbSY, data2->SY, &data2->NbSY, 5, 1, + rotate); + transform_list(data1->VY, data1->NbVY, data2->VY, &data2->NbVY, 5, 3, + rotate); + transform_list(data1->TY, data1->NbTY, data2->TY, &data2->NbTY, 5, 9, + rotate); data2->setName(data1->getName() + "_LatLon"); data2->setFileName(data1->getName() + "_LatLon.pos"); diff --git a/demos/gen_earth_100km.geo b/demos/gen_earth_100km.geo index 53a6d7e4a7a63d82136a39d958c21b505676ed3a..37de0611c03a782277dcbf77b6f9fb0cdc4c63dd 100644 --- a/demos/gen_earth_100km.geo +++ b/demos/gen_earth_100km.geo @@ -1,4 +1,4 @@ -// you need gshhs data file, available on ftp://ftp.soest.hawaii.edu/pwessel/gshhs/gshhs_1.6.zip +// you need gshhs data file, available on ftp://ftp.soest.hawaii.edu/ Field[1] = Box; Field[1].VIn = 2000; Field[1].VOut = 100000;