slimFields.cpp 3.46 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
#ifdef HAVE_GMSH
lambrechts's avatar
lambrechts committed
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#include "GModel.h"
#include "Field.h"
#include "slimFields.h"

#ifdef HAVE_PROJ
#include "proj_api.h"
class utmField : public Field
{
  int iField, zone;
  projPJ pj_utm, pj_latlon;
  public:
  std::string getDescription()
  {
    return "Evaluate Field[IField] in Universal Transverse Mercator coordinates.\n";
  }
  utmField()
  {
    iField = 1;
    zone = 0;
    options["IField"] = new FieldOptionInt
      (iField, "Index of the field to evaluate");
    options["Zone"] = new FieldOptionInt
      (zone, "Zone of the UTM projection", & update_needed);
    pj_utm = pj_latlon = NULL;
  }
  const char *getName()
  {
    return "utm";
  }
  virtual ~utmField()
  {
    if(pj_utm) pj_free(pj_utm);
    if(pj_latlon) pj_free(pj_latlon);
  }
  double operator() (double x, double y, double z, GEntity *ge=0)
  {
    Field *field = GModel::current()->getFields()->get(iField);
    if(!field || iField == id) return MAX_LC;
    double lon = x;
    double lat = y;
    if (update_needed){
      if(pj_utm) pj_free(pj_utm);
      if(pj_latlon) pj_free(pj_latlon);
      std::ostringstream convertor;
      convertor << "+proj=utm +ellps=WGS84 +zone=" << zone << " +south";
      const char *utmSyst = convertor.str().c_str();
      if ( !(pj_utm = pj_init_plus(utmSyst)) ) {
        printf("Error init UTM");
        exit(1); }
      if ( !(pj_latlon = pj_init_plus("+proj=latlong +ellps=WGS84")) )  {
        printf("Error init latLon");
        exit(1); }
      update_needed = false;
    }
    int er=0;
    if ( (er = pj_transform( pj_utm, pj_latlon, 1, 1, &lon, &lat, NULL )) != 0) {
      printf ("Transform failed: %s", pj_strerrno (er));
      exit(1); }
    return (*field)(lon, lat, 0);
  }
};
#endif

class xy2d2LonLatField : public Field
{
  int iField;
  double R, phiP, thetaP;
 public:
  std::string getDescription()
  {
    return "Evaluate Field[IField] in geographic coordinates (longitude, latitude) from x and y projected on a plane\n\n";
  }
  xy2d2LonLatField()
  {
    iField = 1;
    options["IField"] = new FieldOptionInt
      (iField, "Index of the field to evaluate.");
    R = 6371e3;
    phiP = 0;
    thetaP = 0;

    options["Radius"] = new FieldOptionDouble
      (R, "radius of the sphere");
    options["Phi"] = new FieldOptionDouble
      (phiP, "longitude of the projection point (in degrees)");
    options["Theta"] = new FieldOptionDouble
      (thetaP, "latitude of the projection point (in degrees)");
  }
  const char *getName()
  {
    return "xy2d2LonLat";
  }
  double operator() (double x2d, double y2d, double z2d, GEntity *ge=0)
  {
    Field *field = GModel::current()->getFields()->get(iField);
    if(!field || iField == id) return MAX_LC;
    double phi = phiP*M_PI/180;
    double theta = thetaP*M_PI/180;
    double pOx = cos(theta)*cos(phi)*R;
    double pOy = cos(theta)*sin(phi)*R;
    double pOz = sin(theta)*R;
    double pPhiX = -sin(phi);
    double pPhiY = cos(phi);
    double pPhiZ = 0;
    double pThetaX = -sin(theta)*cos(phi);
    double pThetaY = -sin(theta)*sin(phi);
    double pThetaZ = cos(theta);
    
    double x = pPhiX * x2d + pThetaX * y2d + pOx;
		double y = pPhiY * x2d + pThetaY * y2d + pOy;
		double z = pPhiZ * x2d + pThetaZ * y2d + pOz;
    return (*field)(atan2(y, x), asin(z / R), 0);
  }
};

void slimRegisterField(GModel *m)
{
  FieldManager &fields = *m->getFields();
#ifdef HAVE_PROJ
  fields.map_type_name["utm"] = new FieldFactoryT<utmField>();
#endif
  fields.map_type_name["xy2d2LonLat"] = new FieldFactoryT<xy2d2LonLatField>();
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
125
#endif