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