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