src/libengrid/surfaceprojection.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2012 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 "surfaceprojection.h"
  24
  25 SurfaceProjection::SurfaceProjection(CadInterface *cad_interface)
  26 {
  27   m_CadInterface = cad_interface;
  28 }
  29
  30 void SurfaceProjection::setForegroundGrid(vtkUnstructuredGrid *grid)
  31 {
  32   m_FGrid = grid;
  33   m_FPart.trackGrid(m_FGrid);
  34 }
  35
  36 vec3_t SurfaceProjection::projectNode(vec3_t x, vtkIdType id_node, bool, vec3_t v, bool strict_direction, bool allow_search)
  37 {
  38   vec3_t n = v;
  39   if (n.abs() < 1e-3) {
  40     if (id_node == -1) {
  41       EG_BUG;
  42     }
  43     n = m_FPart.globalNormal(id_node);
  44   }
  45   if (!checkVector(x)) {
  46     EG_BUG;
  47   }
  48   if (!checkVector(n)) {
  49     cout << "vector defect (id_node=" << id_node << ")" << endl;
  50     return x;
  51     EG_BUG;
  52   }
  53
  54   vec3_t x_proj = x;
  55   m_LastNormal = n;
  56   m_LastRadius = 1e10;
  57
  58   vec3_t x_hit1, n_hit1, x_hit2, n_hit2;
  59   double r_hit1, r_hit2;
  60   CadInterface::HitType hit_type1, hit_type2;
  61
  62   hit_type1 = m_CadInterface->shootRay(x, n, x_hit1, n_hit1, r_hit1);
  63   if (hit_type1 == CadInterface::Miss && !strict_direction) {
  64     n *= -1;
  65     hit_type1 = m_CadInterface->shootRay(x, n, x_hit1, n_hit1, r_hit1);
  66   }
  67   if (hit_type1 == CadInterface::Miss) {
  68     m_Failed = true;
  69     return x;
  70   }
  71   m_Failed = false;
  72   n *= -1;
  73   x_proj = x_hit1;
  74   m_LastNormal = n_hit1;
  75   m_LastRadius = r_hit1;
  76   if (!strict_direction) {
  77     hit_type2 = m_CadInterface->shootRay(x_hit1, n, x_hit2, n_hit2, r_hit2);
  78     if (hit_type2 != CadInterface::Miss) {
  79       if ((x - x_hit2).abs() < (x - x_hit1).abs()) {
  80         x_proj = x_hit2;
  81         m_LastNormal = n_hit2;
  82         m_LastRadius = r_hit2;
  83       }
  84     }
  85   }
  86
  87   if (id_node != -1 && allow_search) {
  88     EG_VTKDCN(vtkDoubleArray, cl, m_FGrid, "node_meshdensity_desired");
  89     double L = 0.5*cl->GetValue(id_node);
  90     if ((x - x_proj).abs() > L) {
  91       vec3_t x_old;
  92       m_FGrid->GetPoint(id_node, x_old.data());
  93       double w = 0.1;
  94       vec3_t x_corr = w*x_proj + (1-w)*x_old;
  95       m_FGrid->GetPoints()->SetPoint(id_node, x_corr.data());
  96       x_corr = projectNode(x_corr, id_node, true, m_FPart.globalNormal(id_node), false, false);
  97       m_FGrid->GetPoints()->SetPoint(id_node, x_old.data());
  98       if ((x_corr - x_proj).abs() > L) {
  99         m_Failed = true;
 100         return x;
 101       }
 102     }
 103   }
 104   /*
 105   {
 106     vec3_t x_old;
 107     vec3_t xp = x_proj;
 108     double scal, w=0;
 109     int count = 0;
 110     do {
 111       if (count >= 10) {
 112         m_Failed = true;
 113         return x;
 114       }
 115       vec3_t n1 = m_FPart.globalNormal(id_node);
 116       m_FGrid->GetPoint(id_node, x_old.data());
 117       m_FGrid->GetPoints()->SetPoint(id_node, x_proj.data());
 118       vec3_t n2 = m_FPart.globalNormal(id_node);
 119       m_FGrid->GetPoints()->SetPoint(id_node, x_old.data());
 120       scal = n1*n2;
 121       if (scal < 0.5) {
 122         w = min(1.0, w + 0.1);
 123         x_proj = w*x + (1-w)*xp;
 124       }
 125       ++count;
 126     } while (scal < 0.5);
 127   }
 128   */
 129
 130   return x_proj;
 131 }
 132
 133 vec3_t SurfaceProjection::snapNode(vec3_t x, vtkIdType, bool)
 134 {
 135   vec3_t x_snap = m_CadInterface->snap(x);
 136   m_LastNormal = m_CadInterface->getLastNormal();
 137   m_LastRadius = m_CadInterface->getLastRadius();
 138   return x_snap;
 139 }
 140
 141
 142 double SurfaceProjection::getRadius(vtkIdType id_node)
 143 {
 144   vec3_t x;
 145   m_FGrid->GetPoint(id_node, x.data());
 146   m_ForceRay = true;
 147   snapNode(x, id_node, false);
 148   m_ForceRay = false;
 149   return m_LastRadius;
 150 }
 151