src/libengrid/utilities.cpp

   1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   2 // +                                                                      +
   3 // + This file is part of enGrid.                                         +
   4 // +                                                                      +
   5 // + Copyright 2008-2014 enGits GmbH                                      +
   6 // +                                                                      +
   7 // + enGrid is free software: you can redistribute it and/or modify       +
   8 // + it under the terms of the GNU General Public License as published by +
   9 // + the Free Software Foundation, either version 3 of the License, or    +
  10 // + (at your option) any later version.                                  +
  11 // +                                                                      +
  12 // + enGrid is distributed in the hope that it will be useful,            +
  13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  15 // + GNU General Public License for more details.                         +
  16 // +                                                                      +
  17 // + You should have received a copy of the GNU General Public License    +
  18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  19 // +                                                                      +
  20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  21 #include "utilities.h"
  22
  23 #include <cmath>
  24
  25 #include <QApplication>
  26 #include <QPushButton>
  27 #include <QHBoxLayout>
  28 #include <QFormLayout>
  29 #include <QtDebug>
  30 #include <QMainWindow>
  31 #include <QHBoxLayout>
  32 #include <QString>
  33 #include <QStringList>
  34 #include <QDir>
  35 #include <QtDebug>
  36 #include <QFileInfo>
  37 // #include <QFileDialogArgs>
  38
  39 #include "error.h"
  40 #include "engrid.h"
  41 #include "math/mathvector.h"
  42 #include "math/smallsquarematrix.h"
  43 #include "math/linsolve.h"
  44
  45 #include "vtkIntArray.h"
  46 #include "vtkCellData.h"
  47 #include "vtkDataSet.h"
  48 #include "egvtkobject.h"
  49 #include "vtkCell.h"
  50
  51 double toDouble(QString str) {
  52   str.replace(QString(QObject::tr(",")), QString(QObject::tr(".")));
  53   return(str.toDouble());
  54 }
  55
  56 QString toString(double x, QString separator) {
  57   QString ret;
  58   ret.setNum(x);
  59   ret.replace(QString(QObject::tr(".")), separator);
  60   return(ret);
  61 }
  62
  63 QString addSuffix(QString filename, QString suffix, bool remove_old_suffix) {
  64   QFileInfo fileinfo(filename);
  65   if (remove_old_suffix) {
  66     return fileinfo.completeBaseName() + QObject::tr(".") + suffix;
  67   } else {
  68     if (fileinfo.suffix().toLower() == suffix) {
  69       return fileinfo.absoluteFilePath();
  70     } else {
  71       return fileinfo.absoluteFilePath() + QObject::tr(".") + suffix;
  72     }
  73   }
  74 }
  75
  76 Qt::CheckState int2CheckState(int a) {
  77   if (a == 0) return(Qt::Unchecked);
  78   if (a == 1) return(Qt::PartiallyChecked);
  79   else return(Qt::Checked);
  80 }
  81
  82 int CheckState2int(Qt::CheckState a) {
  83   if (a == Qt::Unchecked) return(0);
  84   if (a == Qt::PartiallyChecked) return(1);
  85   else return(2);
  86 }
  87
  88 QString vector2csv(QVector <double> vector) {
  89   QString csv;
  90   for (int i = 0; i < vector.size(); i++) {
  91     csv += QString::number(vector[i]);
  92     if (i < vector.size() - 1) csv += QObject::tr(";");
  93   }
  94   return csv;
  95 }
  96
  97 QVector <double> csv2vector(QString csv)
  98 {
  99   QVector <double> vector;
 100   if (csv.isEmpty()) return vector;
 101   QStringList list = csv.split(QObject::tr(";"));
 102   foreach(QString str, list) {
 103     vector.push_back(str.toDouble());
 104   }
 105   return vector;
 106 }
 107
 108 mat3_t rotationMatrix_X(double a_rad)
 109 {
 110   mat3_t Rx;
 111   Rx[0][0] = 1; Rx[0][1] = 0;           Rx[0][2] = 0;
 112   Rx[1][0] = 0; Rx[1][1] = cos(a_rad);  Rx[1][2] = sin(a_rad);
 113   Rx[2][0] = 0; Rx[2][1] = -sin(a_rad); Rx[2][2] = cos(a_rad);
 114   return Rx;
 115 }
 116
 117 mat3_t rotationMatrix_Y(double a_rad)
 118 {
 119   mat3_t Ry;
 120   Ry[0][0] = cos(a_rad); Ry[0][1] = 0; Ry[0][2] = -sin(a_rad);
 121   Ry[1][0] = 0;          Ry[1][1] = 1; Ry[1][2] = 0;
 122   Ry[2][0] = sin(a_rad); Ry[2][1] = 0; Ry[2][2] = cos(a_rad);
 123   return Ry;
 124 }
 125
 126 mat3_t rotationMatrix_Z(double a_rad)
 127 {
 128   mat3_t Rz;
 129   Rz[0][0] = cos(a_rad);  Rz[0][1] = sin(a_rad); Rz[0][2] = 0;
 130   Rz[1][0] = -sin(a_rad); Rz[1][1] = cos(a_rad); Rz[1][2] = 0;
 131   Rz[2][0] = 0;           Rz[2][1] = 0;          Rz[2][2] = 1;
 132   return Rz;
 133 }
 134
 135 mat3_t rotateRelativeZXZ(double angle_1_rad, double angle_2_rad, double angle_3_rad)
 136 {
 137   return rotationMatrix_Z(angle_3_rad)*rotationMatrix_X(angle_2_rad)*rotationMatrix_Z(angle_1_rad);
 138 }
 139
 140 mat3_t rotateAbsoluteZXY(double angle_1_rad, double angle_2_rad, double angle_3_rad)
 141 {
 142   return rotationMatrix_Z(angle_1_rad)*rotationMatrix_X(angle_2_rad)*rotationMatrix_Y(angle_3_rad);
 143 }
 144
 145 double getGamma(vec3_t V) {
 146   return V[0] == 0.0 && V[1] == 0.0 && V[2] == 0.0 ? 0.0 : atan2(sqrt(V[0]*V[0] + V[1]*V[1]), V[2]);
 147 }
 148
 149 double getPhi(vec3_t V) {
 150   return V[0] == 0.0 && V[1] == 0.0 ? 0.0 : atan2(V[1], V[0]);
 151 }
 152
 153 QString getDirectory(QWidget * parent, const QString & caption, const QString & selected)
 154 {
 155   qDebug() << "selected=" << selected;
 156   QFileInfo fileinfo(selected);
 157   QString dir = fileinfo.absolutePath();
 158   qDebug() << "dir=" << dir;
 159
 160   // create a qt dialog
 161   QFileDialog dialog(parent, caption, dir);//(args);
 162   dialog.setFileMode(QFileDialog::Directory);
 163 #if QT_VERSION >= 0x040500
 164   dialog.setOption(QFileDialog::ShowDirsOnly, true);
 165 #endif
 166   dialog.selectFile(selected);
 167
 168   if (dialog.exec() == QDialog::Accepted) {
 169     return dialog.selectedFiles().value(0);
 170   }
 171   return QString();
 172 }
 173
 174 int cout_grid(ostream &stream, vtkUnstructuredGrid *grid, bool npoints, bool ncells, bool points, bool cells) {
 175   stream << "=============" << endl;
 176   if (npoints) stream << "grid->GetNumberOfPoints()=" << grid->GetNumberOfPoints() << endl;
 177   if (ncells) stream << "grid->GetNumberOfCells()=" << grid->GetNumberOfCells() << endl;
 178   if (points) {
 179     for (vtkIdType i = 0; i < grid->GetNumberOfPoints(); ++i) {
 180       vec3_t x;
 181       grid->GetPoint(i, x.data());
 182       stream << "Vertex " << i << " = " << x << endl;
 183     }
 184   }
 185   if (cells) {
 186     EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 187     for (vtkIdType i = 0; i < grid->GetNumberOfCells(); ++i) {
 188       vtkCell *C = (vtkCell *) vtkCell::New();
 189       C = grid->GetCell(i);
 190       EG_GET_CELL(i, grid);
 191       stream << "Cell " << i << " = ";
 192       for (int j = 0; j < num_pts; j++)
 193         stream << pts[j] << " ";
 194       stream << "boundary_code=" << cell_code->GetValue(i);
 195       stream << endl;
 196     }
 197   }
 198   stream << "=============" << endl;
 199   return 0;
 200 }
 201
 202 ///////////////////////////////////////////
 203 //Warning: Untested
 204 int addCell(vtkUnstructuredGrid* a_grid, vtkIdType A, vtkIdType B, vtkIdType C, int bc) {
 205   vtkIdType npts = 3;
 206   vtkIdType pts[3];
 207   pts[0] = A;
 208   pts[1] = B;
 209   pts[2] = C;
 210   vtkIdType newCellId = a_grid->InsertNextCell(VTK_TRIANGLE, npts, pts);
 211   EG_VTKDCC(vtkIntArray, cell_code, a_grid, "cell_code");
 212   cell_code->SetValue(newCellId, bc);
 213   return(0);
 214 }
 215
 216 ///////////////////////////////////////////
 217
 218 int getShortestSide(vtkIdType a_id_cell, vtkUnstructuredGrid* a_grid) {
 219   EG_GET_CELL(a_id_cell, a_grid);
 220   vec3_t* x = new vec3_t[num_pts];
 221   for (int i = 0; i < num_pts; i++) a_grid->GetPoints()->GetPoint(pts[i], x[i].data());
 222   int id_minlen = 0;
 223   double minlen = (x[1] - x[0]).abs();
 224   for (int i = 1; i < num_pts; i++) {
 225     double len = (x[(i+1)%num_pts] - x[i]).abs();
 226     if (len < minlen) {
 227       minlen = len;
 228       id_minlen = i;
 229     }
 230   }
 231   delete x;
 232   return(id_minlen);
 233 }
 234
 235 int getLongestSide(vtkIdType a_id_cell, vtkUnstructuredGrid* a_grid) {
 236   EG_GET_CELL(a_id_cell, a_grid);
 237   vec3_t* x = new vec3_t[num_pts];
 238   for (int i = 0; i < num_pts; i++) a_grid->GetPoints()->GetPoint(pts[i], x[i].data());
 239   int id_maxlen = 0;
 240   double maxlen = (x[1] - x[0]).abs();
 241   cout << "maxlen=" << maxlen << endl;
 242   for (int i = 1; i < num_pts; i++) {
 243     double len = (x[(i+1)%num_pts] - x[i]).abs();
 244     cout << "len[" << i << "]=" << len << endl;
 245     if (len > maxlen) {
 246       maxlen = len;
 247       id_maxlen = i;
 248     }
 249   }
 250   delete x;
 251   return(id_maxlen);
 252 }
 253
 254 int getSide(vtkIdType a_id_cell, vtkUnstructuredGrid* a_grid, vtkIdType a_id_node1, vtkIdType a_id_node2) {
 255   EG_GET_CELL(a_id_cell, a_grid);
 256   QVector <vtkIdType> edge(2);
 257
 258   int n = 0;
 259   for (int i = 0; i < num_pts; i++) {
 260     if (pts[i] == a_id_node1) { edge[0] = i; n++;}
 261     if (pts[i] == a_id_node2) { edge[1] = i; n++;}
 262   }
 263   if (n != 2) {
 264     EG_BUG;
 265     return(-1);
 266   }
 267   qSort(edge.begin(), edge.end());
 268   if (edge[0] == 0 && edge[1] == num_pts - 1) return(num_pts - 1);
 269   else return(edge[0]);
 270 }
 271 ///////////////////////////////////////////
 272
 273 QString cell2str(vtkIdType id_cell, vtkUnstructuredGrid* grid) {
 274   QString tmp;
 275   tmp.setNum(id_cell);
 276   QString txt = "id_cell=" + tmp;
 277
 278   EG_GET_CELL(id_cell, grid);
 279
 280   txt += " [";
 281   for (int i_pts = 0; i_pts < num_pts; ++i_pts) {
 282     tmp.setNum(pts[i_pts]);
 283     txt += tmp;
 284     if (i_pts < num_pts - 1) {
 285       txt += ",";
 286     }
 287   }
 288   txt += "]";
 289   return(txt);
 290 }
 291
 292
 293
 294 ///////////////////////////////////////////
 295
 296 pair<vtkIdType, vtkIdType> OrderedPair(vtkIdType a, vtkIdType b) {
 297   vtkIdType x = min(a, b);
 298   vtkIdType y = max(a, b);
 299   return(pair<vtkIdType, vtkIdType>(x, y));
 300 }
 301
 302 const char* VertexType2Str(char T) {
 303   if (T == EG_SIMPLE_VERTEX)         return("EG_SIMPLE_VERTEX");
 304   if (T == EG_FIXED_VERTEX)          return("EG_FIXED_VERTEX");
 305   if (T == EG_FEATURE_EDGE_VERTEX)   return("EG_FEATURE_EDGE_VERTEX");
 306   if (T == EG_FEATURE_CORNER_VERTEX) return("EG_FEATURE_CORNER_VERTEX");
 307   if (T == EG_BOUNDARY_EDGE_VERTEX)  return("EG_BOUNDARY_EDGE_VERTEX");
 308   else return("Unknown vertex type");
 309 }
 310
 311 char Str2VertexType(QString S) {
 312   if (S == "EG_SIMPLE_VERTEX")         return(EG_SIMPLE_VERTEX);
 313   if (S == "EG_FIXED_VERTEX")          return(EG_FIXED_VERTEX);
 314   if (S == "EG_FEATURE_EDGE_VERTEX")   return(EG_FEATURE_EDGE_VERTEX);
 315   if (S == "EG_FEATURE_CORNER_VERTEX") return(EG_FEATURE_CORNER_VERTEX);
 316   if (S == "EG_BOUNDARY_EDGE_VERTEX")  return(EG_BOUNDARY_EDGE_VERTEX);
 317   else return((char) - 1);
 318 }
 319
 320 bool checkVector(vec3_t V)
 321 {
 322   for (int i = 0; i < 3; i++) {
 323     if (isnan(V[i])) {
 324       return false;
 325     }
 326     if (isinf(V[i])) {
 327       return false;
 328     }
 329   }
 330   return true;
 331 }
 332
 333 bool checkVector(vec2_t V)
 334 {
 335   for (int i = 0; i < 2; i++) {
 336     if (isnan(V[i])) {
 337       return false;
 338     }
 339     if (isinf(V[i])) {
 340       return false;
 341     }
 342   }
 343   return true;
 344 }
 345
 346 bool checkReal(double v)
 347 {
 348   if (isnan(v)) {
 349     return false;
 350   }
 351   if (isinf(v)) {
 352     return false;
 353   }
 354   return true;
 355 }
 356
 357 QDebug operator<<(QDebug dbg, const vec3_t &v)
 358 {
 359   dbg.nospace() << "(" << v[0] << ", " << v[1] << ", " << v[2] << ")";
 360   return dbg.space();
 361 }
 362
 363 QDebug operator<<(QDebug dbg, const vec2_t &v)
 364 {
 365   dbg.nospace() << "(" << v[0] << ", " << v[1] << ")";
 366   return dbg.space();
 367 }
 368
 369 QDebug operator<<(QDebug dbg, const dcmplx &c)
 370 {
 371   dbg.nospace() << real(c)<<" + "<<imag(c)<<" *i";
 372   return dbg.space();
 373 }
 374
 375 dcmplx complex_pow(dcmplx base, double power)
 376 {
 377   dcmplx i(0,1);
 378   double alpha = atan2(imag(base),real(base));
 379   return pow(abs(base),power)*exp(power*alpha*i);
 380 }
 381
 382 #define SWAP(a,b) do { double tmp = b ; b = a ; a = tmp ; } while(0)
 383
 384 int poly_solve_cubic(double a, double b, double c, double * x0, double * x1, double * x2)
 385 {
 386   double q = (a * a - 3 * b);
 387   double r = (2 * a * a * a - 9 * a * b + 27 * c);
 388
 389   double Q = q / 9;
 390   double R = r / 54;
 391
 392   double Q3 = Q * Q * Q;
 393   double R2 = R * R;
 394
 395   double CR2 = 729 * r * r;
 396   double CQ3 = 2916 * q * q * q;
 397
 398   if (R == 0 && Q == 0)
 399   {
 400     *x0 = - a / 3 ;
 401     *x1 = - a / 3 ;
 402     *x2 = - a / 3 ;
 403     return 3 ;
 404   }
 405   else if (CR2 == CQ3)
 406   {
 407     // this test is actually R2 == Q3, written in a form suitable
 408     // for exact computation with integers
 409
 410     // Due to finite precision some double roots may be missed, and
 411     // considered to be a pair of complex roots z = x +/- epsilon i
 412     // close to the real axis.
 413
 414     double sqrtQ = sqrt (Q);
 415
 416     if (R > 0)
 417     {
 418       *x0 = -2 * sqrtQ  - a / 3;
 419       *x1 = sqrtQ - a / 3;
 420       *x2 = sqrtQ - a / 3;
 421     } else {
 422       *x0 = - sqrtQ  - a / 3;
 423       *x1 = - sqrtQ - a / 3;
 424       *x2 = 2 * sqrtQ - a / 3;
 425     }
 426     return 3;
 427   } else if (CR2 < CQ3) { // equivalent to R2 < Q3
 428     double sqrtQ = sqrt (Q);
 429     double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
 430     double theta = acos (R / sqrtQ3);
 431     double norm = -2 * sqrtQ;
 432     *x0 = norm * cos (theta / 3) - a / 3;
 433     *x1 = norm * cos ((theta + 2.0 * M_PI) / 3) - a / 3;
 434     *x2 = norm * cos ((theta - 2.0 * M_PI) / 3) - a / 3;
 435
 436     // Sort *x0, *x1, *x2 into increasing order
 437
 438     if (*x0 > *x1) {
 439       SWAP(*x0, *x1);
 440     }
 441     if (*x1 > *x2)
 442     {
 443       SWAP(*x1, *x2);
 444       if (*x0 > *x1) {
 445         SWAP(*x0, *x1);
 446       }
 447     }
 448
 449     return 3;
 450   } else {
 451     double sgnR = (R >= 0 ? 1 : -1);
 452     double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0/3.0);
 453     double B = Q / A ;
 454     *x0 = A + B - a / 3;
 455     return 1;
 456   }
 457 }
 458
 459 // a x^2 + b x + c = 0
 460 int poly_solve_quadratic(double a, double b, double c, double * x0, double * x1)
 461 {
 462   if (a == 0) {
 463     if (b == 0) {
 464       return(0);
 465     } else {
 466       *x0 = -c / b;
 467       return(1);
 468     }
 469   } else {
 470     double delta = pow(b, 2) - 4 * a * c;
 471     if (delta < 0) {
 472       return(0);
 473     } else {
 474       *x0 = (-b + sqrt(delta)) / (2 * a);
 475       *x1 = (-b - sqrt(delta)) / (2 * a);
 476       if (*x1 < *x0) {
 477         double tmp = *x0;
 478         *x0 = *x1;
 479         *x1 = tmp;
 480       }
 481     }
 482   }
 483   return(2);
 484 }