From 2a7a5cc60c4ea81a088ac2c7c7295e879a494ed0 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Thu, 20 Mar 2008 10:21:15 +0000
Subject: [PATCH] indentation
---
Mesh/Field.cpp | 1527 +++++++++++++++++++++----------------
Mesh/Field.h | 123 +--
Plugin/GSHHS.cpp | 375 ++++++---
Plugin/GSHHS.h | 8 +-
Plugin/GeoEarthImport.cpp | 153 ++--
Plugin/GeoEarthImport.h | 38 +-
6 files changed, 1348 insertions(+), 876 deletions(-)
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 32555f7e87..99c5ea5fe3 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1,4 +1,4 @@
-// $Id: Field.cpp,v 1.21 2008-03-19 17:26:48 geuzaine Exp $
+// $Id: Field.cpp,v 1.22 2008-03-20 10:21:15 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -44,87 +44,161 @@
extern Context_T CTX;
-class FieldOptionDouble:public FieldOption{
- public:
- double &val;
- FieldOptionType get_type(){return FIELD_OPTION_DOUBLE;};
- FieldOptionDouble(double &_val,bool *_status=NULL):FieldOption(_status),val(_val){};
- double numerical_value()const { return val;};
- void numerical_value(double v){modified(); val=v;};
- void get_text_representation(std::string &v_str){
- std::ostringstream sstream;
- sstream.precision(16);
- sstream<<val;
- v_str=sstream.str();
- }
+class FieldOptionDouble:public FieldOption
+{
+public:
+ double &val;
+ FieldOptionType get_type()
+ {
+ return FIELD_OPTION_DOUBLE;
+ };
+FieldOptionDouble(double &_val, bool * _status = NULL):FieldOption(_status),
+ val(_val) {
+ };
+ double numerical_value() const
+ {
+ return val;
+ };
+ void numerical_value(double v)
+ {
+ modified();
+ val = v;
+ };
+ void get_text_representation(std::string & v_str)
+ {
+ std::ostringstream sstream;
+ sstream.precision(16);
+ sstream << val;
+ v_str = sstream.str();
+ }
};
-class FieldOptionInt:public FieldOption{
- public:
- int &val;
- FieldOptionType get_type(){return FIELD_OPTION_INT;};
- FieldOptionInt(int &_val, bool *_status=NULL):FieldOption(_status),val(_val){};
- double numerical_value()const{ return val;};
- void numerical_value(double v){modified(); val=v;};
- void get_text_representation(std::string &v_str){
- std::ostringstream sstream;
- sstream<<val;
- v_str=sstream.str();
- }
+class FieldOptionInt:public FieldOption
+{
+public:
+ int &val;
+ FieldOptionType get_type()
+ {
+ return FIELD_OPTION_INT;
+ };
+FieldOptionInt(int &_val, bool * _status = NULL):FieldOption(_status),
+ val(_val) {
+ };
+ double numerical_value() const
+ {
+ return val;
+ };
+ void numerical_value(double v)
+ {
+ modified();
+ val = v;
+ };
+ void get_text_representation(std::string & v_str)
+ {
+ std::ostringstream sstream;
+ sstream << val;
+ v_str = sstream.str();
+ }
};
-class FieldOptionList:public FieldOption{
- public:
- std::list<int> &val;
- FieldOptionType get_type(){return FIELD_OPTION_LIST;};
- FieldOptionList(std::list<int> &_val, bool *_status=NULL):FieldOption(_status),val(_val){};
- std::list<int> &list(){modified(); return val; }
- const std::list<int> &list()const { return val; }
- void get_text_representation(std::string &v_str){
- std::ostringstream sstream;
- sstream<<"{";
- for(std::list<int>::iterator it=val.begin();it!=val.end();it++){
- if(it!=val.begin())
- sstream<<", ";
- sstream<<*it;
- }
- sstream<<"}";
- v_str=sstream.str();
- }
+class FieldOptionList:public FieldOption
+{
+public:
+ std::list < int >&val;
+ FieldOptionType get_type()
+ {
+ return FIELD_OPTION_LIST;
+ };
+FieldOptionList(std::list < int >&_val, bool * _status = NULL):FieldOption(_status),
+ val(_val)
+ {
+ };
+ std::list < int >&list()
+ {
+ modified();
+ return val;
+ }
+ const std::list < int >&list() const
+ {
+ return val;
+ }
+ void get_text_representation(std::string & v_str)
+ {
+ std::ostringstream sstream;
+ sstream << "{";
+ for(std::list < int >::iterator it = val.begin(); it != val.end(); it++) {
+ if(it != val.begin())
+ sstream << ", ";
+ sstream << *it;
+ }
+ sstream << "}";
+ v_str = sstream.str();
+ }
};
-class FieldOptionString:public FieldOption{
- public:
- std::string &val;
- virtual FieldOptionType get_type(){return FIELD_OPTION_STRING;};
- FieldOptionString(std::string &_val, bool *_status=NULL):FieldOption(_status),val(_val){};
- std::string &string(){modified(); return val; }
- const std::string &string()const {return val; }
- void get_text_representation(std::string &v_str){
- std::ostringstream sstream;
- sstream<<"\""<<val<<"\"";
- v_str=sstream.str();
- }
+class FieldOptionString:public FieldOption
+{
+public:
+ std::string & val;
+ virtual FieldOptionType get_type()
+ {
+ return FIELD_OPTION_STRING;
+ };
+FieldOptionString(std::string & _val, bool * _status = NULL):FieldOption(_status),
+ val(_val) {
+ };
+ std::string & string() {
+ modified();
+ return val;
+ }
+ const std::string & string() const
+ {
+ return val;
+ }
+ void get_text_representation(std::string & v_str)
+ {
+ std::ostringstream sstream;
+ sstream << "\"" << val << "\"";
+ v_str = sstream.str();
+ }
};
-class FieldOptionBool:public FieldOption{
- public:
- bool &val;
- FieldOptionType get_type(){return FIELD_OPTION_BOOL;};
- FieldOptionBool(bool &_val, bool *_status=NULL):FieldOption(_status),val(_val){};
- double numerical_value()const { return val;};
- void numerical_value(double v){modified(); val=v;};
- void get_text_representation(std::string &v_str){
- std::ostringstream sstream;
- sstream<<val;
- v_str=sstream.str();
- }
+class FieldOptionBool:public FieldOption
+{
+public:
+ bool & val;
+ FieldOptionType get_type()
+ {
+ return FIELD_OPTION_BOOL;
+ };
+FieldOptionBool(bool & _val, bool * _status = NULL):FieldOption(_status),
+ val(_val) {
+ };
+ double numerical_value() const
+ {
+ return val;
+ };
+ void numerical_value(double v)
+ {
+ modified();
+ val = v;
+ };
+ void get_text_representation(std::string & v_str)
+ {
+ std::ostringstream sstream;
+ sstream << val;
+ v_str = sstream.str();
+ }
};
-class FieldOptionPath:public FieldOptionString{
- public:
- FieldOptionType get_type(){return FIELD_OPTION_PATH;};
+class FieldOptionPath:public FieldOptionString
+{
+public:
+ FieldOptionType get_type()
+ {
+ return FIELD_OPTION_PATH;
+ };
};
void FieldManager::reset()
{
- for(std::map<int,Field*>::iterator it = begin(); it != end(); it++){
+ for(std::map < int, Field * >::iterator it = begin(); it != end(); it++) {
delete it->second;
}
clear();
@@ -133,7 +207,7 @@ void FieldManager::reset()
Field *FieldManager::get(int id)
{
iterator it = find(id);
- if(it == end()){
+ if(it == end()) {
return NULL;
}
return it->second;
@@ -141,546 +215,686 @@ Field *FieldManager::get(int id)
Field *FieldManager::new_field(int id, const char *type_name)
{
- if(find(id) != end()){
- Msg(GERROR, "Field id %i is already defined.",id);
- return NULL;
- }
- if(map_type_name.find(type_name) == map_type_name.end()){
- Msg(GERROR, "Unknown field type \"%s\".",type_name);
- return NULL;
- }
- Field *f=(*map_type_name[type_name])();
- if(!f)
- return NULL;
- f->id=id;
- (*this)[id]=f;
+ if(find(id) != end()) {
+ Msg(GERROR, "Field id %i is already defined.", id);
+ return NULL;
+ }
+ if(map_type_name.find(type_name) == map_type_name.end()) {
+ Msg(GERROR, "Unknown field type \"%s\".", type_name);
+ return NULL;
+ }
+ Field *f = (*map_type_name[type_name]) ();
+ if(!f)
+ return NULL;
+ f->id = id;
+ (*this)[id] = f;
return f;
}
-int FieldManager::new_id(){
- int i=0;
- iterator it=begin();
- while(1){
- i++;
- while(it!=end() && it->first<i)it++;
- if(it==end() || it->first!=i)break;
- }
- return std::max(i,1);
+int FieldManager::new_id()
+{
+ int i = 0;
+ iterator it = begin();
+ while(1) {
+ i++;
+ while(it != end() && it->first < i)
+ it++;
+ if(it == end() || it->first != i)
+ break;
+ }
+ return std::max(i, 1);
}
-int FieldManager::max_id(){
- if(!empty())
- return rbegin()->first;
- else return 0;
+
+int FieldManager::max_id()
+{
+ if(!empty())
+ return rbegin()->first;
+ else
+ return 0;
}
-void FieldManager::delete_field(int id){
- iterator it=find(id);
- if(it==end()){
- Msg(GERROR, "Cannot delete field id %i, it does not exist.",id);
- return;
- }
- delete it->second;
- erase(it);
+void FieldManager::delete_field(int id)
+{
+ iterator it = find(id);
+ if(it == end()) {
+ Msg(GERROR, "Cannot delete field id %i, it does not exist.", id);
+ return;
+ }
+ delete it->second;
+ erase(it);
}
// StructuredField
-class StructuredField : public Field{
+class StructuredField:public Field
+{
double o[3], d[3];
int n[3];
double *data;
- bool error_status;
- std::string file_name;
-public :
- StructuredField(){
- options["FileName"]=new FieldOptionString(file_name,&update_needed);
- data=NULL;
- }
- const char *get_name(){
- return "Structured";
- }
- virtual ~StructuredField(){
- if(data) delete []data;
- }
- double operator()(double x, double y, double z){
- if(update_needed){
- error_status=false;
- try{
- std::ifstream input(file_name.c_str());
- if(!input.is_open())throw(1);
- input.exceptions ( std::ifstream::eofbit | std::ifstream::failbit | std::ifstream::badbit );
- input.read((char*)o, 3 * sizeof(double));
- input.read((char*)d, 3 * sizeof(double));
- input.read((char*)n, 3 * sizeof(int));
- int nt = n[0] * n[1] * n[2];
- if(data) delete []data;
- data = new double[nt];
- input.read((char*)data, nt * sizeof(double));
- input.close();
- }
- catch(...){
- error_status=true;
- Msg(GERROR,"Field %i : error reading file %s",this->id,file_name.c_str());
- }
- update_needed=false;
- }
- if(error_status)return MAX_LC;
- //tri-linear
- int id[2][3];
- double xi[3];
- double xyz[3]={x,y,z};
- for(int i = 0; i < 3; i++){
- id[0][i] = (int)floor((xyz[i] - o[i]) / d[i]);
- id[1][i] = id[0][i] + 1;
- id[0][i] = std::max(std::min(id[0][i], n[i] - 1), 0);
- id[1][i] = std::max(std::min(id[1][i], n[i] - 1), 0);
- xi[i] = xyz[i] - (o[i] + id[0][i] * d[i]);
- xi[i] = std::max(std::min(xi[i], 1.), 0.);
- }
- double v = 0;
- for(int i = 0; i < 2; i++)
- for(int j = 0;j < 2; j++)
- for(int k = 0; k < 2; k++){
- v += data[id[i][0] * n[1] * n[2] + id[j][1] * n[2] + id[k][2]]
- * (i * xi[0] + (1 - i) * (1 - xi[0]))
- * (j * xi[1] + (1 - j) * (1 - xi[1]))
- * (k * xi[2] + (1 - k) * (1 - xi[2]));
- }
- return v;
- }
- FieldDialogBox *&dialog_box(){
- static FieldDialogBox *dialogBox=NULL;
- return dialogBox;
- }
+ bool error_status;
+ std::string file_name;
+public:StructuredField()
+ {
+ options["FileName"] = new FieldOptionString(file_name, &update_needed);
+ data = NULL;
+ }
+ const char *get_name()
+ {
+ return "Structured";
+ }
+ virtual ~ StructuredField() {
+ if(data)
+ delete[]data;
+ }
+ double operator() (double x, double y, double z)
+ {
+ if(update_needed) {
+ error_status = false;
+ try {
+ std::ifstream input(file_name.c_str());
+ if(!input.is_open())
+ throw(1);
+ input.
+ exceptions(std::ifstream::eofbit | std::ifstream::failbit | std::
+ ifstream::badbit);
+ input.read((char *)o, 3 * sizeof(double));
+ input.read((char *)d, 3 * sizeof(double));
+ input.read((char *)n, 3 * sizeof(int));
+ int nt = n[0] * n[1] * n[2];
+ if(data)
+ delete[]data;
+ data = new double[nt];
+ input.read((char *)data, nt * sizeof(double));
+ input.close();
+ }
+ catch(...) {
+ error_status = true;
+ Msg(GERROR, "Field %i : error reading file %s", this->id,
+ file_name.c_str());
+ }
+ update_needed = false;
+ }
+ if(error_status)
+ return MAX_LC;
+ //tri-linear
+ int id[2][3];
+ double xi[3];
+ double xyz[3] = { x, y, z };
+ for(int i = 0; i < 3; i++) {
+ id[0][i] = (int)floor((xyz[i] - o[i]) / d[i]);
+ id[1][i] = id[0][i] + 1;
+ id[0][i] = std::max(std::min(id[0][i], n[i] - 1), 0);
+ id[1][i] = std::max(std::min(id[1][i], n[i] - 1), 0);
+ xi[i] = (xyz[i] - (o[i] + id[0][i] * d[i])) / d[i];
+ xi[i] = std::max(std::min(xi[i], 1.), 0.);
+ }
+ double v = 0;
+ for(int i = 0; i < 2; i++)
+ for(int j = 0; j < 2; j++)
+ for(int k = 0; k < 2; k++) {
+ v += data[id[i][0] * n[1] * n[2] + id[j][1] * n[2] + id[k][2]]
+ * (i * xi[0] + (1 - i) * (1 - xi[0]))
+ * (j * xi[1] + (1 - j) * (1 - xi[1]))
+ * (k * xi[2] + (1 - k) * (1 - xi[2]));
+ }
+ return v;
+ }
+ FieldDialogBox *&dialog_box()
+ {
+ static FieldDialogBox *dialogBox = NULL;
+ return dialogBox;
+ }
};
-class LonLatField : public Field{
- int field_id;
+class LonLatField:public Field
+{
+ int field_id;
public:
- LonLatField(){
- field_id=1;
- options["IField"]=new FieldOptionInt(field_id);
- }
- const char *get_name(){return "LonLat";}
- 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);
- }
- FieldDialogBox *&dialog_box(){
- static FieldDialogBox *dialogBox=NULL;
- return dialogBox;
- }
+ LonLatField()
+ {
+ field_id = 1;
+ options["IField"] = new FieldOptionInt(field_id);
+ }
+ const char *get_name()
+ {
+ return "LonLat";
+ }
+ 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);
+ }
+ FieldDialogBox *&dialog_box()
+ {
+ static FieldDialogBox *dialogBox = NULL;
+ return dialogBox;
+ }
};
-class BoxField : public Field{
- double v_in,v_out,x_min,x_max,y_min,y_max,z_min,z_max;
- public:
- BoxField(){
- v_in=v_out=x_min=x_max=y_min=y_max=z_min=z_max=0;
- options["VIn"]=new FieldOptionDouble(v_in);
- options["VOut"]=new FieldOptionDouble(v_out);
- options["XMin"]=new FieldOptionDouble(x_min);
- options["XMax"]=new FieldOptionDouble(x_max);
- options["YMin"]=new FieldOptionDouble(y_min);
- options["YMax"]=new FieldOptionDouble(y_max);
- options["ZMin"]=new FieldOptionDouble(z_min);
- options["ZMax"]=new FieldOptionDouble(z_max);
- }
- const char *get_name(){return "Box";}
- FieldDialogBox *&dialog_box(){
- static FieldDialogBox *dialogBox=NULL;
- return dialogBox;
- }
- double operator()(double x, double y, double z){
- return (x<=x_max && x>=x_min && y<=y_max && y>=y_min && z<=z_max && z>=z_min) ?
- v_in : v_out;
- }
+class BoxField:public Field
+{
+ double v_in, v_out, x_min, x_max, y_min, y_max, z_min, z_max;
+public:
+ BoxField()
+ {
+ v_in = v_out = x_min = x_max = y_min = y_max = z_min = z_max = 0;
+ options["VIn"] = new FieldOptionDouble(v_in);
+ options["VOut"] = new FieldOptionDouble(v_out);
+ options["XMin"] = new FieldOptionDouble(x_min);
+ options["XMax"] = new FieldOptionDouble(x_max);
+ options["YMin"] = new FieldOptionDouble(y_min);
+ options["YMax"] = new FieldOptionDouble(y_max);
+ options["ZMin"] = new FieldOptionDouble(z_min);
+ options["ZMax"] = new FieldOptionDouble(z_max);
+ }
+ const char *get_name()
+ {
+ return "Box";
+ }
+ FieldDialogBox *&dialog_box()
+ {
+ static FieldDialogBox *dialogBox = NULL;
+ return dialogBox;
+ }
+ double operator() (double x, double y, double z)
+ {
+ return (x <= x_max && x >= x_min && y <= y_max && y >= y_min && z <= z_max
+ && z >= z_min) ? v_in : v_out;
+ }
};
-class ThresholdField : public Field{
- int iField;
+class ThresholdField:public Field
+{
+ int iField;
double dmin, dmax, lcmin, lcmax;
- public:
- const char *get_name(){
- return "Threshold";
- }
- ThresholdField(){
- iField=0;
- dmin=1;
- dmax=10;
- lcmin=0.1;
- lcmax=1;
- options["IField"]=new FieldOptionInt(iField);
- options["DistMin"]=new FieldOptionDouble(dmin);
- options["DistMax"]=new FieldOptionDouble(dmax);
- options["LcMin"]=new FieldOptionDouble(lcmin);
- options["LcMax"]=new FieldOptionDouble(lcmax);
- }
- double operator()(double x, double y, double z){
- Field *field=GModel::current()->getFields()->get(iField);
- double r = ((*field)(x, y, z) - dmin) / (dmax - dmin);
- r = std::max(std::min(r, 1.), 0.);
- double lc = lcmin * (1 - r) + lcmax * r;
- return lc;
- }
- FieldDialogBox *&dialog_box(){
- static FieldDialogBox *dialogBox=NULL;
- return dialogBox;
- }
+public:
+ const char *get_name()
+ {
+ return "Threshold";
+ }
+ ThresholdField()
+ {
+ iField = 0;
+ dmin = 1;
+ dmax = 10;
+ lcmin = 0.1;
+ lcmax = 1;
+ options["IField"] = new FieldOptionInt(iField);
+ options["DistMin"] = new FieldOptionDouble(dmin);
+ options["DistMax"] = new FieldOptionDouble(dmax);
+ options["LcMin"] = new FieldOptionDouble(lcmin);
+ options["LcMax"] = new FieldOptionDouble(lcmax);
+ }
+ double operator() (double x, double y, double z)
+ {
+ Field *field = GModel::current()->getFields()->get(iField);
+ double r = ((*field) (x, y, z) - dmin) / (dmax - dmin);
+ r = std::max(std::min(r, 1.), 0.);
+ double lc = lcmin * (1 - r) + lcmax * r;
+ return lc;
+ }
+ FieldDialogBox *&dialog_box()
+ {
+ static FieldDialogBox *dialogBox = NULL;
+ return dialogBox;
+ }
};
-class GradientField : public Field{
- int iField,kind;
- double delta;
- public:
- const char *get_name(){
- return "Gradient";
- }
- GradientField():iField(0),kind(3),delta(CTX.lc/1e4){
- options["IField"]=new FieldOptionInt(iField);
- options["Kind"]=new FieldOptionInt(kind);
- options["Delta"]=new FieldOptionDouble(delta);
- }
- double operator()(double x, double y, double z){
- Field *field=GModel::current()->getFields()->get(iField);
- double gx, gy, gz;
- switch(kind){
- case 0 : /* x */
- return ((*field)(x + delta / 2, y, z) - (*field)(x - delta / 2, y, z)) / delta;
- case 1 : /* y */
- return ((*field)(x, y + delta / 2, z) - (*field)(x, y - delta / 2, z)) / delta;
- case 2 : /* z */
- return ((*field)(x, y, z + delta / 2) - (*field)(x, y, z - delta / 2)) / delta;
- case 3 : /* norm */
- gx = ((*field)(x + delta / 2, y, z) - (*field)(x - delta / 2, y, z)) / delta;
- gy = ((*field)(x, y + delta / 2, z) - (*field)(x, y - delta / 2, z)) / delta;
- gz = ((*field)(x, y, z + delta / 2) - (*field)(x, y, z - delta / 2)) / delta;
- return sqrt(gx * gx + gy * gy + gz * gz);
- default :
- Msg(GERROR, "Field %i : Unknown kind (%i) of gradient.",this->id, kind);
- return MAX_LC;
- }
- }
- FieldDialogBox *&dialog_box(){
- static FieldDialogBox *dialogBox=NULL;
- return dialogBox;
- }
+class GradientField:public Field
+{
+ int iField, kind;
+ double delta;
+public:
+ const char *get_name()
+ {
+ return "Gradient";
+ }
+ GradientField():iField(0), kind(3), delta(CTX.lc / 1e4)
+ {
+ options["IField"] = new FieldOptionInt(iField);
+ options["Kind"] = new FieldOptionInt(kind);
+ options["Delta"] = new FieldOptionDouble(delta);
+ }
+ double operator() (double x, double y, double z)
+ {
+ Field *field = GModel::current()->getFields()->get(iField);
+ double gx, gy, gz;
+ switch (kind) {
+ case 0: /* x */
+ return ((*field) (x + delta / 2, y, z) -
+ (*field) (x - delta / 2, y, z)) / delta;
+ case 1: /* y */
+ return ((*field) (x, y + delta / 2, z) -
+ (*field) (x, y - delta / 2, z)) / delta;
+ case 2: /* z */
+ return ((*field) (x, y, z + delta / 2) -
+ (*field) (x, y, z - delta / 2)) / delta;
+ case 3: /* norm */
+ gx =
+ ((*field) (x + delta / 2, y, z) -
+ (*field) (x - delta / 2, y, z)) / delta;
+ gy =
+ ((*field) (x, y + delta / 2, z) -
+ (*field) (x, y - delta / 2, z)) / delta;
+ gz =
+ ((*field) (x, y, z + delta / 2) -
+ (*field) (x, y, z - delta / 2)) / delta;
+ return sqrt(gx * gx + gy * gy + gz * gz);
+ default:
+ Msg(GERROR, "Field %i : Unknown kind (%i) of gradient.", this->id,
+ kind);
+ return MAX_LC;
+ }
+ }
+ FieldDialogBox *&dialog_box()
+ {
+ static FieldDialogBox *dialogBox = NULL;
+ return dialogBox;
+ }
};
#if defined(HAVE_MATH_EVAL)
-class MathEvalExpression{
- bool error_status;
- std::list<Field*> *list;
- int nvalues;
+class MathEvalExpression
+{
+ bool error_status;
+ std::list < Field * >*list;
+ int nvalues;
char **names;
double *values;
void *eval;
- int *evaluators_id;
- std::string function;
- char *c_str_function;
- public:
- double evaluate(double x,double y, double z){
- if(error_status) return MAX_LC;
- for(int i=0;i<nvalues;i++){
- Field *f;
- switch (evaluators_id[i]){
- case -1: values[i]=x; break;
- case -2: values[i]=y; break;
- case -3: values[i]=z; break;
- default: {
- f=GModel::current()->getFields()->get(evaluators_id[i]);
- values[i] = f ? (*f)(x,y,z) : MAX_LC;
- }
- }
- }
- return evaluator_evaluate(eval, nvalues, names, values);
- }
- MathEvalExpression(){
- eval=NULL;
- values=NULL;
- c_str_function=NULL;
- evaluators_id=NULL;
- }
- bool set_function(const std::string &f){
- free_members();
- error_status=false;
- c_str_function=strdup(f.c_str());
- eval = evaluator_create(c_str_function);
- if(!eval){
- error_status=true;
- return false;
- }
- evaluator_get_variables(eval,&names,&nvalues);
- values = new double[nvalues];
- evaluators_id=new int[nvalues];
- for(int i=0;i<nvalues;i++){
- int id;
- if(!strcmp("x",names[i])) evaluators_id[i]=-1;
- else if(!strcmp("y",names[i])) evaluators_id[i]=-2;
- else if(!strcmp("z",names[i])) evaluators_id[i]=-3;
- else if(sscanf(names[i],"F%i",&id)==1) evaluators_id[i]=id;
- else{
- Msg(GERROR, "Unknown matheval argument \"%s\"\n",names[i]);
- error_status=true;
- return false;
- }
- }
- return true;
- }
- void free_members(){
- if(c_str_function)free(c_str_function);
- if(eval) evaluator_destroy(eval);
- if(values)delete[]values;
- if(evaluators_id) delete evaluators_id;
- }
- ~MathEvalExpression(){
- free_members();
- }
+ int *evaluators_id;
+ std::string function;
+ char *c_str_function;
+public:
+ double evaluate(double x, double y, double z)
+ {
+ if(error_status)
+ return MAX_LC;
+ for(int i = 0; i < nvalues; i++)
+ {
+ Field *f;
+ switch (evaluators_id[i]) {
+ case -1:
+ values[i] = x;
+ break;
+ case -2:values[i] = y;
+ break;
+ case -3:values[i] = z;
+ break;
+ default:
+ {
+ f = GModel::current()->getFields()->get(evaluators_id[i]);
+ values[i] = f ? (*f) (x, y, z) : MAX_LC;
+ }
+ }
+ }
+ return evaluator_evaluate(eval, nvalues, names, values);
+ }
+ MathEvalExpression() {
+ eval = NULL;
+ values = NULL;
+ c_str_function = NULL;
+ evaluators_id = NULL;
+ }
+ bool set_function(const std::string & f)
+ {
+ free_members();
+ error_status = false;
+ c_str_function = strdup(f.c_str());
+ eval = evaluator_create(c_str_function);
+ if(!eval) {
+ error_status = true;
+ return false;
+ }
+ evaluator_get_variables(eval, &names, &nvalues);
+ values = new double[nvalues];
+ evaluators_id = new int[nvalues];
+ for(int i = 0; i < nvalues; i++) {
+ int id;
+ if(!strcmp("x", names[i]))
+ evaluators_id[i] = -1;
+ else if(!strcmp("y", names[i]))
+ evaluators_id[i] = -2;
+ else if(!strcmp("z", names[i]))
+ evaluators_id[i] = -3;
+ else if(sscanf(names[i], "F%i", &id) == 1)
+ evaluators_id[i] = id;
+ else {
+ Msg(GERROR, "Unknown matheval argument \"%s\"\n", names[i]);
+ error_status = true;
+ return false;
+ }
+ }
+ return true;
+ }
+ void free_members()
+ {
+ if(c_str_function)
+ free(c_str_function);
+ if(eval)
+ evaluator_destroy(eval);
+ if(values)
+ delete[]values;
+ if(evaluators_id)
+ delete evaluators_id;
+ }
+ ~MathEvalExpression() {
+ free_members();
+ }
};
-class MathEvalField : public Field{
- MathEvalExpression expr;
- std::string f;
- public:
- MathEvalField(){
- options["F"]=new FieldOptionString(f,&update_needed);
- }
- double operator()(double x, double y, double z){
- if(update_needed){
- if(!expr.set_function(f))
- Msg(GERROR, "Field %i : Invalid matheval expression \"%s\"\n",this->id,f.c_str());
- update_needed=false;
- }
- return expr.evaluate(x,y,z);
- }
- FieldDialogBox *&dialog_box(){
- static FieldDialogBox *dialogBox=NULL;
- return dialogBox;
- }
- const char *get_name(){return "MathEval";}
+class MathEvalField:public Field
+{
+ MathEvalExpression expr;
+ std::string f;
+public:
+ MathEvalField()
+ {
+ options["F"] = new FieldOptionString(f, &update_needed);
+ }
+ double operator() (double x, double y, double z)
+ {
+ if(update_needed) {
+ if(!expr.set_function(f))
+ Msg(GERROR, "Field %i : Invalid matheval expression \"%s\"\n",
+ this->id, f.c_str());
+ update_needed = false;
+ }
+ return expr.evaluate(x, y, z);
+ }
+ FieldDialogBox *&dialog_box()
+ {
+ static FieldDialogBox *dialogBox = NULL;
+ return dialogBox;
+ }
+ const char *get_name()
+ {
+ return "MathEval";
+ }
};
-class ParametricField : public Field{
- MathEvalExpression expr[3];
- std::string f[3];
- int ifield;
- public:
- ParametricField(){
- options["IField"] = new FieldOptionInt(ifield);
- options["FX"]=new FieldOptionString(f[0],&update_needed);
- options["FY"]=new FieldOptionString(f[1],&update_needed);
- options["FZ"]=new FieldOptionString(f[2],&update_needed);
- }
- double operator()(double x, double y, double z){
- if(update_needed){
- for(int i=0;i<3;i++){
- if(!expr[i].set_function(f[i]))
- Msg(GERROR, "Field %i : Invalid matheval expression \"%s\"\n",this->id,f[i].c_str());
- }
- update_needed=false;
- }
- return (*GModel::current()->getFields()->get(ifield))(expr[0].evaluate(x,y,z),expr[1].evaluate(x,y,z),expr[2].evaluate(x,y,z));
- }
- FieldDialogBox *&dialog_box(){
- static FieldDialogBox *dialogBox=NULL;
- return dialogBox;
- }
- const char *get_name(){return "Param";}
+class ParametricField:public Field
+{
+ MathEvalExpression expr[3];
+ std::string f[3];
+ int ifield;
+public:
+ ParametricField()
+ {
+ options["IField"] = new FieldOptionInt(ifield);
+ options["FX"] = new FieldOptionString(f[0], &update_needed);
+ options["FY"] = new FieldOptionString(f[1], &update_needed);
+ options["FZ"] = new FieldOptionString(f[2], &update_needed);
+ }
+ double operator() (double x, double y, double z)
+ {
+ if(update_needed) {
+ for(int i = 0; i < 3; i++) {
+ if(!expr[i].set_function(f[i]))
+ Msg(GERROR, "Field %i : Invalid matheval expression \"%s\"\n",
+ this->id, f[i].c_str());
+ }
+ update_needed = false;
+ }
+ return (*GModel::current()->getFields()->get(ifield)) (expr[0].
+ evaluate(x, y, z),
+ expr[1].evaluate(x,
+ y,
+ z),
+ expr[2].evaluate(x,
+ y,
+ z));
+ }
+ FieldDialogBox *&dialog_box()
+ {
+ static FieldDialogBox *dialogBox = NULL;
+ return dialogBox;
+ }
+ const char *get_name()
+ {
+ return "Param";
+ }
};
#endif
-class PostViewField:public Field{
+class PostViewField:public Field
+{
OctreePost *octree;
- public :
- int view_index;
- double operator()(double x, double y, double z) {
- // FIXME: should test unique view num instead, but that would be slower
- if(view_index < 0 || view_index >= (int)PView::list.size()) return MAX_LC;
- if(update_needed){
- if(octree)delete octree;
- octree = new OctreePost(PView::list[view_index]);
- update_needed=false;
- }
- double l = 0.;
- if(!octree->searchScalar(x, y, z, &l, 0)){
- // uncomment the following to try really hard to find an element
- // around the point
- /*
- double fact[9] = {0.001, 0.0025, 0.005, 0.0075, 0.01, 0.025, 0.05, 0.075, 0.1};
- for(int i = 0; i < 9; i++){
- double eps = CTX.lc * fact[i];
- if(octree->searchScalar(x + eps, y, z, &l, 0)) break;
- if(octree->searchScalar(x - eps, y, z, &l, 0)) break;
- if(octree->searchScalar(x, y + eps, z, &l, 0)) break;
- if(octree->searchScalar(x, y - eps, z, &l, 0)) break;
- if(octree->searchScalar(x, y, z + eps, &l, 0)) break;
- if(octree->searchScalar(x, y, z - eps, &l, 0)) break;
- if(octree->searchScalar(x + eps, y - eps, z - eps, &l, 0)) break;
- if(octree->searchScalar(x + eps, y + eps, z - eps, &l, 0)) break;
- if(octree->searchScalar(x - eps, y - eps, z - eps, &l, 0)) break;
- if(octree->searchScalar(x - eps, y + eps, z - eps, &l, 0)) break;
- if(octree->searchScalar(x + eps, y - eps, z + eps, &l, 0)) break;
- if(octree->searchScalar(x + eps, y + eps, z + eps, &l, 0)) break;
- if(octree->searchScalar(x - eps, y - eps, z + eps, &l, 0)) break;
- if(octree->searchScalar(x - eps, y + eps, z + eps, &l, 0)) break;
- }
- */
- }
- //if(l <= 0) return MAX_LC;
- return l;
- }
- const char *get_name(){ return "PostView"; }
- PostViewField() {
- octree=NULL;
- options["IView"]=new FieldOptionInt(view_index,&update_needed);
- }
- ~PostViewField(){
- if(octree) delete octree;
- }
- FieldDialogBox *&dialog_box(){
- static FieldDialogBox *dialogBox=NULL;
- return dialogBox;
- }
+public:int view_index;
+ double operator() (double x, double y, double z)
+ {
+ // FIXME: should test unique view num instead, but that would be slower
+ if(view_index < 0 || view_index >= (int)PView::list.size())
+ return MAX_LC;
+ if(update_needed)
+ {
+ if(octree)
+ delete octree;
+ octree = new OctreePost(PView::list[view_index]);
+ update_needed = false;
+ }
+ double l = 0.;
+ if(!octree->searchScalar(x, y, z, &l, 0)) {
+ // uncomment the following to try really hard to find an element
+ // around the point
+ /*
+ double fact[9] = {0.001, 0.0025, 0.005, 0.0075, 0.01, 0.025, 0.05, 0.075, 0.1};
+ for(int i = 0; i < 9; i++){
+ double eps = CTX.lc * fact[i];
+ if(octree->searchScalar(x + eps, y, z, &l, 0)) break;
+ if(octree->searchScalar(x - eps, y, z, &l, 0)) break;
+ if(octree->searchScalar(x, y + eps, z, &l, 0)) break;
+ if(octree->searchScalar(x, y - eps, z, &l, 0)) break;
+ if(octree->searchScalar(x, y, z + eps, &l, 0)) break;
+ if(octree->searchScalar(x, y, z - eps, &l, 0)) break;
+ if(octree->searchScalar(x + eps, y - eps, z - eps, &l, 0)) break;
+ if(octree->searchScalar(x + eps, y + eps, z - eps, &l, 0)) break;
+ if(octree->searchScalar(x - eps, y - eps, z - eps, &l, 0)) break;
+ if(octree->searchScalar(x - eps, y + eps, z - eps, &l, 0)) break;
+ if(octree->searchScalar(x + eps, y - eps, z + eps, &l, 0)) break;
+ if(octree->searchScalar(x + eps, y + eps, z + eps, &l, 0)) break;
+ if(octree->searchScalar(x - eps, y - eps, z + eps, &l, 0)) break;
+ if(octree->searchScalar(x - eps, y + eps, z + eps, &l, 0)) break;
+ }
+ */
+ }
+ //if(l <= 0) return MAX_LC;
+ return l;
+ }
+ const char *get_name()
+ {
+ return "PostView";
+ }
+ PostViewField() {
+ octree = NULL;
+ options["IView"] = new FieldOptionInt(view_index, &update_needed);
+ }
+ ~PostViewField() {
+ if(octree)
+ delete octree;
+ }
+ FieldDialogBox *&dialog_box()
+ {
+ static FieldDialogBox *dialogBox = NULL;
+ return dialogBox;
+ }
};
-class MinField : public Field{
- std::list<int> idlist;
- public:
- MinField(){
- options["FieldsList"]=new FieldOptionList(idlist,&update_needed);
- }
- double operator()(double x, double y, double z){
- double v=MAX_LC;
- for(std::list<int>::iterator it=idlist.begin();it!=idlist.end();it++){
- Field *f=(GModel::current()->getFields()->get(*it));
- if(f)
- v=std::min(v,(*f)(x,y,z));
- }
- return v;
- }
- FieldDialogBox *&dialog_box(){
- static FieldDialogBox *dialogBox=NULL;
- return dialogBox;
- }
- const char *get_name(){return "Min";}
+class MinField:public Field
+{
+ std::list < int >idlist;
+public:
+ MinField()
+ {
+ options["FieldsList"] = new FieldOptionList(idlist, &update_needed);
+ }
+ double operator() (double x, double y, double z)
+ {
+ double v = MAX_LC;
+ for(std::list < int >::iterator it = idlist.begin(); it != idlist.end();
+ it++) {
+ Field *f = (GModel::current()->getFields()->get(*it));
+ if(f)
+ v = std::min(v, (*f) (x, y, z));
+ }
+ return v;
+ }
+ FieldDialogBox *&dialog_box()
+ {
+ static FieldDialogBox *dialogBox = NULL;
+ return dialogBox;
+ }
+ const char *get_name()
+ {
+ return "Min";
+ }
};
-class MaxField : public Field{
- std::list<int> idlist;
- public:
- MaxField(){
- options["FieldsList"]=new FieldOptionList(idlist,&update_needed);
- }
- double operator()(double x, double y, double z){
- double v=-MAX_LC;
- for(std::list<int>::iterator it=idlist.begin();it!=idlist.end();it++){
- Field *f=(GModel::current()->getFields()->get(*it));
- if(f)
- v=std::max(v,(*f)(x,y,z));
- }
- return v;
- }
- FieldDialogBox *&dialog_box(){
- static FieldDialogBox *dialogBox=NULL;
- return dialogBox;
- }
- const char *get_name(){return "Max";}
+class MaxField:public Field
+{
+ std::list < int >idlist;
+public:
+ MaxField()
+ {
+ options["FieldsList"] = new FieldOptionList(idlist, &update_needed);
+ }
+ double operator() (double x, double y, double z)
+ {
+ double v = -MAX_LC;
+ for(std::list < int >::iterator it = idlist.begin(); it != idlist.end();
+ it++) {
+ Field *f = (GModel::current()->getFields()->get(*it));
+ if(f)
+ v = std::max(v, (*f) (x, y, z));
+ }
+ return v;
+ }
+ FieldDialogBox *&dialog_box()
+ {
+ static FieldDialogBox *dialogBox = NULL;
+ return dialogBox;
+ }
+ const char *get_name()
+ {
+ return "Max";
+ }
};
+
#ifdef HAVE_ANN
-class AttractorField : public Field{
- ANNkd_tree* kdtree;
+class AttractorField:public Field
+{
+ ANNkd_tree *kdtree;
ANNpointArray zeronodes;
ANNidxArray index;
ANNdistArray dist;
- std::list<int> nodes_id;
- std::list<int> edges_id;
- int n_nodes_by_edge;
-public :
- AttractorField(): kdtree (0), zeronodes(0) {
- index = new ANNidx[1];
- dist = new ANNdist[1];
- options["NodesList"]=new FieldOptionList(nodes_id,&update_needed);
- options["EdgesList"]=new FieldOptionList(edges_id,&update_needed);
- options["NNodesByEdge"]=new FieldOptionInt(n_nodes_by_edge,&update_needed);
- n_nodes_by_edge=20;
- }
- ~AttractorField(){
- if(kdtree) delete kdtree;
- if(zeronodes) annDeallocPts(zeronodes);
- delete [] index;
- delete [] dist;
- }
- const char *get_name(){
- return "Attractor";
- }
- virtual double operator()(double X, double Y, double Z){
- if(update_needed){
- if(zeronodes){
- annDeallocPts(zeronodes);
- delete kdtree;
- }
- int totpoints = nodes_id.size()+n_nodes_by_edge*edges_id.size();
- if (totpoints)
- zeronodes = annAllocPts(totpoints, 4);
- int k = 0;
- for(std::list <int>::iterator it = nodes_id.begin(); it != nodes_id.end(); ++it){
- Vertex *v = FindPoint(*it);
- if(v){
- zeronodes[k][0]=v->Pos.X;
- zeronodes[k][1]=v->Pos.Y;
- zeronodes[k++][2]=v->Pos.Z;
- }else{
- GVertex *gv = GModel::current()->getVertexByTag(*it);
- if(gv){
- zeronodes[k][0]=gv->x();
- zeronodes[k][1]=gv->y();
- zeronodes[k++][2]=gv->z();
- }
- }
- }
- for(std::list <int>::iterator it = edges_id.begin(); it != edges_id.end(); ++it){
- Curve *c = FindCurve(*it);
- if(c){
- for(int i = 0; i < n_nodes_by_edge ; i++){
- double u = (double)i / (n_nodes_by_edge - 1);
- Vertex V = InterpolateCurve(c, u, 0);
- zeronodes[k][0]=V.Pos.X;
- zeronodes[k][1]=V.Pos.Y;
- zeronodes[k++][2]=V.Pos.Z;
- }
- }else{
- GEdge *ge = GModel::current()->getEdgeByTag(*it);
- if(ge){
- for(int i = 0; i < n_nodes_by_edge; i++){
- double u = (double)i / (n_nodes_by_edge - 1);
- Range<double> b = ge->parBounds(0);
- double t = b.low() + u * (b.high() - b.low());
- GPoint gp = ge->point(t);
- zeronodes[k][0]=gp.x();
- zeronodes[k][1]=gp.y();
- zeronodes[k++][2]=gp.z();
- }
- }
- }
- }
- kdtree = new ANNkd_tree(zeronodes, totpoints, 3);
- update_needed=false;
- }
- double xyz[3] = {X, Y, Z};
- kdtree->annkSearch(xyz, 1, index, dist);
- return sqrt(dist[0]);
- }
- FieldDialogBox *&dialog_box(){
- static FieldDialogBox *dialogBox=NULL;
- return dialogBox;
- }
+ std::list < int >nodes_id;
+ std::list < int >edges_id;
+ int n_nodes_by_edge;
+public:AttractorField():kdtree(0), zeronodes(0)
+ {
+ index = new ANNidx[1];
+ dist = new ANNdist[1];
+ options["NodesList"] = new FieldOptionList(nodes_id, &update_needed);
+ options["EdgesList"] = new FieldOptionList(edges_id, &update_needed);
+ options["NNodesByEdge"] =
+ new FieldOptionInt(n_nodes_by_edge, &update_needed);
+ n_nodes_by_edge = 20;
+ }
+ ~AttractorField()
+ {
+ if(kdtree)
+ delete kdtree;
+ if(zeronodes)
+ annDeallocPts(zeronodes);
+ delete[]index;
+ delete[]dist;
+ }
+ const char *get_name()
+ {
+ return "Attractor";
+ }
+ virtual double operator() (double X, double Y, double Z)
+ {
+ if(update_needed) {
+ if(zeronodes) {
+ annDeallocPts(zeronodes);
+ delete kdtree;
+ }
+ int totpoints = nodes_id.size() + n_nodes_by_edge * edges_id.size();
+ if(totpoints)
+ zeronodes = annAllocPts(totpoints, 4);
+ int k = 0;
+ for(std::list < int >::iterator it = nodes_id.begin();
+ it != nodes_id.end(); ++it) {
+ Vertex *v = FindPoint(*it);
+ if(v) {
+ zeronodes[k][0] = v->Pos.X;
+ zeronodes[k][1] = v->Pos.Y;
+ zeronodes[k++][2] = v->Pos.Z;
+ }
+ else {
+ GVertex *gv = GModel::current()->getVertexByTag(*it);
+ if(gv) {
+ zeronodes[k][0] = gv->x();
+ zeronodes[k][1] = gv->y();
+ zeronodes[k++][2] = gv->z();
+ }
+ }
+ }
+ for(std::list < int >::iterator it = edges_id.begin();
+ it != edges_id.end(); ++it) {
+ Curve *c = FindCurve(*it);
+ if(c) {
+ for(int i = 0; i < n_nodes_by_edge; i++) {
+ double u = (double)i / (n_nodes_by_edge - 1);
+ Vertex V = InterpolateCurve(c, u, 0);
+ zeronodes[k][0] = V.Pos.X;
+ zeronodes[k][1] = V.Pos.Y;
+ zeronodes[k++][2] = V.Pos.Z;
+ }
+ }
+ else {
+ GEdge *ge = GModel::current()->getEdgeByTag(*it);
+ if(ge) {
+ for(int i = 0; i < n_nodes_by_edge; i++) {
+ double u = (double)i / (n_nodes_by_edge - 1);
+ Range < double >b = ge->parBounds(0);
+ double t = b.low() + u * (b.high() - b.low());
+ GPoint gp = ge->point(t);
+ zeronodes[k][0] = gp.x();
+ zeronodes[k][1] = gp.y();
+ zeronodes[k++][2] = gp.z();
+ }
+ }
+ }
+ }
+ kdtree = new ANNkd_tree(zeronodes, totpoints, 3);
+ update_needed = false;
+ }
+ double xyz[3] = { X, Y, Z };
+ kdtree->annkSearch(xyz, 1, index, dist);
+ return sqrt(dist[0]);
+ }
+ FieldDialogBox *&dialog_box()
+ {
+ static FieldDialogBox *dialogBox = NULL;
+ return dialogBox;
+ }
};
#endif
#if 0
-void addMapLc (std::map<MVertex*, double> &maplc, MVertex *v, double l)
+void addMapLc(std::map < MVertex *, double >&maplc, MVertex * v, double l)
{
- std::map<MVertex*, double>::iterator it = maplc.find(v);
- if(it == maplc.end()) maplc[v] = l;
- else if(it->second > l) it->second = l;
+ std::map < MVertex *, double >::iterator it = maplc.find(v);
+ if(it == maplc.end())
+ maplc[v] = l;
+ else if(it->second > l)
+ it->second = l;
}
#endif
@@ -697,83 +911,83 @@ public:
virtual double operator()(double X, double Y, double Z) ;
virtual void eval(double X, double Y, double Z, double &l, double &lpt, double &dist) ;
};*/
-AttractorField_1DMesh::AttractorField_1DMesh(GModel *m, double dmax, double dmin,
- double lcmax)
- : _dmax(dmax), _dmin(dmin), _lcmax(lcmax)
+AttractorField_1DMesh::AttractorField_1DMesh(GModel * m, double dmax,
+ double dmin, double lcmax)
+:_dmax(dmax), _dmin(dmin), _lcmax(lcmax)
{
GModel::eiter it = m->firstEdge();
- std::map<MVertex*, double> maplc;
+ std::map < MVertex *, double >maplc;
- while(it != m->lastEdge()){
+ while(it != m->lastEdge()) {
MVertex *first = (*it)->getBeginVertex()->mesh_vertices[0];
- for(unsigned int i = 1; i <= (*it)->mesh_vertices.size(); ++i){
- MVertex *last = (i == (*it)->mesh_vertices.size()) ?
- (*it)->getEndVertex()->mesh_vertices[0] : (*it)->mesh_vertices[i];
+ for(unsigned int i = 1; i <= (*it)->mesh_vertices.size(); ++i) {
+ MVertex *last = (i == (*it)->mesh_vertices.size())?
+ (*it)->getEndVertex()->mesh_vertices[0] : (*it)->mesh_vertices[i];
double l = sqrt((first->x() - last->x()) * (first->x() - last->x()) +
- (first->y() - last->y()) * (first->y() - last->y()) +
- (first->z() - last->z()) * (first->z() - last->z()));
+ (first->y() - last->y()) * (first->y() - last->y()) +
+ (first->z() - last->z()) * (first->z() - last->z()));
addMapLc(maplc, first, l);
addMapLc(maplc, last, l);
first = last;
}
- }
-
- std::map<MVertex*, double>::iterator itm = maplc.begin();
-
- while(itm != maplc.end()){
+ }
+
+ std::map < MVertex *, double >::iterator itm = maplc.begin();
+
+ while(itm != maplc.end()) {
addPoint(itm->first->x(), itm->first->y(), itm->first->z());
lcs.push_back(itm->second);
}
}
-AttractorField_1DMesh::AttractorField_1DMesh(GFace *gf, double dmax, double dmin,
- double lcmax)
- : _dmax(dmax), _dmin(dmin), _lcmax(lcmax)
+AttractorField_1DMesh::AttractorField_1DMesh(GFace * gf, double dmax,
+ double dmin, double lcmax)
+:_dmax(dmax), _dmin(dmin), _lcmax(lcmax)
{
- std::list<GEdge*> edges = gf->edges();
- std::list<GEdge*>::iterator it = edges.begin();
+ std::list < GEdge * >edges = gf->edges();
+ std::list < GEdge * >::iterator it = edges.begin();
- std::map<MVertex*, double> maplc;
- std::map<MVertex*, double> maplc2;
+ std::map < MVertex *, double >maplc;
+ std::map < MVertex *, double >maplc2;
- while(it != edges.end()){
+ while(it != edges.end()) {
MVertex *first = (*it)->getBeginVertex()->mesh_vertices[0];
- for(unsigned int i = 1; i <= (*it)->mesh_vertices.size(); ++i){
- MVertex *last = (i == (*it)->mesh_vertices.size()) ?
- (*it)->getEndVertex()->mesh_vertices[0] : (*it)->mesh_vertices[i];
+ for(unsigned int i = 1; i <= (*it)->mesh_vertices.size(); ++i) {
+ MVertex *last = (i == (*it)->mesh_vertices.size())?
+ (*it)->getEndVertex()->mesh_vertices[0] : (*it)->mesh_vertices[i];
double l = sqrt((first->x() - last->x()) * (first->x() - last->x()) +
- (first->y() - last->y()) * (first->y() - last->y()) +
- (first->z() - last->z()) * (first->z() - last->z()));
+ (first->y() - last->y()) * (first->y() - last->y()) +
+ (first->z() - last->z()) * (first->z() - last->z()));
double t;
- last->getParameter(0,t);
- double l2 = LC_MVertex_PNTS(last->onWhat(),t,0);
+ last->getParameter(0, t);
+ double l2 = LC_MVertex_PNTS(last->onWhat(), t, 0);
addMapLc(maplc, first, l);
addMapLc(maplc, last, l);
addMapLc(maplc2, first, l2);
addMapLc(maplc2, last, l2);
- first = last;
+ first = last;
}
++it;
- }
-
- std::map<MVertex*, double>::iterator itm = maplc.begin();
- while(itm != maplc.end()){
+ }
+
+ std::map < MVertex *, double >::iterator itm = maplc.begin();
+ while(itm != maplc.end()) {
addPoint(itm->first->x(), itm->first->y(), itm->first->z());
lcs.push_back(itm->second);
++itm;
}
- itm = maplc2.begin();
- while(itm != maplc2.end()){
+ itm = maplc2.begin();
+ while(itm != maplc2.end()) {
lcs2.push_back(itm->second);
++itm;
}
}
-double AttractorField_1DMesh::operator()(double X, double Y, double Z)
+double AttractorField_1DMesh::operator() (double X, double Y, double Z)
{
#ifdef HAVE_ANN
- double xyz[3] = {X, Y, Z};
+ double xyz[3] = { X, Y, Z };
kdtree->annkSearch(xyz, maxpts, index, dist);
double d = sqrt(dist[0]);
double lcmin = lcs[index[0]];
@@ -782,57 +996,65 @@ double AttractorField_1DMesh::operator()(double X, double Y, double Z)
double lc = lcmin * (1 - r) + _lcmax * r;
return lc;
#else
- Msg(GERROR,"GMSH should be compiled with ANN in order to enable attractors");
+ Msg(GERROR,
+ "GMSH should be compiled with ANN in order to enable attractors");
return 1.e22;
#endif
}
-void AttractorField_1DMesh::eval(double X, double Y, double Z, double &lcmin, double &lcpt, double &d)
+void AttractorField_1DMesh::eval(double X, double Y, double Z, double &lcmin,
+ double &lcpt, double &d)
{
#ifdef HAVE_ANN
- double xyz[3] = {X, Y, Z};
+ double xyz[3] = { X, Y, Z };
kdtree->annkSearch(xyz, maxpts, index, dist);
d = sqrt(dist[0]);
lcmin = lcs[index[0]];
- lcpt = lcs2[index[0]] * CTX.mesh.lc_factor;
+ lcpt = lcs2[index[0]] * CTX.mesh.lc_factor;
#else
- Msg(GERROR,"GMSH should be compiled with ANN in order to enable attractors");
+ Msg(GERROR,
+ "GMSH should be compiled with ANN in order to enable attractors");
#endif
}
#endif
-template<class F>
-class FieldFactoryT:public FieldFactory{
- public:
- Field *operator()(){return new F;};
+template < class F > class FieldFactoryT:public FieldFactory {
+public:
+ Field * operator()() {
+ return new F;
+ };
};
-template<class F>
-Field *field_factory(){return new F();};
+template < class F > Field * field_factory()
+{
+ return new F();
+};
-FieldManager::FieldManager(){
- map_type_name["Structured"]=new FieldFactoryT<StructuredField>();
- map_type_name["Threshold"]=new FieldFactoryT<ThresholdField>();
- map_type_name["Box"]=new FieldFactoryT<BoxField>();
- map_type_name["LonLat"]=new FieldFactoryT<LonLatField>();
+FieldManager::FieldManager()
+{
+ map_type_name["Structured"] = new FieldFactoryT < StructuredField > ();
+ map_type_name["Threshold"] = new FieldFactoryT < ThresholdField > ();
+ map_type_name["Box"] = new FieldFactoryT < BoxField > ();
+ map_type_name["LonLat"] = new FieldFactoryT < LonLatField > ();
#if defined(HAVE_MATH_EVAL)
- map_type_name["Param"]=new FieldFactoryT<ParametricField>();
- map_type_name["MathEval"]=new FieldFactoryT<MathEvalField>();
+ map_type_name["Param"] = new FieldFactoryT < ParametricField > ();
+ map_type_name["MathEval"] = new FieldFactoryT < MathEvalField > ();
#endif
#if defined(HAVE_ANN)
- map_type_name["Attractor"]=new FieldFactoryT<AttractorField>();
+ map_type_name["Attractor"] = new FieldFactoryT < AttractorField > ();
#endif
- map_type_name["PostView"]=new FieldFactoryT<PostViewField>();
- map_type_name["Gradient"]=new FieldFactoryT<GradientField>();
- map_type_name["Min"]=new FieldFactoryT<MinField>();
- map_type_name["Max"]=new FieldFactoryT<MaxField>();
- background_field=-1;
+ map_type_name["PostView"] = new FieldFactoryT < PostViewField > ();
+ map_type_name["Gradient"] = new FieldFactoryT < GradientField > ();
+ map_type_name["Min"] = new FieldFactoryT < MinField > ();
+ map_type_name["Max"] = new FieldFactoryT < MaxField > ();
+ background_field = -1;
}
-static void evaluate(Field *field, List_T *list1, int nbElm1, int nbNod,
- int nbComp, int comp )
+static void evaluate(Field * field, List_T * list1, int nbElm1, int nbNod,
+ int nbComp, int comp)
{
- if(!nbElm1) return;
+ if(!nbElm1)
+ return;
int nb = List_Nbr(list1) / nbElm1;
for(int i = 0; i < List_Nbr(list1); i += nb) {
double *x = (double *)List_Pointer_Fast(list1, i);
@@ -840,31 +1062,34 @@ static void evaluate(Field *field, List_T *list1, int nbElm1, int nbNod,
double *z = (double *)List_Pointer_Fast(list1, i + 2 * nbNod);
for(int j = 0; j < nbNod; j++) {
// store data from the main view into v
- double *val1 = (double *)List_Pointer_Fast(list1,
- i + 3 * nbNod +
- nbNod * nbComp * 0 + nbComp * j);
- val1[comp] = (*field)(x[j],y[j],z[j]);
+ double *val1 = (double *)List_Pointer_Fast(list1,
+ i + 3 * nbNod +
+ nbNod * nbComp * 0 +
+ nbComp * j);
+ val1[comp] = (*field) (x[j], y[j], z[j]);
}
}
}
-Field::Field(){
+Field::Field()
+{
}
-void Field::put_on_view(PView *view,int comp)
+void Field::put_on_view(PView * view, int comp)
{
- PViewDataList *data = dynamic_cast<PViewDataList*>(view->getData());
- if(!data) return ;
+ PViewDataList *data = dynamic_cast < PViewDataList * >(view->getData());
+ if(!data)
+ return;
evaluate(this, data->SP, data->NbSP, 1, 1, 0);
evaluate(this, data->SL, data->NbSL, 2, 1, 0);
evaluate(this, data->ST, data->NbST, 3, 1, 0);
- evaluate(this, data->SQ, data->NbSQ, 4, 1, 0);
+ evaluate(this, data->SQ, data->NbSQ, 4, 1, 0);
evaluate(this, data->SS, data->NbSS, 4, 1, 0);
- evaluate(this, data->SH, data->NbSH, 8, 1, 0);
+ evaluate(this, data->SH, data->NbSH, 8, 1, 0);
evaluate(this, data->SI, data->NbSI, 6, 1, 0);
evaluate(this, data->SY, data->NbSY, 5, 1, 0);
- for(int cc = 0; cc < 3; cc++){
- if(comp < 0 || comp == cc){
+ for(int cc = 0; cc < 3; cc++) {
+ if(comp < 0 || comp == cc) {
evaluate(this, data->VP, data->NbVP, 1, 3, cc);
evaluate(this, data->VL, data->NbVL, 2, 3, cc);
evaluate(this, data->VT, data->NbVT, 3, 3, cc);
@@ -875,11 +1100,11 @@ void Field::put_on_view(PView *view,int comp)
evaluate(this, data->VY, data->NbVY, 5, 3, cc);
}
}
- for(int cc = 0; cc < 9; cc++){
- if(comp < 0 || comp == cc){
+ for(int cc = 0; cc < 9; cc++) {
+ if(comp < 0 || comp == cc) {
evaluate(this, data->TP, data->NbTP, 1, 9, cc);
evaluate(this, data->TL, data->NbTL, 2, 9, cc);
- evaluate(this, data->TT, data->NbTT, 3, 9, cc);
+ evaluate(this, data->TT, data->NbTT, 3, 9, cc);
evaluate(this, data->TQ, data->NbTQ, 4, 9, cc);
evaluate(this, data->TS, data->NbTS, 4, 9, cc);
evaluate(this, data->TH, data->NbTH, 8, 9, cc);
@@ -891,11 +1116,11 @@ void Field::put_on_view(PView *view,int comp)
view->setChanged(true);
}
-void FieldManager::set_background_mesh(int iView){
- int id=new_id();
- Field *f=new_field(id,"PostView");
- f->options["IView"]->numerical_value(iView);
- (*this)[id]=f;
- background_field=id;
+void FieldManager::set_background_mesh(int iView)
+{
+ int id = new_id();
+ Field *f = new_field(id, "PostView");
+ f->options["IView"]->numerical_value(iView);
+ (*this)[id] = f;
+ background_field = id;
}
-
diff --git a/Mesh/Field.h b/Mesh/Field.h
index 2e3f24657b..ab26804c6a 100644
--- a/Mesh/Field.h
+++ b/Mesh/Field.h
@@ -26,59 +26,94 @@
#include "PView.h"
class Field;
-typedef enum {FIELD_OPTION_DOUBLE=0,FIELD_OPTION_INT,FIELD_OPTION_STRING,FIELD_OPTION_PATH,FIELD_OPTION_BOOL,FIELD_OPTION_LIST}FieldOptionType;
-class FieldOption{
- protected:
- bool *status;
- inline void modified(){if(status)*status=true;}
- public:
- FieldOption(bool *_status):status(_status){};
- virtual FieldOptionType get_type()=0;
- virtual void get_text_representation(std::string &v_str)=0;
- virtual void numerical_value(double val){throw (1);}
- virtual double numerical_value()const {throw (1);}
- virtual const std::list<int> & list()const {throw (1);}
- virtual std::list<int> & list(){throw (1);}
- virtual const std::string & string()const{throw (1);}
- virtual std::string & string(){throw (1);}
+typedef enum
+{ FIELD_OPTION_DOUBLE =
+ 0, FIELD_OPTION_INT, FIELD_OPTION_STRING, FIELD_OPTION_PATH,
+ FIELD_OPTION_BOOL, FIELD_OPTION_LIST } FieldOptionType;
+class FieldOption
+{
+protected:
+ bool * status;
+ inline void modified()
+ {
+ if(status)
+ *status = true;
+ }
+public:
+ FieldOption(bool * _status):status(_status)
+ {
+ };
+ virtual FieldOptionType get_type() = 0;
+ virtual void get_text_representation(std::string & v_str) = 0;
+ virtual void numerical_value(double val)
+ {
+ throw(1);
+ }
+ virtual double numerical_value() const
+ {
+ throw(1);
+ }
+ virtual const std::list < int >&list() const
+ {
+ throw(1);
+ }
+ virtual std::list < int >&list()
+ {
+ throw(1);
+ }
+ virtual const std::string & string() const
+ {
+ throw(1);
+ }
+ virtual std::string & string()
+ {
+ throw(1);
+ }
};
class FieldDialogBox;
-class Field{
- struct lstr{
- bool operator() (const char* s1, const char* s2 ) const{
- return strcmp(s1,s2)<0;
- }
- };
+class Field
+{
+ struct lstr
+ {
+ bool operator() (const char *s1, const char *s2)const
+ {
+ return strcmp(s1, s2) < 0;
+ }
+ };
public:
- int id;
- std::map<const char *, FieldOption*,lstr> options;
- virtual double operator()(double x, double y, double z) = 0;
- virtual ~Field(){}
- bool update_needed;
- Field();
- virtual const char *get_name()=0;
- virtual FieldDialogBox *&dialog_box()=0;
- void put_on_view(PView *view,int comp=-1);
+ int id;
+ std::map < const char *, FieldOption *, lstr > options;
+ virtual double operator() (double x, double y, double z) = 0;
+ virtual ~ Field()
+ {
+ }
+ bool update_needed;
+ Field();
+ virtual const char *get_name() = 0;
+ virtual FieldDialogBox *&dialog_box() = 0;
+ void put_on_view(PView * view, int comp = -1);
};
-class FieldFactory{
- public:
- virtual Field *operator()()=0;
+class FieldFactory
+{
+public:
+ virtual Field * operator() () = 0;
};
-class FieldManager:public std::map<int, Field*>{
- public:
- std::map<const std::string,FieldFactory*> map_type_name;
+class FieldManager:public std::map < int, Field * >
+{
+public:
+ std::map < const std::string, FieldFactory * >map_type_name;
void reset();
Field *get(int id);
- Field *new_field(int id, const char *type_name);
- void delete_field(int id);
- int new_id();
- int max_id();
- FieldManager();
- int background_field;
- /* compatibility with -bgm */
- void set_background_mesh(int iView);
+ Field *new_field(int id, const char *type_name);
+ void delete_field(int id);
+ int new_id();
+ int max_id();
+ FieldManager();
+ int background_field;
+ /* compatibility with -bgm */
+ void set_background_mesh(int iView);
};
#endif
diff --git a/Plugin/GSHHS.cpp b/Plugin/GSHHS.cpp
index 3d824a1ff0..26d9766ad1 100644
--- a/Plugin/GSHHS.cpp
+++ b/Plugin/GSHHS.cpp
@@ -1,4 +1,4 @@
-// $Id: GSHHS.cpp,v 1.2 2008-03-19 17:26:54 geuzaine Exp $
+// $Id: GSHHS.cpp,v 1.3 2008-03-20 10:21:15 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -29,7 +29,8 @@
#include "GeoEarthImport.h"
extern Context_T CTX;
-/* $Id: GSHHS.cpp,v 1.2 2008-03-19 17:26:54 geuzaine Exp $
+// ************** GSHHS **************
+/* $Id: GSHHS.cpp,v 1.3 2008-03-20 10:21:15 remacle Exp $
*
* PROGRAM: gshhs.c
* AUTHOR: Paul Wessel (pwessel@hawaii.edu)
@@ -57,46 +58,238 @@ extern Context_T CTX;
* GNU General Public License for more details.
*
* Contact info: www.soest.hawaii.edu/pwessel */
-
/* For byte swapping on little-endian systems (GSHHS is bigendian) */
#define swabi2(i2) (((i2) >> 8) + (((i2) & 255) << 8))
#define swabi4(i4) (((i4) >> 24) + (((i4) >> 8) & 65280) + (((i4) & 65280) << 8) + (((i4) & 255) << 24))
-#define _POSIX_SOURCE 1 /* GSHHS code is POSIX compliant */
-
+#define _POSIX_SOURCE 1 /* GSHHS code is POSIX compliant */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
-
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
-
-#ifndef SEEK_CUR /* For really ancient systems */
+#ifndef SEEK_CUR /* For really ancient systems */
#define SEEK_CUR 1
#endif
+ struct GSHHS
+{ /* Global Self-consistent Hierarchical High-resolution Shorelines */
+ int id; /* Unique polygon id number, starting at 0 */
+ int n; /* Number of points in this polygon */
+ int level; /* 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake */
+ int west, east, south, north; /* min/max extent in micro-degrees */
+ int area; /* Area of polygon in 1/10 km^2 */
+ int version; /* Version of GSHHS polygon (3 is latest and first with this item) */
+ short int greenwich; /* Greenwich is 1 if Greenwich is crossed */
+ short int source; /* 0 = CIA WDBII, 1 = WVS */
+};
-struct GSHHS { /* Global Self-consistent Hierarchical High-resolution Shorelines */
- int id; /* Unique polygon id number, starting at 0 */
- int n; /* Number of points in this polygon */
- int level; /* 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake */
- int west, east, south, north; /* min/max extent in micro-degrees */
- int area; /* Area of polygon in 1/10 km^2 */
- int version; /* Version of GSHHS polygon (3 is latest and first with this item) */
- short int greenwich; /* Greenwich is 1 if Greenwich is crossed */
- short int source; /* 0 = CIA WDBII, 1 = WVS */
+struct POINT
+{ /* Each lon, lat pair is stored in micro-degrees in 4-byte integer format */
+ int x;
+ int y;
};
-struct POINT { /* Each lon, lat pair is stored in micro-degrees in 4-byte integer format */
- int x;
- int y;
+void import_gshhs(FILE * fp, GeoEarthImport & geo_import)
+{
+ double w, e, s, n;
+ char source;
+ int k, max_east = 270000000, info, n_read, flip;
+ struct POINT p;
+ struct GSHHS h;
+ while(1) {
+ n_read = fread((void *)&h, (size_t) sizeof(struct GSHHS), (size_t) 1, fp);
+ if(n_read != 1)
+ break;
+ flip = (!(h.level > 0 && h.level < 5)); /* Take as sign that byte-swabbing is needed */
+ if(flip) {
+ h.id = swabi4((unsigned int)h.id);
+ h.n = swabi4((unsigned int)h.n);
+ h.level = swabi4((unsigned int)h.level);
+ h.west = swabi4((unsigned int)h.west);
+ h.east = swabi4((unsigned int)h.east);
+ h.south = swabi4((unsigned int)h.south);
+ h.north = swabi4((unsigned int)h.north);
+ h.area = swabi4((unsigned int)h.area);
+ h.version = swabi4((unsigned int)h.version);
+ h.greenwich = swabi2((unsigned int)h.greenwich);
+ h.source = swabi2((unsigned int)h.source);
+ }
+ w = h.west * 1.0e-6; /* Convert from microdegrees to degrees */
+ e = h.east * 1.0e-6;
+ s = h.south * 1.0e-6;
+ n = h.north * 1.0e-6;
+ source = (h.source == 1) ? 'W' : 'C'; /* Either WVS or CIA (WDBII) pedigree */
+ if(h.level != 1) { // Skip data (lake,island,pond)
+ fseek(fp, (long)(h.n * sizeof(struct POINT)), SEEK_CUR);
+ continue;
+ }
+ for(k = 0; k < h.n; k++) {
+ if(fread((void *)&p, (size_t) sizeof(struct POINT), (size_t) 1, fp) !=
+ 1) {
+ Msg(GERROR,
+ "gshhs: Error reading gshhs file for polygon %d, point %d.\n",
+ h.id, k);
+ return;
+ }
+ if(flip) {
+ p.x = swabi4((unsigned int)p.x);
+ p.y = swabi4((unsigned int)p.y);
+ }
+ SPoint2 ll(((h.greenwich
+ && p.x >
+ max_east) ? p.x * 1.0e-6 -
+ 360.0 : p.x * 1.0e-6) * M_PI / 180,
+ (p.y * 1.0e-6) * M_PI / 180);
+ geo_import.add_point_lon_lat(ll);
+ }
+ geo_import.end_loop();
+ }
+ geo_import.end_surface();
+}
+
+
+// ************** UTM **************
+#include "math.h"
+#include "stdio.h"
+class UTM
+{
+ public:double x, y;
+ int zone;
+ UTM():x(0), y(0), zone(-1)
+ {
+ };
+ UTM(double _x, double _y, int _zone):x(_x), y(_y), zone(_zone)
+ {
+ };
+};
+class UTMConverter
+{
+ double meridionalarc(double lat, double lon);
+ double a, b, n, n2, n3, n4, n5, e, e2, e1, e12, e13, e14, J1, J2, J3, J4,
+ Ap, Bp, Cp, Dp, Ep, e4, e6, ep, ep2, ep4, k0, mu_fact;
+ public: void latlon2utm(double lat, double lon, UTM & utm);
+ void utm2latlon(UTM & utm, double &lat, double &lon);
+ UTMConverter(double _a = 6378137, double _b = 6356752.3142) {
+ /* see http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.HTM */
+ a = _a; /* Rayon Equatorial */
+ b = _b; /* Rayon Polaire */
+ n = (a - b) / (a + b);
+ n2 = n * n;
+ n3 = n * n * n;
+ n4 = n * n * n * n;
+ n5 = n * n * n * n * n;
+ e = sqrt(1 - b * b / a / a);
+ e2 = e * e;
+ e1 = (1 - sqrt(1 - e2)) / (1 + sqrt(1 - e2));
+ e12 = e1 * e1;
+ e13 = e1 * e1 * e1;
+ e14 = e1 * e1 * e1 * e1;
+ J1 = (3 * e1 / 2 - 27 * e13 / 32);
+ J2 = (21 * e12 / 16 - 55 * e14 / 32);
+ J3 = 151 * e13 / 96;
+ J4 = 1097 * e14 / 512;
+ Ap = a * (1 - n + (5. / 4.) * (n2 - n3) + (81. / 64.) * (n4 - n5));
+ Bp =
+ -3 * a * n / 2 * (1 - n + (7. / 8.) * (n2 - n3) +
+ (55. / 64.) * (n4 - n5));
+ Cp = 14 * a * n2 / 16 * (1 - n + (3. / 4) * (n2 - n3));
+ Dp = -35 * a * n3 / 48 * (1 - n + 11. / 16. * (n2 - n3));
+ Ep = +315 * a * n4 / 51 * (1 - n);
+ e4 = e2 * e2;
+ e6 = e2 * e2 * e2;
+ ep = e * a / b;
+ ep2 = ep * ep;
+ ep4 = ep2 * ep2;
+ k0 = 0.9996;
+ mu_fact = 1 / (k0 * a * (1 - e2 / 4 - 3 * e4 / 64 - 5 * e6 / 256));
+ }
};
+
+double UTMConverter::meridionalarc(double lat, double lon)
+{
+ return Ap * lat + Bp * sin(2 * lat) + Cp * sin(4 * lat) +
+ Dp * sin(6 * lat) + Ep;
+}
+
+void UTMConverter::utm2latlon(UTM & utm, double &lat, double &lon)
+{
+ double mu = utm.y * mu_fact;
+ double fp =
+ mu + J1 * sin(2 * mu) + J2 * sin(4 * mu) + J3 * sin(6 * mu) +
+ J4 * sin(8 * mu);
+ double cfp = cos(fp);
+ double cfp2 = cfp * cfp;
+ double sfp = sin(fp);
+ double sfp2 = sfp * sfp;
+ double c1 = ep2 * cfp2;
+ double c12 = c1 * c1;
+ double t1 = sfp2 / cfp2;
+ double t12 = t1 * t1;
+ double r1 = a * (1 - e2) / pow((1 - e2 * sfp2), 1.5);
+ double n1 = a / sqrt(1 - e2 * sfp2);
+ double d = (utm.x - 5e5) / (n1 * k0);
+ double d2 = d * d;
+ double d3 = d2 * d;
+ double d4 = d2 * d2;
+ double d5 = d4 * d;
+ double d6 = d4 * d2;
+ lat =
+ fp - n1 * sfp / cfp / r1 * (d2 / 2 -
+ (5 + 3 * t1 + 10 * c1 - 4 * c12 -
+ 9 * ep2) * d4 / 24 + (61 + 90 * t1 +
+ 298 * c1 + 45 * t12 -
+ 3 * c12 -
+ 252 * ep2) * d6 / 720);
+ lon =
+ ((utm.zone - 0.5) / 30 - 1) * M_PI + (d - (1 + 2 * t1 + c1) * d3 / 6 +
+ (5 - 2 * c1 + 28 * t1 - 3 * c12 +
+ 8 * ep2 +
+ 24 * t12) * d5 / 120) / cfp;
+}
+
+void UTMConverter::latlon2utm(double lat, double lon, UTM & utm)
+{
+ if(utm.zone == -1)
+ utm.zone = (int)ceil((lon / M_PI + 1) * 30);
+ double S = meridionalarc(lat, lon);
+ double slat = sin(lat);
+ double clat = cos(lat);
+ double slat2 = slat * slat;
+ double clat2 = clat * clat;
+ double clat3 = clat2 * clat;
+ double clat4 = clat3 * clat;
+ double tlat2 = slat2 / clat2;
+ double nu = a / sqrt(1 - e * e * slat2);
+ double p = lon - ((utm.zone - 0.5) / 30 - 1) * M_PI;
+ double p2 = p * p;
+ double p3 = p * p2;
+ double p4 = p2 * p2;
+ utm.y =
+ S * k0 + k0 * nu * slat * clat / 2 * p2 +
+ k0 * nu * slat * clat3 / 24 * (5 - tlat2 + 9 * ep2 * clat2 +
+ 4 * ep4 * clat4) * p4;
+ utm.x =
+ k0 * nu * clat * p + (k0 * nu * clat3 / 6) * (1 - tlat2 +
+ ep2 * clat2) * p3 + 5e5;
+}
+
+// ************** MAIN PLUGIN **************
StringXNumber GSHHSOptions_Number[] = {
- {GMSH_FULLRC, "iField", NULL, -1.}
+ {GMSH_FULLRC, "iField", NULL, -1.},
+ {GMSH_FULLRC, "UTMEquatorialRadius", NULL, 6378137},
+ {GMSH_FULLRC, "UTMPolarRadius", NULL, 6356752.3142},
+ {GMSH_FULLRC, "UTMScale", NULL, 1},
+ {GMSH_FULLRC, "UTMShiftX", NULL, 0},
+ {GMSH_FULLRC, "UTMShiftY", NULL, 0},
+ {GMSH_FULLRC, "UTMZone", NULL, 0}
+
};
StringXString GSHHSOptions_String[] = {
- {GMSH_FULLRC, "InFileName",NULL,"gshhs_c.b"},
- {GMSH_FULLRC, "OutFileName",NULL,"earth.geo"}
+ {GMSH_FULLRC, "InFileName", NULL, "gshhs_c.b"},
+ {GMSH_FULLRC, "OutFileName", NULL, "earth.geo"},
+ {GMSH_FULLRC, "Format", NULL, "gshhs"},
+ {GMSH_FULLRC, "Coordinate", NULL, "cartesian"}
};
extern "C"
@@ -113,18 +306,18 @@ void GMSH_GSHHSPlugin::getName(char *name) const
}
void GMSH_GSHHSPlugin::getInfos(char *author, char *copyright,
- char *help_text) const
+ char *help_text) const
{
- strcpy(author, "J. Lambrechts (jonathanlambrechts@gmail.org)");
+ strcpy(author, "J. Lambrechts (jonathanlambrechts@gmail.com)");
strcpy(copyright, "GPL");
- strcpy(help_text,
- "Plugin(GSHHS) import GSHHS data.\n");
+ strcpy(help_text, "Plugin(GSHHS) import GSHHS data.\n");
}
int GMSH_GSHHSPlugin::getNbOptions() const
{
return sizeof(GSHHSOptions_Number) / sizeof(StringXNumber);
}
+
int GMSH_GSHHSPlugin::getNbOptionsStr() const
{
return sizeof(GSHHSOptions_String) / sizeof(StringXString);
@@ -134,6 +327,7 @@ StringXNumber *GMSH_GSHHSPlugin::getOption(int iopt)
{
return &GSHHSOptions_Number[iopt];
}
+
StringXString *GMSH_GSHHSPlugin::getOptionStr(int iopt)
{
return &GSHHSOptions_String[iopt];
@@ -144,74 +338,67 @@ void GMSH_GSHHSPlugin::catchErrorMessage(char *errorMessage) const
strcpy(errorMessage, "GSHHS failed...");
}
-PView *GMSH_GSHHSPlugin::execute(PView *v)
+PView *GMSH_GSHHSPlugin::execute(PView * v)
{
int iField = (int)GSHHSOptions_Number[0].def;
- char *filename=(char*)GSHHSOptions_String[0].def;
- char *outfilename=(char*)GSHHSOptions_String[1].def;
- Field *field=NULL;
- GeoEarthImport geo_import(outfilename);
- if(iField!=-1){
- field = GModel::current()->getFields()->get(iField);
- if(!field){
- Msg(GERROR, "Field[%d] does not exist", iField);
- return NULL;
- }else{
- geo_import.set_size_field(field);
- }
- }
- FILE *fp;
- if ((fp = fopen (filename, "rb")) == NULL ) {
- Msg(GERROR, "gshhs: Could not find file %s.\n", filename);
- return NULL;
- }
- double w, e, s, n;
- char source;
- int k, max_east = 270000000, info, n_read, flip;
- struct POINT p;
- struct GSHHS h;
- while(1){
- n_read = fread ((void *)&h, (size_t)sizeof (struct GSHHS), (size_t)1, fp);
- if(n_read!=1)break;
- flip = (! (h.level > 0 && h.level < 5)); /* Take as sign that byte-swabbing is needed */
- if (flip) {
- h.id = swabi4 ((unsigned int)h.id);
- h.n = swabi4 ((unsigned int)h.n);
- h.level = swabi4 ((unsigned int)h.level);
- h.west = swabi4 ((unsigned int)h.west);
- h.east = swabi4 ((unsigned int)h.east);
- h.south = swabi4 ((unsigned int)h.south);
- h.north = swabi4 ((unsigned int)h.north);
- h.area = swabi4 ((unsigned int)h.area);
- h.version = swabi4 ((unsigned int)h.version);
- h.greenwich = swabi2 ((unsigned int)h.greenwich);
- h.source = swabi2 ((unsigned int)h.source);
- }
- w = h.west * 1.0e-6; /* Convert from microdegrees to degrees */
- e = h.east * 1.0e-6;
- s = h.south * 1.0e-6;
- n = h.north * 1.0e-6;
- source = (h.source == 1) ? 'W' : 'C'; /* Either WVS or CIA (WDBII) pedigree */
- if ( h.level!=1) { // Skip data (lake,island,pond)
- fseek (fp, (long)(h.n * sizeof(struct POINT)), SEEK_CUR);
- continue;
- }
- for (k = 0; k < h.n; k++) {
- if (fread ((void *)&p, (size_t)sizeof(struct POINT), (size_t)1, fp) != 1) {
- Msg(GERROR, "gshhs: Error reading gshhs file for polygon %d, point %d.\n",h.id, k);
- return NULL;
- }
- if (flip) {
- p.x = swabi4 ((unsigned int)p.x);
- p.y = swabi4 ((unsigned int)p.y);
- }
- SPoint2 ll(
- ((h.greenwich && p.x > max_east) ? p.x * 1.0e-6 - 360.0 : p.x * 1.0e-6)*M_PI/180,
- (p.y * 1.0e-6)*M_PI/180);
- geo_import.add_point_lon_lat(ll);
- }
- geo_import.end_loop();
- }
- geo_import.end_surface();
+ char *filename = (char *)GSHHSOptions_String[0].def;
+ char *outfilename = (char *)GSHHSOptions_String[1].def;
+ std::string format(GSHHSOptions_String[2].def);
+ std::string coordinate_name(GSHHSOptions_String[3].def);
+ double scale = GSHHSOptions_Number[3].def;
+ double shiftx = GSHHSOptions_Number[4].def;
+ double shifty = GSHHSOptions_Number[5].def;
+ if(coordinate_name != "cartesian" && coordinate_name != "utm") {
+ Msg(GERROR, "gshhs: unknown coordinate system %s.\n",
+ coordinate_name.c_str());
+ return NULL;
+ }
+ Field *field = NULL;
+ GeoEarthImport geo_import(outfilename);
+ if(iField != -1) {
+ field = GModel::current()->getFields()->get(iField);
+ if(!field) {
+ Msg(GERROR, "Field[%d] does not exist", iField);
+ return NULL;
+ }
+ else {
+ geo_import.set_size_field(field);
+ }
+ }
+ FILE *fp;
+ if((fp = fopen(filename, "rb")) == NULL) {
+ Msg(GERROR, "gshhs: Could not open file %s.\n", filename);
+ return NULL;
+ }
+ if(format == "gshhs") {
+ import_gshhs(fp, geo_import);
+ }
+ else if(format == "loops2") {
+ int npoints_in_loop;
+ UTMConverter utmc((double)GSHHSOptions_Number[1].def,
+ (double)GSHHSOptions_Number[2].def);
+ UTM utm;
+ while(fscanf(fp, "%d", &npoints_in_loop) == 1) {
+ utm.zone = (int)GSHHSOptions_Number[6].def;
+ for(int i = 0; i < npoints_in_loop; i++) {
+ if(fscanf(fp, "%le %le", &utm.x, &utm.y) != 2) {
+ Msg(GERROR, "gshhs: Error reading loops2 file \'%s\'.\n",
+ filename);
+ return NULL;
+ }
+ double lat, lon;
+ utm.x = utm.x * scale + shiftx;
+ utm.y = utm.y * scale + shifty;
+ utmc.utm2latlon(utm, lat, lon);
+ geo_import.add_point_lon_lat(SPoint2(lon, lat));
+ }
+ geo_import.end_loop();
+ }
+ geo_import.end_surface();
+ }
+ else {
+ Msg(GERROR, "gshhs: Unknown format %s.\n", format.c_str());
+ return NULL;
+ }
return NULL;
}
diff --git a/Plugin/GSHHS.h b/Plugin/GSHHS.h
index 0191f64869..11b6145d35 100644
--- a/Plugin/GSHHS.h
+++ b/Plugin/GSHHS.h
@@ -27,16 +27,16 @@ extern "C"
GMSH_Plugin *GMSH_RegisterGSHHSPlugin();
}
-class GMSH_GSHHSPlugin : public GMSH_Post_Plugin
+class GMSH_GSHHSPlugin:public GMSH_Post_Plugin
{
- public:
+public:
void getName(char *name) const;
void getInfos(char *author, char *copyright, char *help_text) const;
void catchErrorMessage(char *errorMessage) const;
int getNbOptions() const;
int getNbOptionsStr() const;
- StringXNumber* getOption(int iopt);
- StringXString* getOptionStr(int iopt);
+ StringXNumber *getOption(int iopt);
+ StringXString *getOptionStr(int iopt);
PView *execute(PView *);
};
diff --git a/Plugin/GeoEarthImport.cpp b/Plugin/GeoEarthImport.cpp
index 9fd72ac328..dde759acf0 100644
--- a/Plugin/GeoEarthImport.cpp
+++ b/Plugin/GeoEarthImport.cpp
@@ -1,76 +1,101 @@
#include "GeoEarthImport.h"
-GeoEarthImport::~GeoEarthImport(){
- file.close();
+GeoEarthImport::~GeoEarthImport()
+{
+ file << "Euclidian Coordinates;";
+ file.close();
}
-GeoEarthImport::GeoEarthImport(const std::string _filename){
- filename=_filename;
- file.open(filename.c_str());
- size_field=NULL;
- loop_buff.precision(16);
- std::ostringstream buff;
- il=ip=ill=is=ifi=0;
- buff<<"IP = newp;\n";
- buff<<"IL = newl;\n";
- buff<<"ILL = newll;\n";
- buff<<"IS = news;\n";
- buff<<"IFI = newf;\n";
- buff<<"Point ( IP + "<<ip++<<" ) = {0, 0, 0 , 1};\n";
- buff<<"Point ( IP + "<<ip++<<" ) = {0, 0, 6.371e6 , 1};\n";
- buff<<"PolarSphere ( IS + "<<is++<<" ) = {IP , IP+1};\n";
- file<<buff.str();
- new_surface();
- new_loop();
+GeoEarthImport::GeoEarthImport(const std::string _filename)
+{
+ filename = _filename;
+ file.open(filename.c_str());
+ size_field = NULL;
+ loop_buff.precision(16);
+ std::ostringstream buff;
+ il = ip = ill = is = ifi = 0;
+ buff << "IP = newp;\n";
+ buff << "IL = newl;\n";
+ buff << "ILL = newll;\n";
+ buff << "IS = news;\n";
+ buff << "IFI = newf;\n";
+ buff << "Point ( IP + " << ip++ << " ) = {0, 0, 0 , 1};\n";
+ buff << "Point ( IP + " << ip++ << " ) = {0, 0, 6.371e6 , 1};\n";
+ buff << "PolarSphere ( IS + " << is++ << " ) = {IP , IP+1};\n";
+ file << buff.str();
+ new_surface();
+ new_loop();
}
-void GeoEarthImport::set_size_field(Field *_size_field){
- size_field=_size_field;
+
+void GeoEarthImport::set_size_field(Field * _size_field)
+{
+ size_field = _size_field;
}
-void GeoEarthImport::add_point(const SPoint3 &point){
- SPoint3 midpoint = point;
- SPoint2 stereo(-point.x()/(1+point.z()),-point.y()/(1+point.z()));
- midpoint+=lastpoint;
- midpoint/=2;
- if( ip==first_point_in_loop || !size_field || point.distance(lastpoint)*6361e3 > (*size_field)(midpoint[0],midpoint[1],midpoint[2])){
- loop_buff<<"Point ( IP + "<<ip++<<" ) = {"<<stereo.x()<<", "<<stereo.y()<<", "<<0<<", 1};\n";
- lastpoint=point;
- }
+void GeoEarthImport::add_point(const SPoint3 & point)
+{
+ SPoint3 midpoint = point;
+ SPoint2 stereo(-point.x() / (1 + point.z()), -point.y() / (1 + point.z()));
+ midpoint += lastpoint;
+ midpoint /= 2;
+ if(ip == first_point_in_loop || !size_field
+ || point.distance(lastpoint) * 6361e3 > (*size_field) (midpoint[0],
+ midpoint[1],
+ midpoint[2])) {
+ loop_buff << "Point ( IP + " << ip++ << " ) = {" << stereo.
+ x() << ", " << stereo.y() << ", " << 0 << ", 1};\n";
+ lastpoint = point;
+ }
}
-void GeoEarthImport::add_point_lon_lat(const SPoint2 &ll){
- SPoint3 xyz( cos(ll.y())*cos(ll.x()), cos(ll.y())*sin(ll.x()), sin(ll.y()));
- add_point(xyz);
+void GeoEarthImport::add_point_lon_lat(const SPoint2 & ll)
+{
+ SPoint3 xyz(cos(ll.y()) * cos(ll.x()), cos(ll.y()) * sin(ll.x()),
+ sin(ll.y()));
+ add_point(xyz);
}
-void GeoEarthImport::new_loop(){
- loop_buff.str("");
- first_point_in_loop=ip;
+
+void GeoEarthImport::new_loop()
+{
+ loop_buff.str("");
+ first_point_in_loop = ip;
}
-void GeoEarthImport::end_loop(){
- if(ip-first_point_in_loop>3){
- loop_buff<<"BSpline ( IL + "<<il++<<" ) = { IP + "<<first_point_in_loop<<" : IP + "<<ip-1<<", IP + "<<first_point_in_loop<<" };\n";
- loop_buff<<"Line Loop ( ILL + "<<ill++<<" ) = { IL + "<<il-1<<" };";
- file<<loop_buff.str();
- if(!empty_surface)
- surface_buff<<", ";
- surface_buff<< "ILL + "<<ill-1;
- empty_surface=false;
- }else{
- ip=first_point_in_loop;
- }
- new_loop();
+
+void GeoEarthImport::end_loop()
+{
+ if(ip - first_point_in_loop > 3) {
+ loop_buff << "BSpline ( IL + " << il++ << " ) = { IP + " <<
+ first_point_in_loop << " : IP + " << ip -
+ 1 << ", IP + " << first_point_in_loop << " };\n";
+ loop_buff << "Line Loop ( ILL + " << ill++ << " ) = { IL + " << il -
+ 1 << " };";
+ file << loop_buff.str();
+ if(!empty_surface)
+ surface_buff << ", ";
+ surface_buff << "ILL + " << ill - 1;
+ empty_surface = false;
+ }
+ else {
+ ip = first_point_in_loop;
+ }
+ new_loop();
}
-void GeoEarthImport::new_surface(){
- surface_buff.str("");
- surface_buff<<"Plane Surface( IS + "<<is++<<" ) = { ";
- first_point_in_surface=ip;
- bool empty_surface=true;
+
+void GeoEarthImport::new_surface()
+{
+ surface_buff.str("");
+ surface_buff << "Plane Surface( IS + " << is++ << " ) = { ";
+ first_point_in_surface = ip;
+ bool empty_surface = true;
}
-void GeoEarthImport::end_surface(){
- if(!empty_surface){
- surface_buff<<"};\n";
- surface_buff<<"Field [ IFI + "<<ifi<<"] = Attractor;\n";
- surface_buff<<"Field [ IFI + "<<ifi<<"].NodesList = { IP + "<<first_point_in_surface <<" : IP + "<<ip-1<<" };";
- ifi++;
- file<<surface_buff.str();
- }
- new_surface();
+
+void GeoEarthImport::end_surface()
+{
+ if(!empty_surface) {
+ surface_buff << "};\n";
+ surface_buff << "Field [ IFI + " << ifi << "] = Attractor;\n";
+ surface_buff << "Field [ IFI + " << ifi << "].NodesList = { IP + " <<
+ first_point_in_surface << " : IP + " << ip - 1 << " };";
+ ifi++;
+ file << surface_buff.str();
+ }
+ new_surface();
}
/*static void projectLatLon(Point3D &psphere,Point &pplan){
diff --git a/Plugin/GeoEarthImport.h b/Plugin/GeoEarthImport.h
index bc65a3fd56..c027097fb8 100644
--- a/Plugin/GeoEarthImport.h
+++ b/Plugin/GeoEarthImport.h
@@ -5,25 +5,25 @@
#include <sstream>
#include <fstream>
#include "GeoStringInterface.h"
-class GeoEarthImport{
- Field *size_field;
- std::ostringstream loop_buff,surface_buff;;
- std::string filename;
- std::ofstream file;
- int il,ip,is,ill,ifi;
- int first_point_in_loop,first_point_in_surface;
- bool empty_surface;
- SPoint3 lastpoint;
- void new_surface();
- void new_loop();
+class GeoEarthImport
+{
+ Field *size_field;
+ std::ostringstream loop_buff, surface_buff;;
+ std::string filename;
+ std::ofstream file;
+ int il, ip, is, ill, ifi;
+ int first_point_in_loop, first_point_in_surface;
+ bool empty_surface;
+ SPoint3 lastpoint;
+ void new_surface();
+ void new_loop();
- public :
- GeoEarthImport(const std::string _filename);
- ~GeoEarthImport();
- void set_size_field(Field *_size_field);
- void add_point(const SPoint3 &point);
- void add_point_lon_lat(const SPoint2 &ll);
- void end_loop();
- void end_surface();
+ public: GeoEarthImport(const std::string _filename);
+ ~GeoEarthImport();
+ void set_size_field(Field * _size_field);
+ void add_point(const SPoint3 & point);
+ void add_point_lon_lat(const SPoint2 & ll);
+ void end_loop();
+ void end_surface();
};
#endif
--
GitLab