src/libengrid/orthogonalityoptimiser.cpp

   1 #include "orthogonalityoptimiser.h"
   2 #include "surfaceprojection.h"
   3 #include "guimainwindow.h"
   4
   5 OrthogonalityOptimiser::OrthogonalityOptimiser()
   6 {
   7   m_StrictPrismChecking = true;
   8 }
   9
  10 double OrthogonalityOptimiser::faceError(const face_t &face, vec3_t x)
  11 {
  12   vec3_t n = face.c1 - face.c2;
  13   n.normalise();
  14   vec3_t xm = x;
  15   foreach (vtkIdType id_node, face.nodes) {
  16     if (id_node != m_IdNode) {
  17       vec3_t xn;
  18       m_Grid->GetPoint(id_node, xn.data());
  19       xm += xn;
  20     }
  21   }
  22   xm *= 1.0/face.nodes.size();
  23   vec3_t v = x - xm;
  24   if (!checkVector(v)) {
  25     EG_BUG;
  26   }
  27   return sqr(v*n);
  28 }
  29
  30 double OrthogonalityOptimiser::func(vec3_t x)
  31 {
  32   if (m_IdNode == 427) {
  33     //cout << "stop" << endl;
  34   }
  35   double err = 0;
  36   foreach (int i_face, m_Node2Face[m_IdNode]) {
  37     double e = faceError(m_Faces[i_face], x);
  38     if (isnan(e)) {
  39       EG_BUG;
  40     }
  41     if (isinf(e)) {
  42       EG_BUG;
  43     }
  44     err += e;
  45   }
  46   return err;
  47 }
  48
  49 void OrthogonalityOptimiser::getFace(vtkIdType id_cell, int subidx, QSet<vtkIdType> &nodes)
  50 {
  51   nodes.clear();
  52   vtkIdType N_pts, *pts;
  53   m_Grid->GetCellPoints(id_cell, N_pts, pts);
  54   vtkIdType type_cell = m_Grid->GetCellType(id_cell);
  55   if (type_cell == VTK_TETRA) {
  56     if (subidx == 0) {
  57       nodes.insert(pts[2]);
  58       nodes.insert(pts[1]);
  59       nodes.insert(pts[0]);
  60     }
  61     if (subidx == 1) {
  62       nodes.insert(pts[1]);
  63       nodes.insert(pts[3]);
  64       nodes.insert(pts[0]);
  65     }
  66     if (subidx == 2) {
  67       nodes.insert(pts[3]);
  68       nodes.insert(pts[2]);
  69       nodes.insert(pts[0]);
  70     }
  71     if (subidx == 3) {
  72       nodes.insert(pts[2]);
  73       nodes.insert(pts[3]);
  74       nodes.insert(pts[1]);
  75     }
  76   }
  77   if (type_cell == VTK_WEDGE) {
  78     if (subidx == 0) {
  79       nodes.insert(pts[0]);
  80       nodes.insert(pts[1]);
  81       nodes.insert(pts[2]);
  82     }
  83     if (subidx == 1) {
  84       nodes.insert(pts[3]);
  85       nodes.insert(pts[5]);
  86       nodes.insert(pts[4]);
  87     }
  88     if (subidx == 2) {
  89       nodes.insert(pts[3]);
  90       nodes.insert(pts[4]);
  91       nodes.insert(pts[1]);
  92       nodes.insert(pts[0]);
  93     }
  94     if (subidx == 3) {
  95       nodes.insert(pts[1]);
  96       nodes.insert(pts[4]);
  97       nodes.insert(pts[2]);
  98       nodes.insert(pts[5]);
  99     }
 100     if (subidx == 4) {
 101       nodes.insert(pts[0]);
 102       nodes.insert(pts[2]);
 103       nodes.insert(pts[5]);
 104       nodes.insert(pts[3]);
 105     }
 106   }
 107   if (type_cell == VTK_HEXAHEDRON) {
 108     if (subidx == 0) {
 109       nodes.insert(pts[0]);
 110       nodes.insert(pts[3]);
 111       nodes.insert(pts[2]);
 112       nodes.insert(pts[1]);
 113     }
 114     if (subidx == 1) {
 115       nodes.insert(pts[4]);
 116       nodes.insert(pts[5]);
 117       nodes.insert(pts[6]);
 118       nodes.insert(pts[7]);
 119     }
 120     if (subidx == 2) {
 121       nodes.insert(pts[0]);
 122       nodes.insert(pts[1]);
 123       nodes.insert(pts[5]);
 124       nodes.insert(pts[4]);
 125     }
 126     if (subidx == 3) {
 127       nodes.insert(pts[3]);
 128       nodes.insert(pts[7]);
 129       nodes.insert(pts[6]);
 130       nodes.insert(pts[2]);
 131     }
 132     if (subidx == 4) {
 133       nodes.insert(pts[0]);
 134       nodes.insert(pts[4]);
 135       nodes.insert(pts[7]);
 136       nodes.insert(pts[3]);
 137     }
 138     if (subidx == 5) {
 139       nodes.insert(pts[1]);
 140       nodes.insert(pts[2]);
 141       nodes.insert(pts[6]);
 142       nodes.insert(pts[5]);
 143     }
 144   }
 145 }
 146
 147 int OrthogonalityOptimiser::getNumberOfFaces(vtkIdType id_cell)
 148 {
 149   vtkIdType cell_type = m_Grid->GetCellType(id_cell);
 150   int N = 0;
 151   if      (cell_type == VTK_TETRA)      N = 4;
 152   else if (cell_type == VTK_PYRAMID)    N = 5;
 153   else if (cell_type == VTK_WEDGE)      N = 5;
 154   else if (cell_type == VTK_HEXAHEDRON) N = 6;
 155   return N;
 156 }
 157
 158 void OrthogonalityOptimiser::buildNode2Face()
 159 {
 160   l2g_t cells = m_Part.getCells();
 161   m_Node2Face.resize(m_Grid->GetNumberOfPoints());
 162   QList<face_t> faces;
 163   foreach (vtkIdType id_cell1, cells) {
 164     if (!isSurface(id_cell1, m_Grid)) {
 165       for (int i = 0; i < m_Part.c2cGSize(id_cell1); ++i) {
 166         vtkIdType id_cell2 = m_Part.c2cGG(id_cell1, i);
 167         if (id_cell1 < id_cell2 && !isSurface(id_cell2, m_Grid)) {
 168           for (int j = 0; j < getNumberOfFaces(id_cell1); ++j) {
 169             QSet<vtkIdType> nodes1;
 170             getFace(id_cell1, j, nodes1);
 171             for (int k = 0; k < getNumberOfFaces(id_cell2); ++k) {
 172               QSet<vtkIdType> nodes2;
 173               getFace(id_cell2, j, nodes2);
 174               if (nodes1 == nodes2) {
 175                 face_t F;
 176                 F.c1 = cellCentre(m_Grid, id_cell1);
 177                 F.c2 = cellCentre(m_Grid, id_cell2);
 178                 F.nodes.resize(nodes1.size());
 179                 qCopy(nodes1.begin(), nodes1.end(), F.nodes.begin());
 180                 faces.append(F);
 181               }
 182             }
 183           }
 184         }
 185       }
 186     }
 187   }
 188   m_Faces.resize(faces.size());
 189   qCopy(faces.begin(), faces.end(), m_Faces.begin());
 190   for (int i = 0; i < m_Faces.size(); ++i) {
 191     foreach (vtkIdType id_node, m_Faces[i].nodes) {
 192       m_Node2Face[id_node].append(i);
 193     }
 194   }
 195 }
 196
 197 void OrthogonalityOptimiser::fixBoundaryNodes()
 198 {
 199   l2g_t cells = m_Part.getCells();
 200   m_Fixed.fill(false, m_Grid->GetNumberOfPoints());
 201   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 202   QVector<QSet<int> > n2bc(m_Grid->GetNumberOfPoints());
 203   foreach (vtkIdType id_cell, cells) {
 204     if (isSurface(id_cell, m_Grid)) {
 205       vtkIdType N_pts, *pts;
 206       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 207       for (int i = 0; i < N_pts; ++i) {
 208         n2bc[pts[i]].insert(cell_code->GetValue(id_cell));
 209       }
 210     }
 211   }
 212   m_Projection.fill(NULL, m_Grid->GetNumberOfPoints());
 213   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 214     if (n2bc[id_node].size() > 1) {
 215       m_Fixed[id_node] = true;
 216     } else if (n2bc[id_node].size() == 1) {
 217       GuiMainWindow::pointer()->getSurfProj(n2bc[id_node].toList().first());
 218     }
 219     if (m_Node2Face[id_node].size() == 0) {
 220       m_Fixed[id_node] = true;
 221     }
 222   }
 223 }
 224
 225 vec3_t OrthogonalityOptimiser::newPosition(vtkIdType id_node)
 226 {
 227   vec3_t x_old;
 228   m_Grid->GetPoint(id_node, x_old.data());
 229   if (m_Node2Face[id_node].size() == 0) {
 230     return x_old;
 231   }
 232   vec3_t x_new(0,0,0);
 233   foreach (int i_face, m_Node2Face[id_node]) {
 234     const face_t &face = m_Faces[i_face];
 235     vec3_t n = face.c1 - face.c2;
 236     n.normalise();
 237     vec3_t xm(0,0,0);
 238     foreach (vtkIdType id_node, face.nodes) {
 239       vec3_t xn;
 240       m_Grid->GetPoint(id_node, xn.data());
 241       xm += xn;
 242     }
 243     xm *= 1.0/face.nodes.size();
 244     vec3_t v = x_old - xm;
 245     if (!checkVector(v)) {
 246       EG_BUG;
 247     }
 248     vec3_t Dx = (v*n)*n;
 249     x_new += x_old - Dx;
 250   }
 251   x_new *= (1.0/m_Node2Face[id_node].size());
 252   return x_new;
 253 }
 254
 255 void OrthogonalityOptimiser::operate()
 256 {
 257   computeNormals();
 258   buildNode2Face();
 259   fixBoundaryNodes();
 260   l2g_t nodes = m_Part.getNodes();
 261   int N = 0;
 262   foreach (vtkIdType id_node, nodes) {
 263     if (!m_Fixed[id_node]) {
 264       m_IdNode = id_node;
 265       vec3_t x_old;
 266       m_Grid->GetPoint(id_node, x_old.data());
 267       //vec3_t x_new = optimise(x_old);
 268       vec3_t x_new = newPosition(id_node);
 269       if (checkVector(x_new)) {
 270         if (m_Projection[id_node]) {
 271           x_new = m_Projection[id_node]->projectRestricted(x_new, id_node);
 272         }
 273         if (checkVector(x_new)) {
 274           vec3_t Dx = x_new - x_old;
 275           if (moveNode(id_node, Dx)) ++N;
 276         }
 277       }
 278     }
 279   }
 280   cout << N << " nodes optimised" << endl;
 281 }