src/utilities.cpp

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