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