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       vtkIdType npts = C->GetNumberOfPoints();
 191       vtkIdType* pts;
 192       grid->GetCellPoints(i, npts, pts);
 193       stream << "Cell " << i << " = ";
 194       for (int j = 0; j < npts; j++) stream << pts[j] << " ";
 195       stream << "boundary_code=" << cell_code->GetValue(i);
 196       stream << endl;
 197     }
 198   }
 199   stream << "=============" << endl;
 200   return 0;
 201 }
 202
 203 ///////////////////////////////////////////
 204 //Warning: Untested
 205 int addCell(vtkUnstructuredGrid* a_grid, vtkIdType A, vtkIdType B, vtkIdType C, int bc) {
 206   vtkIdType npts = 3;
 207   vtkIdType pts[3];
 208   pts[0] = A;
 209   pts[1] = B;
 210   pts[2] = C;
 211   vtkIdType newCellId = a_grid->InsertNextCell(VTK_TRIANGLE, npts, pts);
 212   EG_VTKDCC(vtkIntArray, cell_code, a_grid, "cell_code");
 213   cell_code->SetValue(newCellId, bc);
 214   return(0);
 215 }
 216
 217 ///////////////////////////////////////////
 218
 219 int getShortestSide(vtkIdType a_id_cell, vtkUnstructuredGrid* a_grid) {
 220   vtkIdType N_pts, *pts;
 221   a_grid->GetCellPoints(a_id_cell, N_pts, pts);
 222   vec3_t* x = new vec3_t[N_pts];
 223   for (int i = 0; i < N_pts; i++) a_grid->GetPoints()->GetPoint(pts[i], x[i].data());
 224   int id_minlen = 0;
 225   double minlen = (x[1] - x[0]).abs();
 226   for (int i = 1; i < N_pts; i++) {
 227     double len = (x[(i+1)%N_pts] - x[i]).abs();
 228     if (len < minlen) {
 229       minlen = len;
 230       id_minlen = i;
 231     }
 232   }
 233   delete x;
 234   return(id_minlen);
 235 }
 236
 237 int getLongestSide(vtkIdType a_id_cell, vtkUnstructuredGrid* a_grid) {
 238   vtkIdType N_pts, *pts;
 239   a_grid->GetCellPoints(a_id_cell, N_pts, pts);
 240   vec3_t* x = new vec3_t[N_pts];
 241   for (int i = 0; i < N_pts; i++) a_grid->GetPoints()->GetPoint(pts[i], x[i].data());
 242   int id_maxlen = 0;
 243   double maxlen = (x[1] - x[0]).abs();
 244   cout << "maxlen=" << maxlen << endl;
 245   for (int i = 1; i < N_pts; i++) {
 246     double len = (x[(i+1)%N_pts] - x[i]).abs();
 247     cout << "len[" << i << "]=" << len << endl;
 248     if (len > maxlen) {
 249       maxlen = len;
 250       id_maxlen = i;
 251     }
 252   }
 253   delete x;
 254   return(id_maxlen);
 255 }
 256
 257 int getSide(vtkIdType a_id_cell, vtkUnstructuredGrid* a_grid, vtkIdType a_id_node1, vtkIdType a_id_node2) {
 258   vtkIdType N_pts, *pts;
 259   a_grid->GetCellPoints(a_id_cell, N_pts, pts);
 260   QVector <vtkIdType> edge(2);
 261
 262   int n = 0;
 263   for (int i = 0; i < N_pts; i++) {
 264     if (pts[i] == a_id_node1) { edge[0] = i; n++;}
 265     if (pts[i] == a_id_node2) { edge[1] = i; n++;}
 266   }
 267   if (n != 2) {
 268     EG_BUG;
 269     return(-1);
 270   }
 271   qSort(edge.begin(), edge.end());
 272   if (edge[0] == 0 && edge[1] == N_pts - 1) return(N_pts - 1);
 273   else return(edge[0]);
 274 }
 275 ///////////////////////////////////////////
 276
 277 QString cell2str(vtkIdType id_cell, vtkUnstructuredGrid* grid) {
 278   QString tmp;
 279   tmp.setNum(id_cell);
 280   QString txt = "id_cell=" + tmp;
 281
 282   vtkIdType N_pts, *pts;
 283   grid->GetCellPoints(id_cell, N_pts, pts);
 284
 285   txt += " [";
 286   for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 287     tmp.setNum(pts[i_pts]);
 288     txt += tmp;
 289     if (i_pts < N_pts - 1) {
 290       txt += ",";
 291     }
 292   }
 293   txt += "]";
 294   return(txt);
 295 }
 296
 297
 298
 299 ///////////////////////////////////////////
 300
 301 pair<vtkIdType, vtkIdType> OrderedPair(vtkIdType a, vtkIdType b) {
 302   vtkIdType x = min(a, b);
 303   vtkIdType y = max(a, b);
 304   return(pair<vtkIdType, vtkIdType>(x, y));
 305 }
 306
 307 const char* VertexType2Str(char T) {
 308   if (T == EG_SIMPLE_VERTEX)         return("EG_SIMPLE_VERTEX");
 309   if (T == EG_FIXED_VERTEX)          return("EG_FIXED_VERTEX");
 310   if (T == EG_FEATURE_EDGE_VERTEX)   return("EG_FEATURE_EDGE_VERTEX");
 311   if (T == EG_FEATURE_CORNER_VERTEX) return("EG_FEATURE_CORNER_VERTEX");
 312   if (T == EG_BOUNDARY_EDGE_VERTEX)  return("EG_BOUNDARY_EDGE_VERTEX");
 313   else return("Unknown vertex type");
 314 }
 315
 316 char Str2VertexType(QString S) {
 317   if (S == "EG_SIMPLE_VERTEX")         return(EG_SIMPLE_VERTEX);
 318   if (S == "EG_FIXED_VERTEX")          return(EG_FIXED_VERTEX);
 319   if (S == "EG_FEATURE_EDGE_VERTEX")   return(EG_FEATURE_EDGE_VERTEX);
 320   if (S == "EG_FEATURE_CORNER_VERTEX") return(EG_FEATURE_CORNER_VERTEX);
 321   if (S == "EG_BOUNDARY_EDGE_VERTEX")  return(EG_BOUNDARY_EDGE_VERTEX);
 322   else return((char) - 1);
 323 }
 324
 325 bool checkVector(vec3_t V)
 326 {
 327   for (int i = 0; i < 3; i++) {
 328     if (isnan(V[i])) {
 329       return false;
 330     }
 331     if (isinf(V[i])) {
 332       return false;
 333     }
 334   }
 335   return true;
 336 }
 337
 338 bool checkVector(vec2_t V)
 339 {
 340   for (int i = 0; i < 2; i++) {
 341     if (isnan(V[i])) {
 342       return false;
 343     }
 344     if (isinf(V[i])) {
 345       return false;
 346     }
 347   }
 348   return true;
 349 }
 350
 351 bool checkReal(double v)
 352 {
 353   if (isnan(v)) {
 354     return false;
 355   }
 356   if (isinf(v)) {
 357     return false;
 358   }
 359   return true;
 360 }
 361
 362 QDebug operator<<(QDebug dbg, const vec3_t &v)
 363 {
 364   dbg.nospace() << "(" << v[0] << ", " << v[1] << ", " << v[2] << ")";
 365   return dbg.space();
 366 }
 367
 368 QDebug operator<<(QDebug dbg, const vec2_t &v)
 369 {
 370   dbg.nospace() << "(" << v[0] << ", " << v[1] << ")";
 371   return dbg.space();
 372 }
 373
 374 QDebug operator<<(QDebug dbg, const dcmplx &c)
 375 {
 376   dbg.nospace() << real(c)<<" + "<<imag(c)<<" *i";
 377   return dbg.space();
 378 }
 379
 380 dcmplx complex_pow(dcmplx base, double power)
 381 {
 382   dcmplx i(0,1);
 383   double alpha = atan2(imag(base),real(base));
 384   return pow(abs(base),power)*exp(power*alpha*i);
 385 }
 386
 387 #define SWAP(a,b) do { double tmp = b ; b = a ; a = tmp ; } while(0)
 388
 389 int poly_solve_cubic(double a, double b, double c, double * x0, double * x1, double * x2)
 390 {
 391   double q = (a * a - 3 * b);
 392   double r = (2 * a * a * a - 9 * a * b + 27 * c);
 393
 394   double Q = q / 9;
 395   double R = r / 54;
 396
 397   double Q3 = Q * Q * Q;
 398   double R2 = R * R;
 399
 400   double CR2 = 729 * r * r;
 401   double CQ3 = 2916 * q * q * q;
 402
 403   if (R == 0 && Q == 0)
 404   {
 405     *x0 = - a / 3 ;
 406     *x1 = - a / 3 ;
 407     *x2 = - a / 3 ;
 408     return 3 ;
 409   }
 410   else if (CR2 == CQ3)
 411   {
 412     // this test is actually R2 == Q3, written in a form suitable
 413     // for exact computation with integers
 414
 415     // Due to finite precision some double roots may be missed, and
 416     // considered to be a pair of complex roots z = x +/- epsilon i
 417     // close to the real axis.
 418
 419     double sqrtQ = sqrt (Q);
 420
 421     if (R > 0)
 422     {
 423       *x0 = -2 * sqrtQ  - a / 3;
 424       *x1 = sqrtQ - a / 3;
 425       *x2 = sqrtQ - a / 3;
 426     } else {
 427       *x0 = - sqrtQ  - a / 3;
 428       *x1 = - sqrtQ - a / 3;
 429       *x2 = 2 * sqrtQ - a / 3;
 430     }
 431     return 3;
 432   } else if (CR2 < CQ3) { // equivalent to R2 < Q3
 433     double sqrtQ = sqrt (Q);
 434     double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
 435     double theta = acos (R / sqrtQ3);
 436     double norm = -2 * sqrtQ;
 437     *x0 = norm * cos (theta / 3) - a / 3;
 438     *x1 = norm * cos ((theta + 2.0 * M_PI) / 3) - a / 3;
 439     *x2 = norm * cos ((theta - 2.0 * M_PI) / 3) - a / 3;
 440
 441     // Sort *x0, *x1, *x2 into increasing order
 442
 443     if (*x0 > *x1) {
 444       SWAP(*x0, *x1);
 445     }
 446     if (*x1 > *x2)
 447     {
 448       SWAP(*x1, *x2);
 449       if (*x0 > *x1) {
 450         SWAP(*x0, *x1);
 451       }
 452     }
 453
 454     return 3;
 455   } else {
 456     double sgnR = (R >= 0 ? 1 : -1);
 457     double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0/3.0);
 458     double B = Q / A ;
 459     *x0 = A + B - a / 3;
 460     return 1;
 461   }
 462 }
 463
 464 // a x^2 + b x + c = 0
 465 int poly_solve_quadratic(double a, double b, double c, double * x0, double * x1)
 466 {
 467   if (a == 0) {
 468     if (b == 0) {
 469       return(0);
 470     } else {
 471       *x0 = -c / b;
 472       return(1);
 473     }
 474   } else {
 475     double delta = pow(b, 2) - 4 * a * c;
 476     if (delta < 0) {
 477       return(0);
 478     } else {
 479       *x0 = (-b + sqrt(delta)) / (2 * a);
 480       *x1 = (-b - sqrt(delta)) / (2 * a);
 481       if (*x1 < *x0) {
 482         double tmp = *x0;
 483         *x0 = *x1;
 484         *x1 = tmp;
 485       }
 486     }
 487   }
 488   return(2);
 489 }