src/libengrid/utilities.cpp

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