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