createspecialmapping.h

   1 #ifndef CREATESPECIALMAPPING_H
   2 #define CREATESPECIALMAPPING_H
   3
   4 #include <vtkUnstructuredGrid.h>
   5 #include <vtkPolyData.h>
   6 #include <QSet>
   7 #include <QVector>
   8 #include "egvtkobject.h"
   9 #include "operation.h"
  10 #include "vertexmeshdensity.h"
  11
  12 #include "geometrytools.h"
  13 using namespace GeometryTools;
  14
  15 #include <cmath>
  16 using namespace std;
  17
  18 class CreateSpecialMapping : public Operation {
  19   public:
  20     CreateSpecialMapping();
  21     int Process();
  22     void operate(){};
  23
  24     int N_SmoothIterations;
  25     int maxiter_density;
  26
  27     double Convergence_meshdensity;
  28
  29     bool insert_FP;
  30     bool insert_EP;
  31     bool remove_FP;
  32     bool remove_EP;
  33
  34     bool DoSwap;
  35     bool DoLaplaceSmoothing;
  36
  37     int N_inserted_FP;
  38     int N_inserted_EP;
  39     int N_removed_FP;
  40     int N_removed_EP;
  41
  42     int N_points;
  43     int N_cells;
  44     int N_newpoints;
  45     int N_newcells;
  46     int m_total_N_newpoints;
  47     int m_total_N_newcells;
  48
  49     QSet<int> m_bcs;
  50     QVector<vtkIdType> m_AllCells;
  51     QVector<vtkIdType> m_SelectedCells;
  52     vtkUnstructuredGrid* m_grid;
  53     vtkIdType m_newNodeId;
  54
  55     QMap< pair<vtkIdType,vtkIdType>, vtkIdType> edge_map;
  56     QVector <stencil_t> StencilVector;
  57     QSet <vtkIdType> m_SelectedNodes;
  58
  59     QVector <int> hitlist;//Elements to be terminated (0=keep alive, 1=field agent to eliminate, 2=border agent to eliminate)
  60     QVector <int> offset;//offset caused by terminated elements
  61
  62 /*    double SV_value;
  63     double FV_value;
  64     double FEV_value;
  65     double BEV_value;*/
  66
  67     QVector <VertexMeshDensity> VMDvector;//Vertices of Mass destruction
  68
  69     QMap <vtkIdType,bool> marked_cells;
  70     QMap <vtkIdType,bool> marked_nodes;
  71
  72     void SetInput(QSet<int> a_bcs,vtkUnstructuredGrid* a_grid)
  73     {
  74       m_bcs=a_bcs;
  75       m_grid=a_grid;
  76     };
  77
  78     void SetVertexMeshDensityVector(QVector <VertexMeshDensity> a_VMDvector){VMDvector=a_VMDvector;};
  79
  80 /*    void Set_SV_value(double V){SV_value=V;};
  81     void Set_FV_value(double V){FV_value=V;};
  82     void Set_FEV_value(double V){FEV_value=V;};
  83     void Set_BEV_value(double V){BEV_value=V;};*/
  84
  85     void SetConvergence_meshdensity(double C){Convergence_meshdensity=C;};
  86     void Set_insert_FP(bool B){insert_FP=B;};
  87     void Set_insert_EP(bool B){insert_EP=B;};
  88     void Set_remove_FP(bool B){remove_FP=B;};
  89     void Set_remove_EP(bool B){remove_EP=B;};
  90
  91     VertexMeshDensity getVMD(vtkIdType node, char VertexType);
  92
  93   //utilities
  94   public:
  95     double Um(vtkIdType D) {
  96       double ret=0;
  97       vtkIdType N_pts, *pts;
  98       m_grid->GetCellPoints(D, N_pts, pts);
  99       for(int i=0;i<N_pts;i++)
 100       {
 101         vec3_t A,B;
 102         m_grid->GetPoints()->GetPoint(pts[i], A.data());
 103         m_grid->GetPoints()->GetPoint(pts[(i+1)%N_pts], B.data());
 104         ret+=(B-A).abs();
 105       }
 106       return(ret);
 107     }
 108     double A_U(vtkIdType D) { // area of the circumscribed circle of the triangle
 109       vtkIdType N_pts, *pts;
 110       m_grid->GetCellPoints(D, N_pts, pts);
 111       vec3_t A,B,C;
 112       m_grid->GetPoints()->GetPoint(pts[0], A.data());
 113       m_grid->GetPoints()->GetPoint(pts[1], B.data());
 114       m_grid->GetPoints()->GetPoint(pts[2], C.data());
 115       double a=(C-B).abs();
 116       double alpha=angle((B-A),(C-A));
 117       double R=a/(2*sin(alpha));
 118       return(M_PI*R*R);
 119     }
 120     double A_D(vtkIdType D) { // triangle area
 121       return(cellVA(m_grid,D));
 122     }
 123     double DN(int i,vtkIdType D) { //triangle neighbours
 124       return(c2c[D][i]);
 125     }
 126     double nk(vtkIdType P) {
 127       return(n2n[P].size());
 128     }
 129     double G_k(vtkIdType node) {
 130       EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 131       return(1.0/node_meshdensity->GetValue(node));
 132     };
 133     double DK(int i,vtkIdType D) { // triangle nodes
 134       vtkIdType N_pts, *pts;
 135       m_grid->GetCellPoints(D, N_pts, pts);
 136       return(pts[i]);
 137     }
 138     vtkIdType KK(int i,vtkIdType j,vtkIdType K) {//i=1 or 2, j=node2, K=node1
 139       if(i==1) return(K);
 140       else return(j);
 141     }
 142     double L_k(vtkIdType j,vtkIdType K)// node1 K, node2 j
 143     {
 144         vec3_t A;
 145         vec3_t B;
 146         m_grid->GetPoints()->GetPoint(K, A.data());
 147         m_grid->GetPoints()->GetPoint(j, B.data());
 148         return((B-A).abs());
 149     }
 150     double Q_L(vtkIdType D)
 151     {
 152       // Um(D)/sum(G_k(DK(i,D)),i,1,3)
 153       double denom_sum=0;
 154       for(int i=0;i<3;i++)
 155       {
 156         denom_sum += G_k(DK(i,D));
 157       }
 158       /*if(DebugLevel>0) cout<<"D="<<D<<" Um(D)="<<Um(D)<<" denom_sum="<<denom_sum<<endl;*/
 159       return(Um(D)/denom_sum);
 160     }
 161     double Q_L1(vtkIdType P)
 162     {
 163       // [2*sum(L_k(i~),i,1,nk(P))]/[sum(G_k(KK(1,i~))+G_k(KK(2,i~)),i,1,nk(P))]
 164       double num_sum=0;
 165       double denom_sum=0;
 166       foreach(vtkIdType j,n2n[P])
 167       {
 168         num_sum += 2*L_k(j,P);
 169         denom_sum += G_k(KK(1,j,P))+G_k(KK(2,j,P));
 170       }
 171       return(num_sum/denom_sum);
 172     }
 173     double Q_L2(vtkIdType P)
 174     {
 175
 176       // min([2*L_k(i~)]/[G_k(KK(1,i~))+G_k(KK(2,i~))])
 177       QVector <double> V;
 178       double num,denom;
 179       foreach(vtkIdType j,n2n[P])
 180       {
 181         num = 2*L_k(j,P);
 182         denom = G_k(KK(1,j,P))+G_k(KK(2,j,P));
 183         V.push_back(num/denom);
 184       }
 185       qSort(V.begin(),V.end());
 186       return(V[0]);
 187     }
 188     double T_min(int w)
 189     {
 190       // sum([A_U(i)]/[A_D(i)^w]*[G_k(i)^(2*(w-1))],i,1,Nd)
 191       int N_cells=m_grid->GetNumberOfCells();
 192       double T=0;
 193       for(int i=0;i<N_cells;i++)
 194       {
 195         T += A_U(i)/pow(A_D(i),w)*pow(G_k(i),2*(w-1));
 196       }
 197       return(T);
 198     }
 199     bool insert_fieldpoint(vtkIdType D)
 200     {
 201       double Fred1=1.0/sqrt(3);
 202       double Qmin=1.1;//1.189;
 203       double total=0;
 204       for(int i=0;i<3;i++)
 205       {
 206         vtkIdType cell=DN(i,D);
 207         if(cell!=-1) total += Q_L(cell);
 208       }
 209 /*      cout<<"Q_L(D)>1.0/Fred1="<<Q_L(D)<<">"<<1.0/Fred1<<endl;
 210       cout<<"total>3*Qmin="<<total<<">"<<3*Qmin<<endl;*/
 211       return ( Q_L(D)>1.0/Fred1 && total>3*Qmin );
 212     }
 213     bool insert_edgepoint(vtkIdType j,vtkIdType K)// node1 K, node2 j
 214     {
 215 /*      cout<<"j="<<j<<endl;
 216       cout<<"K="<<K<<endl;
 217       cout<<"0.5*G_k(K)="<<0.5*G_k(K)<<endl;
 218       cout<<"L_k(j,K)="<<L_k(j,K)<<endl;
 219       cout<<"1*G_k(K)="<<1*G_k(K)<<endl;
 220       cout<<"return ( 0.5*G_k(K)<L_k(j,K) && L_k(j,K)<1*G_k(K) );"<<endl;
 221       return ( 0.5*G_k(K)<L_k(j,K) && L_k(j,K)<1*G_k(K) );*/
 222
 223       bool result=L_k(j,K)>0.5*(G_k(j)+G_k(K));
 224       if(DebugLevel>0 && result){
 225         cout<<"j="<<j<<endl;
 226         cout<<"K="<<K<<endl;
 227         cout<<"G_k(j)="<<G_k(j)<<endl;
 228         cout<<"G_k(K)="<<G_k(K)<<endl;
 229         cout<<"0.5*(G_k(j)+G_k(K))="<<0.5*(G_k(j)+G_k(K))<<endl;
 230         cout<<"L_k(j,K)="<<L_k(j,K)<<endl;
 231       }
 232       return ( result );
 233     }
 234     bool remove_fieldpoint(vtkIdType P)
 235     {
 236       double QL1max=0.8;
 237       double QL2max=0.5;
 238       bool result = Q_L1(P)<QL1max && Q_L2(P)<QL2max;
 239       if(DebugLevel>0 && result)
 240       {
 241         cout<<"Q_L1(P)<QL1max="<< Q_L1(P)<< "<" << QL1max<<endl;
 242         cout<<"Q_L2(P)<QL2max="<< Q_L2(P)<< "<" << QL2max<<endl;
 243       }
 244       return ( result );
 245     }
 246     bool remove_edgepoint(vtkIdType P)
 247     {
 248       return ( 0.5*G_k(P)<CurrentVertexAvgDist(P,n2n,m_grid) && CurrentVertexAvgDist(P,n2n,m_grid)<1*G_k(P) );
 249     }
 250     int insert_FP_counter();
 251     int insert_EP_counter();
 252     int remove_FP_counter();
 253     int remove_EP_counter();
 254
 255     int insert_FP_actor(vtkUnstructuredGrid* grid_tmp);
 256     int insert_EP_actor(vtkUnstructuredGrid* grid_tmp);
 257     int remove_FP_actor(vtkUnstructuredGrid* grid_tmp);
 258     int remove_EP_actor(vtkUnstructuredGrid* grid_tmp);
 259
 260     int insert_FP_all();
 261     int insert_EP_all();
 262     int remove_FP_all();
 263     int remove_EP_all();
 264
 265     int FullEdit();
 266     int UpdateDesiredMeshDensity();
 267 //     int UpdateNodeType();
 268 //     bool DeletePoint_2(vtkUnstructuredGrid *src, vtkIdType DeadNode);
 269 //     vtkIdType FindSnapPoint(vtkUnstructuredGrid *src, vtkIdType DeadNode,QSet <vtkIdType> & DeadCells,QSet <vtkIdType> & MutatedCells,QSet <vtkIdType> & MutilatedCells);
 270
 271     int remove_EP_all_2();
 272     int remove_FP_all_2();
 273
 274     int SwapFunction();
 275     int SmoothFunction();
 276
 277 };
 278 //end of CreateSpecialMapping class
 279
 280 // #define VTK_SIMPLE_VERTEX 0
 281 // #define VTK_FIXED_VERTEX 1
 282 // #define VTK_FEATURE_EDGE_VERTEX 2
 283 // #define VTK_BOUNDARY_EDGE_VERTEX 3
 284
 285 #endif