Skip to content
Snippets Groups Projects
Commit f426ce95 authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

*** empty log message ***

parent a5be432d
No related branches found
No related tags found
No related merge requests found
// $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 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
// //
...@@ -453,7 +453,7 @@ public: ...@@ -453,7 +453,7 @@ public:
double operator() (double x, double y, double z) double operator() (double x, double y, double z)
{ {
return (*GModel::current()->getFields()->get(field_id)) 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() FieldDialogBox *&dialog_box()
{ {
......
...@@ -273,7 +273,7 @@ class GeoEarthImport ...@@ -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 * PROGRAM: gshhs.c
* AUTHOR: Paul Wessel (pwessel@hawaii.edu) * AUTHOR: Paul Wessel (pwessel@hawaii.edu)
...@@ -378,6 +378,8 @@ void import_gshhs(FILE * fp, GeoEarthImport & geo_import) ...@@ -378,6 +378,8 @@ void import_gshhs(FILE * fp, GeoEarthImport & geo_import)
p.x = swabi4((unsigned int)p.x); p.x = swabi4((unsigned int)p.x);
p.y = swabi4((unsigned int)p.y); 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( geo_import.add_point(SPoint3(
M_PI / 180*((h.greenwich && p.x > max_east) ? p.x * 1.0e-6 - 360.0 M_PI / 180*((h.greenwich && p.x > max_east) ? p.x * 1.0e-6 - 360.0
: p.x * 1.0e-6), : p.x * 1.0e-6),
......
// $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 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
// //
...@@ -59,8 +59,7 @@ void GMSH_TransformLatLonPlugin::getInfos(char *author, char *copyright, ...@@ -59,8 +59,7 @@ void GMSH_TransformLatLonPlugin::getInfos(char *author, char *copyright,
"If `SwapOrientation' is set, the orientation of the\n" "If `SwapOrientation' is set, the orientation of the\n"
"elements is reversed. If `iView' < 0, the plugin\n" "elements is reversed. If `iView' < 0, the plugin\n"
"is run on the current view.\n" "is run on the current view.\n"
"\n" "\n" "Plugin(TransformLatLon) is executed in-place.\n");
"Plugin(TransformLatLon) is executed in-place.\n");
} }
int GMSH_TransformLatLonPlugin::getNbOptions() const int GMSH_TransformLatLonPlugin::getNbOptions() const
...@@ -79,65 +78,98 @@ void GMSH_TransformLatLonPlugin::catchErrorMessage(char *errorMessage) const ...@@ -79,65 +78,98 @@ void GMSH_TransformLatLonPlugin::catchErrorMessage(char *errorMessage) const
} }
// TransformLatLonation // TransformLatLonation
static void vector_transform(double *vin, double *vout, double x,
static void transform(double v[3], double *x, double *y, double *z,double rotate) double y, double z)
{ {
double r = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); double theta = y * M_PI / 180;
*y=asin(v[2]/r)*180/M_PI; double phi = x * M_PI / 180;
*x=atan2(v[1],v[0])*180/M_PI+rotate; vout[0] = -sin(phi) * vin[0] + cos(phi) * vin[1];
while(*x>180)*x-=360; vout[1] =
*z=0; -sin(theta) * (cos(phi) * vin[0] + sin(phi) * vin[1]) +
} cos(theta) * vin[2];
vout[2] =
static void transform_list(List_T *list, int nbList, List_T *list2, int *nbList2, cos(theta) * (cos(phi) * vin[0] + sin(phi) * vin[1]) +
int nbVert, int nbComp,double rotate) sin(theta) * vin[2];
}
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; int nb = List_Nbr(list) / nbList;
double *x2=new double[nbVert]; double *ll_out=new double[nbVert*3];
double *y2=new double[nbVert]; double *lon=ll_out;
double *z2=new double[nbVert]; 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) { for(int i = 0; i < List_Nbr(list); i += nb) {
double *x = (double *)List_Pointer_Fast(list, i); double *x = (double *)List_Pointer_Fast(list, i);
double *y = (double *)List_Pointer_Fast(list, i + nbVert); double *y = (double *)List_Pointer_Fast(list, i + nbVert);
double *z = (double *)List_Pointer_Fast(list, i + 2 * nbVert); double *z = (double *)List_Pointer_Fast(list, i + 2 * nbVert);
for(int j = 0; j < nbVert; j++) { for(int j = 0; j < nbVert; j++) {
double v[3] = { x[j], y[j], z[j] }; r[j] = sqrt(x[j]*x[j] + y[j]*y[j]+ z[j]*z[j]);
transform(v, &x2[j], &y2[j], &z2[j],rotate); 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 // check the periodicity
bool cross = false; bool cross = false;
for(int j=1;j<nbVert;j++) cross|=(fabs(x2[j]-x2[0])>180); for(int j = 1; j < nbVert; j++)
if(cross) for(int j=0;j<nbVert;j++) if(x2[j]<0)x2[j]+=360; cross |= (fabs(lon[j] - lon[0]) > 180);
for(int j = 0; j < nbVert; j++)List_Add(list2,&x2[j]); if(cross)
for(int j = 0; j < nbVert; j++)List_Add(list2,&y2[j]); for(int j = 0; j < nbVert; j++)
for(int j = 0; j < nbVert; j++)List_Add(list2,&z2[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 k = 0; k < nbVert; k++) {
double d; for(int l = 0; l < nbComp; l++)
for(int l = 0; l < nbComp; l++){ List_Read(list, i + 3 * nbVert + nbComp * k + l, &din[l]);
List_Read(list, i + 3 * nbVert + nbComp * k + l, &d); if(nbComp == 1)
List_Add(list2, &d); 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)++; (*nbList2)++;
if(cross) { if(cross) {
for(int j=0;j<nbVert;j++) x2[j]-=360; for(int j = 0; j < nbVert; j++)
for(int j = 0; j < nbVert; j++)List_Add(list2,&x2[j]); lon[j] -= 360;
for(int j = 0; j < nbVert; j++)List_Add(list2,&y2[j]); for(int j = 0; j < nbVert*3; j++)
for(int j = 0; j < nbVert; j++)List_Add(list2,&z2[j]); List_Add(list2, &ll_out[j]);
for(int k = 0; k < nbVert; k++) { for(int k = 0; k < nbVert; k++) {
double d; for(int l = 0; l < nbComp; l++)
for(int l = 0; l < nbComp; l++){ List_Read(list, i + 3 * nbVert + nbComp * k + l, &din[l]);
List_Read(list, i + 3 * nbVert + nbComp * k + l, &d); if(nbComp == 1){
List_Add(list2, &d); 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)++; (*nbList2)++;
} }
} }
delete []x2; delete[]dout;
delete []y2; delete[]din;
delete []z2; delete[]ll_out;
} }
PView *GMSH_TransformLatLonPlugin::execute(PView * v) PView *GMSH_TransformLatLonPlugin::execute(PView * v)
...@@ -146,45 +178,72 @@ PView *GMSH_TransformLatLonPlugin::execute(PView *v) ...@@ -146,45 +178,72 @@ PView *GMSH_TransformLatLonPlugin::execute(PView *v)
double rotate = (double)TransformLatLonOptions_Number[1].def; double rotate = (double)TransformLatLonOptions_Number[1].def;
PView *v1 = getView(iView, v); PView *v1 = getView(iView, v);
if(!v1) return v; if(!v1)
return v;
PView *v2 = new PView(true); PView *v2 = new PView(true);
PViewDataList *data2 = getDataList(v2); PViewDataList *data2 = getDataList(v2);
if(!data2) return v; if(!data2)
return v;
PViewDataList *data1 = getDataList(v1); PViewDataList *data1 = getDataList(v1);
if(!data1) return v; if(!data1)
return v;
// points // points
transform_list(data1->SP, data1->NbSP, data2->SP, &data2->NbSP, 1, 1, rotate); transform_list(data1->SP, data1->NbSP, data2->SP, &data2->NbSP, 1, 1,
transform_list(data1->VP, data1->NbVP, data2->VP, &data2->NbVP, 1, 3, rotate); rotate);
transform_list(data1->TP, data1->NbTP, data2->TP, &data2->NbTP, 1, 9, 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 // lines
transform_list(data1->SL, data1->NbSL, data2->SL, &data2->NbSL, 2, 1, rotate); transform_list(data1->SL, data1->NbSL, data2->SL, &data2->NbSL, 2, 1,
transform_list(data1->VL, data1->NbVL, data2->VL, &data2->NbVL, 2, 3, rotate); rotate);
transform_list(data1->TL, data1->NbTL, data2->TL, &data2->NbTL, 2, 9, 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 // triangles
transform_list(data1->ST, data1->NbST, data2->ST, &data2->NbST, 3, 1, rotate); transform_list(data1->ST, data1->NbST, data2->ST, &data2->NbST, 3, 1,
transform_list(data1->VT, data1->NbVT, data2->VT, &data2->NbVT, 3, 3, rotate); rotate);
transform_list(data1->TT, data1->NbTT, data2->TT, &data2->NbTT, 3, 9, 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 // quadrangles
transform_list(data1->SQ, data1->NbSQ, data2->SQ, &data2->NbSQ, 4, 1, rotate); transform_list(data1->SQ, data1->NbSQ, data2->SQ, &data2->NbSQ, 4, 1,
transform_list(data1->VQ, data1->NbVQ, data2->VQ, &data2->NbVQ, 4, 3, rotate); rotate);
transform_list(data1->TQ, data1->NbTQ, data2->TQ, &data2->NbTQ, 4, 9, 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 // tets
transform_list(data1->SS, data1->NbSS, data2->SS, &data2->NbSS, 4, 1, rotate); transform_list(data1->SS, data1->NbSS, data2->SS, &data2->NbSS, 4, 1,
transform_list(data1->VS, data1->NbVS, data2->VS, &data2->NbVS, 4, 3, rotate); rotate);
transform_list(data1->TS, data1->NbTS, data2->TS, &data2->NbTS, 4, 9, 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 // hexas
transform_list(data1->SH, data1->NbSH, data2->SH, &data2->NbSH, 8, 1, rotate); transform_list(data1->SH, data1->NbSH, data2->SH, &data2->NbSH, 8, 1,
transform_list(data1->VH, data1->NbVH, data2->VH, &data2->NbVH, 8, 3, rotate); rotate);
transform_list(data1->TH, data1->NbTH, data2->TH, &data2->NbTH, 8, 9, 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 // prisms
transform_list(data1->SI, data1->NbSI, data2->SI, &data2->NbSI, 6, 1, rotate); transform_list(data1->SI, data1->NbSI, data2->SI, &data2->NbSI, 6, 1,
transform_list(data1->VI, data1->NbVI, data2->VI, &data2->NbVI, 6, 3, rotate); rotate);
transform_list(data1->TI, data1->NbTI, data2->TI, &data2->NbTI, 6, 9, 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 // pyramids
transform_list(data1->SY, data1->NbSY, data2->SY, &data2->NbSY, 5, 1, rotate); transform_list(data1->SY, data1->NbSY, data2->SY, &data2->NbSY, 5, 1,
transform_list(data1->VY, data1->NbVY, data2->VY, &data2->NbVY, 5, 3, rotate); rotate);
transform_list(data1->TY, data1->NbTY, data2->TY, &data2->NbTY, 5, 9, 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->setName(data1->getName() + "_LatLon");
data2->setFileName(data1->getName() + "_LatLon.pos"); data2->setFileName(data1->getName() + "_LatLon.pos");
......
// 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] = Box;
Field[1].VIn = 2000; Field[1].VIn = 2000;
Field[1].VOut = 100000; Field[1].VOut = 100000;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment