src/mesh/zone.cc

   1 /*****************************************************************************/
   2 /*              8888888         88888888         88888888                    */
   3 /*            8                8                8                            */
   4 /*           8                 8                8                            */
   5 /*           8                  88888888         88888888                    */
   6 /*           8      8888                8                8                   */
   7 /*            8       8                 8                8                   */
   8 /*              888888         888888888        888888888                    */
   9 /*                                                                           */
  10 /*       A Two-Dimensional General Purpose Semiconductor Simulator.          */
  11 /*                                                                           */
  12 /*  GSS 0.4x                                                                 */
  13 /*  Last update: May 29, 2006                                                */
  14 /*                                                                           */
  15 /*  Gong Ding                                                                */
  16 /*  gdiso@ustc.edu                                                           */
  17 /*  NINT, No.69 P.O.Box, Xi'an City, China                                   */
  18 /*                                                                           */
  19 /*****************************************************************************/
  20
  21 #include <list>
  22 #include <stack>
  23 #include <queue>
  24 #include <algorithm>
  25 #include <math.h>
  26 #include "zone.h"
  27 #include "geom.h"
  28 #include "log.h"
  29 using namespace std;
  30
  31 /* ----------------------------------------------------------------------------
  32  * ZONE::build_tricell_data:  This function compute circle center, and the distance between
  33  * the circle center and each edge. if circle center locates out of the triangle, negtive the
  34  * corresponding distance.
  35  */
  36 int ZONE::build_tricell_data()
  37 {
  38   double PI=3.14159265358979323846;
  39   for(int i=0;i<datri.size();i++)
  40   {
  41     double x1,y1,x2,y2,x3,y3;
  42     double a,b,c;
  43     x1 = danode[datri[i].node[0]].x;
  44     y1 = danode[datri[i].node[0]].y;
  45     x2 = danode[datri[i].node[1]].x;
  46     y2 = danode[datri[i].node[1]].y;
  47     x3 = danode[datri[i].node[2]].x;
  48     y3 = danode[datri[i].node[2]].y;
  49     datri[i].edge_len[0] = a = Distance(x3,y3,x2,y2);
  50     datri[i].edge_len[1] = b = Distance(x1,y1,x3,y3);
  51     datri[i].edge_len[2] = c = Distance(x1,y1,x2,y2);
  52     datri[i].angle[0] = acos((b*b+c*c-a*a)/(2*b*c));
  53     datri[i].angle[1] = acos((c*c+a*a-b*b)/(2*c*a));
  54     datri[i].angle[2] = acos((a*a+b*b-c*c)/(2*a*b));
  55     datri[i].area = TriArea(x1,y1,x2,y2,x3,y3);
  56     if(CircleCenter(x1,y1,x2,y2,x3,y3,datri[i].xc,datri[i].yc))
  57     {
  58       sprintf(log_buf,"warning! CircleCenter: the three nodes of Tri %d of zone %d coplanar.\n",i,zone_index);
  59       GSS_LOG();
  60     }
  61     datri[i].d[0] = DistancePointLine(x2,y2,x3,y3,datri[i].xc,datri[i].yc);
  62     datri[i].d[1] = DistancePointLine(x1,y1,x3,y3,datri[i].xc,datri[i].yc);
  63     datri[i].d[2] = DistancePointLine(x1,y1,x2,y2,datri[i].xc,datri[i].yc);
  64     if(CompAngle(x3,y3,x1,y1,x2,y2)>PI/2) datri[i].d[0]*=-1;
  65     if(CompAngle(x1,y1,x2,y2,x3,y3)>PI/2) datri[i].d[1]*=-1;
  66     if(CompAngle(x2,y2,x3,y3,x1,y1)>PI/2) datri[i].d[2]*=-1;
  67     datri[i].s[0] = 0.5*datri[i].edge_len[0]*datri[i].d[0];
  68     datri[i].s[1] = 0.5*datri[i].edge_len[1]*datri[i].d[1];
  69     datri[i].s[2] = 0.5*datri[i].edge_len[2]*datri[i].d[2];
  70     for(int j=0;j<dasegment.size();j++)
  71     {
  72       if(dasegment[j].get_edge_bc_index(datri[i].node[0],datri[i].node[1]))
  73          datri[i].bc[2] = dasegment[j].bc_index;
  74       if(dasegment[j].get_edge_bc_index(datri[i].node[1],datri[i].node[2]))
  75          datri[i].bc[0] = dasegment[j].bc_index;
  76       if(dasegment[j].get_edge_bc_index(datri[i].node[2],datri[i].node[0]))
  77          datri[i].bc[1] = dasegment[j].bc_index;
  78     }
  79   }
  80   return 0;
  81 }
  82
  83 /* ----------------------------------------------------------------------------
  84  * ZONE::setup_voronoicell:  This function setup voronoicell data structure
  85  * the adjacency infomation will be build here.
  86  * if the Triangle is count clock wise ordered, this function works well
  87  */
  88
  89 typedef struct
  90 {
  91   int node1,node2;
  92   int flag1,flag2;
  93 }
  94 ATABLE;
  95
  96
  97 int ZONE::setup_voronoicell()
  98 {
  99   int adjtabsize = 0;
 100
 101   vector<ATABLE> **zone_atable;
 102   zone_atable = new vector<ATABLE>*[danode.size()];
 103   for(int i=0;i<danode.size();i++)
 104     zone_atable[i] = new vector<ATABLE>;
 105
 106   ATABLE tmp_atable;
 107   tmp_atable.flag1=0; tmp_atable.flag2=0;
 108
 109   davcell.clear();
 110   davcell.Init(danode.size());   //init davcell, cell num = node num
 111
 112   //insert triangle into each node's adjacent table
 113   for(int i=0;i<datri.size();i++)
 114   {
 115
 116     //triangle
 117     tmp_atable.node1 =  datri[i].node[1];
 118     tmp_atable.node2 =  datri[i].node[2];
 119     zone_atable[datri[i].node[0]]->push_back(tmp_atable);
 120
 121     tmp_atable.node1 =  datri[i].node[2];
 122     tmp_atable.node2 =  datri[i].node[0];
 123     zone_atable[datri[i].node[1]]->push_back(tmp_atable);
 124
 125     tmp_atable.node1 =  datri[i].node[0];
 126     tmp_atable.node2 =  datri[i].node[1];
 127     zone_atable[datri[i].node[2]]->push_back(tmp_atable);
 128
 129   }
 130
 131
 132   stack<int> clockw;
 133   queue<int> cclockw;
 134   for(int i=0;i<davcell.size();i++)
 135   {
 136     vector<ATABLE>::iterator pt =  zone_atable[i]->begin();
 137     int v1,v2;
 138
 139     v1 = pt->node1;
 140     v2 = pt->node2;
 141     //search it's neighbor triangles. begin from zone_atable[0]
 142     //clock wise search.push the nodes into a stack
 143     clockw.push(v1); pt->flag1=1;
 144     for(int num=0;num<zone_atable[i]->size();num++)
 145       for(pt=zone_atable[i]->begin();pt!=zone_atable[i]->end();pt++)
 146         if(pt->node2==v1&&pt->flag1==0&&pt->flag2==0)
 147         {
 148
 149           v1=pt->node1;         pt->flag1 =1;   pt->flag2 =1;
 150           if(v1!=v2) clockw.push(v1);
 151
 152         }
 153     //count clock wise search, push into a queue
 154     pt =  zone_atable[i]->begin();
 155     pt->flag2 = 1;
 156     cclockw.push(v2);
 157
 158     for(int num=0;num<zone_atable[i]->size();num++)
 159       for(pt=zone_atable[i]->begin();pt!=zone_atable[i]->end();pt++)
 160         if(pt->node1==v2&&pt->flag1==0&&pt->flag2==0)
 161         {
 162           v2=pt->node2;         pt->flag1 =1;   pt->flag2 =1;
 163           cclockw.push(v2);
 164         }
 165     //alloc memory for each vcell, set pointers,cp bc_index form danode
 166     VoronoiCell *pCell = davcell.GetPointer(i);
 167     int clockwsize = clockw.size();
 168     int cclockwsize = cclockw.size();
 169     int neighbor = clockw.size() + cclockw.size();
 170
 171     pCell->nb_num = neighbor;
 172     pCell->bc_index = danode[i].bc_index;
 173     pCell->nb_array = new int[neighbor];
 174     pCell->elen  = new double[neighbor];
 175     pCell->ilen  = new double[neighbor];
 176     pCell->angle = new double[neighbor];
 177     pCell->celledge = new int[neighbor];
 178     pCell->inb_array = new int[neighbor];
 179
 180     for(int j=0;j<pCell->nb_num;j++)  pCell->celledge[j] = -1;
 181
 182     //now read the neighbors of this node in count clock order
 183     for(int j=0;j<clockwsize;j++)
 184     {
 185       pCell->nb_array[j] = clockw.top();
 186       clockw.pop();
 187     }
 188     for(int j=clockwsize;j<cclockwsize+clockwsize;j++)
 189     {
 190       pCell->nb_array[j] = cclockw.front();
 191       cclockw.pop();
 192     }
 193
 194     adjtabsize += neighbor;
 195
 196   }
 197
 198   //adjtabsize equal to vedge*2
 199   davedge.clear();
 200   davedge.resize(adjtabsize/2);
 201
 202
 203   //compute invert index of neighbour nodes
 204   //i == davcell[j].nb_array[davcell[i].inb_array[j]]
 205   for(int i=0;i<davcell.size();i++)
 206   {
 207     VoronoiCell *pCell = davcell.GetPointer(i);
 208     for(int j=0;j<davcell[i].nb_num;j++)
 209       for(int k=0;k<davcell[davcell[i].nb_array[j]].nb_num;k++)
 210       {
 211         if(i==davcell[davcell[i].nb_array[j]].nb_array[k])
 212           pCell->inb_array[j]=k;
 213       }
 214   }
 215
 216   //not needed any more, free them.
 217   for(int i=0;i<davcell.size();i++)
 218     delete zone_atable[i];
 219   delete [] zone_atable ;
 220   return 0;
 221
 222 }
 223
 224 /* ----------------------------------------------------------------------------
 225  * ZONE::setup_voronoicell:  This function setup voronoicell data structure
 226  * the adjacency infomation will be build here.
 227  * if the Triangle is not count clock wise ordered, but the zone is rectangular.
 228  * please use this  function.
 229  */
 230 /*
 231 int ZONE::setup_voronoicell()
 232 {
 233         int adjtabsize=0;
 234         list<int>::iterator pt;
 235         list<int>**pneighborlist;
 236         davcell.clear();
 237         davcell.Init(danode.size());   //init davcell, cell num = node num
 238
 239         //compute max triangle, it may be useful
 240         double max_angle = 0;
 241         for(int i=0;i<datri.size();i++)
 242         {
 243                  double angle= MaxTriAngle(danode[datri[i].A].x,danode[datri[i].A].y,
 244                                            danode[datri[i].B].x,danode[datri[i].B].y,
 245                                            danode[datri[i].C].x,danode[datri[i].C].y);
 246                  if(max_angle<angle) max_angle=angle;
 247                  if(datri[i].D!=-1)
 248                         angle= MaxTriAngle(danode[datri[i].A].x,danode[datri[i].A].y,
 249                                            danode[datri[i].C].x,danode[datri[i].C].y,
 250                                            danode[datri[i].D].x,danode[datri[i].D].y);
 251                  if(max_angle<angle) max_angle=angle;
 252         }
 253         //alloc list. each voronoi cell will have it's own neighborlist.
 254         //my god ,tens of thounds of lists i should alloc.
 255         //stupid code
 256
 257         pneighborlist=new list<int>* [davcell.size()];
 258         for(int i=0;i<davcell.size();i++)       pneighborlist[i] = new list<int>;
 259         //first for each cell, insert all neighbor nodes into the each node's neighborlist.
 260         for(int i=0;i<datri.size();i++)
 261         {
 262
 263                 //triangle
 264                 pneighborlist[datri[i].A]->push_front(datri[i].C);
 265                 pneighborlist[datri[i].A]->push_front(datri[i].B);
 266                 pneighborlist[datri[i].B]->push_front(datri[i].A);
 267                 pneighborlist[datri[i].B]->push_front(datri[i].C);
 268                 pneighborlist[datri[i].C]->push_front(datri[i].B);
 269                 pneighborlist[datri[i].C]->push_front(datri[i].A);
 270
 271         }
 272
 273         //then remove the overlap nodes in pneighborlist
 274         for(int i=0;i<davcell.size();i++)
 275           for(pt=pneighborlist[i]->begin();pt!=pneighborlist[i]->end();)
 276             if(*pt==i||count(pneighborlist[i]->begin(),pneighborlist[i]->end(),*pt)>1 )
 277                 pt = pneighborlist[i]->erase(pt);
 278             else
 279                 pt++;
 280
 281         //alloc memory for each vcell, set pointers,cp bc_index form danode
 282         VoronoiCell *pCell = davcell.GetPointer(0);
 283         for(int i=0;i<davcell.size();i++,pCell++)
 284         {
 285                 pCell->nb_num = pneighborlist[i]->size();
 286                 pCell->nb_array = new int[pneighborlist[i]->size()];
 287                 pCell->elen  = new double[pneighborlist[i]->size()];
 288                 pCell->ilen  = new double[pneighborlist[i]->size()];
 289                 pCell->angle = new double[pneighborlist[i]->size()];
 290                 pCell->celledge = new int[pneighborlist[i]->size()];
 291                 pCell->inb_array = new int[pneighborlist[i]->size()];
 292
 293                 adjtabsize += pneighborlist[i]->size();
 294         }
 295         //adjtabsize equal to vedge*2
 296         davedge.clear();
 297         davedge.resize(adjtabsize/2);
 298
 299         // use angle to horizontal line to reorder the neighbour nodes as count clocked
 300         double x1,y1,x2,y2,x3,y3,r,angle;
 301         map<double,int>  adjmap;
 302         map<double,int> ::iterator padjmap;
 303         list<int> NODElist;
 304         list<int>::iterator pNODElist, pstart;
 305
 306         for(int i=0;i<davcell.size();i++)
 307         {       //compute each angle to horizontal line
 308                 //use container map to order them
 309                 for(pt=pneighborlist[i]->begin();pt!=pneighborlist[i]->end();pt++)
 310                 {
 311                         x2 = danode[i].x;               y2 = danode[i].y;
 312                         x1 = danode[*pt].x;             y1 = danode[*pt].y;
 313                         r  = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 314                         x3 = danode[i].x+r;             y3 = danode[i].y;
 315                         angle = CompAngle(x1,y1,x2,y2,x3,y3);
 316                         adjmap.insert(map<double,int>::value_type(angle,*pt));
 317                 }
 318                 int    *pnb_index = davcell.GetPointer(i)->nb_array;
 319                 double *pangle = davcell.GetPointer(i)->angle;
 320
 321                 for(padjmap = adjmap.begin();padjmap!=adjmap.end();padjmap++)
 322                 {
 323                         *pnb_index++ = padjmap->second;
 324                         *pangle++ = padjmap->first;
 325                 }
 326                 adjmap.clear();
 327
 328                 //boundary cell may have problem. the nb_array may be dissevered by bd nodes
 329                 //i have to reorder it to make sure that bd nodes are located at first and last
 330                 //positions of nb_array array
 331                 int do_reorder = 0;
 332                 if( davcell[i].bc_index != 0) //it is a bd cell
 333                 {       for(int j=0;j<davcell[i].nb_num-1;j++)
 334                            if(fabs(davcell[i].angle[j]-davcell[i].angle[j+1])>max_angle)
 335                                  {do_reorder=j+1;break;}
 336
 337                         // the angle between line i,j and line i,j+1 are larger than max_angle
 338                         // denote that j and j+1 are bd nodes.
 339                         // j+1 will be set at the beginning of nb_array ,j set as last value
 340                 }
 341                 if(do_reorder)
 342                 {
 343
 344                         //push all the neighbors into a list
 345                         for(int j=0;j<davcell[i].nb_num;j++)
 346                                 NODElist.push_back(davcell[i].nb_array[j]);
 347                         pNODElist=NODElist.begin();
 348                         while((*pNODElist)!=davcell[i].nb_array[do_reorder]) pNODElist++;
 349                         //find do_reorder point
 350                         pstart = pNODElist;
 351                         pnb_index = davcell.GetPointer(i)->nb_array;
 352                         //change the order
 353                         for(;pNODElist!=NODElist.end();pNODElist++)
 354                                 *pnb_index++ = *pNODElist;
 355                         for(pNODElist=NODElist.begin();pNODElist!=pstart;pNODElist++)
 356                                 *pnb_index++ = *pNODElist;
 357                         //recoder finished
 358                         NODElist.clear();
 359                         do_reorder=0;
 360                         //recompute angle
 361                         pangle = davcell.GetPointer(i)->angle;
 362                         for(int j=0;j<davcell[i].nb_num;j++)
 363                         {
 364                                 x2 = danode[i].x;               y2 = danode[i].y;
 365                                 Node nb_node = danode[davcell[i].nb_array[j]];
 366                                 x1 = nb_node.x;
 367                                 y1 = nb_node.y;
 368                                 r  = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 369                                 x3 = x2+r;
 370                                 y3 = y2;
 371                                 *pangle++ = CompAngle(x1,y1,x2,y2,x3,y3);
 372                         }
 373                  }
 374
 375
 376         }
 377         //compute invert index of neighbour nodes
 378         //davcell[i] == davcell[davcell[j].nb_array[davcell[i].inb_array[j]]]
 379         for(int i=0;i<davcell.size();i++)
 380         {
 381                 VoronoiCell *pCell = davcell.GetPointer(i);
 382                 for(int j=0;j<davcell[i].nb_num;j++)
 383                   for(int k=0;k<davcell[davcell[i].nb_array[j]].nb_num;k++)
 384                   {
 385                         if(i==davcell[davcell[i].nb_array[j]].nb_array[k])
 386                                 pCell->inb_array[j]=k;
 387                   }
 388         }
 389
 390         //not needed any more, free them.
 391         for(int i=0;i<davcell.size();i++)
 392                 delete pneighborlist[i];
 393         delete [] pneighborlist ;
 394         return 0;
 395
 396 }
 397 */
 398
 399
 400 /* ----------------------------------------------------------------------------
 401  * ZONE::setup_voronoicell: the remain work is eary to do. ilen,elen and area
 402  * will be setted here.
 403  */
 404 int ZONE::build_voronoicell_data()
 405 {
 406   //copy some data from danode
 407   for(int i=0;i<davcell.size();i++)
 408   {
 409     VoronoiCell *pCell = davcell.GetPointer(i);
 410     pCell->x = danode[i].x;
 411     pCell->y = danode[i].y;
 412     pCell->bc_index = danode[i].bc_index;
 413   }
 414
 415   //comput angle: the angle from horizontal line to it's neighbor node
 416   double x1,y1,x2,y2,x3,y3,r,angle;
 417   for(int i=0;i<davcell.size();i++)
 418   {
 419     double *pangle = davcell.GetPointer(i)->angle;
 420     for(int j=0;j<davcell[i].nb_num;j++)
 421     {
 422       x2 = danode[i].x;         y2 = danode[i].y;
 423       Node nb_node = danode[davcell[i].nb_array[j]];
 424       x1 = nb_node.x;
 425       y1 = nb_node.y;
 426       r  = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 427       x3 = x2+r;
 428       y3 = y2;
 429       *pangle++ = CompAngle(x1,y1,x2,y2,x3,y3);
 430     }
 431   }
 432
 433   //comput ilen: length between node i and it's neighbor
 434   for(int i=0;i<davcell.size();i++)
 435   {
 436     VoronoiCell *pCell = davcell.GetPointer(i);
 437     for(int j=0;j<davcell[i].nb_num;j++)
 438     {
 439       Node n = danode[davcell[i].nb_array[j]];
 440       pCell->ilen[j]=Distance(danode[i].x,danode[i].y,n.x,n.y);
 441     }
 442   }
 443
 444   //next to comput elen , length of Voronoi Cell's Edge
 445   //at the same time ,davedge will be set
 446   list<double> listxc,listyc; //two list to store location of CircleCenter
 447   list<double>::iterator pxc,pyc;
 448   double xc,yc,xc_next,yc_next;
 449   vector<VEdge>:: iterator pEdge = davedge.begin();
 450   int nedge = 0;
 451   for(int i=0;i<davcell.size();i++)
 452   {
 453     //fill the list with CircleCenter
 454     for(int j=0;j<davcell[i].nb_num-1;j++)
 455     {
 456       int error=CircleCenter(danode[i].x,danode[i].y,
 457                    danode[davcell[i].nb_array[j]].x,
 458                    danode[davcell[i].nb_array[j]].y,
 459                    danode[davcell[i].nb_array[j+1]].x,
 460                    danode[davcell[i].nb_array[j+1]].y,
 461                    xc,yc);
 462       if(error)
 463       {
 464         sprintf(log_buf,"warning! CircleCenter: the three nodes of vcell %d of zone %d coplanar.\n",i,zone_index);
 465         GSS_LOG();
 466       }
 467       listxc.push_back(xc);
 468       listyc.push_back(yc);
 469     }
 470     int nb_first= davcell[i].nb_array[0];
 471     int nb_last = davcell[i].nb_array[davcell[i].nb_num-1];
 472     //if all the 3 points are bd ndoes
 473     if (davcell[i].bc_index&&davcell[nb_first].bc_index&&
 474         davcell[nb_last].bc_index)
 475     {
 476       xc = (danode[i].x+danode[nb_first].x)/2.0;
 477       yc = (danode[i].y+danode[nb_first].y)/2.0;
 478       listxc.push_front(xc);
 479       listyc.push_front(yc);
 480       xc = (danode[i].x+danode[nb_last].x)/2.0;
 481       yc = (danode[i].y+danode[nb_last].y)/2.0;
 482       listxc.push_back(xc);
 483       listyc.push_back(yc);
 484     }
 485     else //compute CircleCenter of the last triangle
 486     {
 487       int error=CircleCenter(danode[i].x,danode[i].y,
 488                    danode[nb_first].x,
 489                    danode[nb_first].y,
 490                    danode[nb_last].x,
 491                    danode[nb_last].y,
 492                    xc,yc);
 493       if(error)
 494       {
 495         sprintf(log_buf,"warning! CircleCenter: the three nodes of vcell %d of zone %d coplanar.\n",i,zone_index);
 496         GSS_LOG();
 497       }
 498       listxc.push_front(xc);
 499       listyc.push_front(yc);
 500       listxc.push_back(xc);
 501       listyc.push_back(yc);
 502     }
 503
 504     pxc=listxc.begin();
 505     pyc=listyc.begin();
 506     VoronoiCell *pCell = davcell.GetPointer(i);
 507
 508     for(int j=0;j<davcell[i].nb_num;j++)
 509     {
 510       xc = *pxc++;       yc = *pyc++;
 511       xc_next = *pxc;    yc_next = *pyc;
 512       pCell->elen[j] = Distance(xc,yc,xc_next,yc_next);
 513
 514       // each vedge belongs to 2 vcells
 515
 516       int i_index = davcell[i].inb_array[j];
 517       if(davcell[i].celledge[j]==-1 && davcell[davcell[i].nb_array[j]].celledge[i_index]==-1)
 518       {       //this edge not registered. set data to it
 519         pEdge->cell1 = i;
 520         pEdge->cell2 = davcell[i].nb_array[j];
 521         pEdge->c1_nb = j;
 522         pEdge->c2_nb = davcell[i].inb_array[j];
 523
 524         pEdge->length = pCell->elen[j];
 525         pEdge->x1 = xc;       pEdge->y1 = yc;
 526         pEdge->x2 = xc_next;  pEdge->y2 = yc_next;
 527         davcell.GetPointer(i)->celledge[j] = nedge;
 528         davcell.GetPointer(davcell[i].nb_array[j])->celledge[i_index] = nedge;
 529         pEdge++; nedge++;
 530       }
 531     }
 532
 533     listxc.clear();
 534     listyc.clear();
 535   }
 536
 537   // compute each cell's volume and zone area
 538   field_area_vcell = 0;
 539   for(int i=0;i<davcell.size();i++)
 540   {
 541     VoronoiCell *pCell = davcell.GetPointer(i);
 542     pCell->area = 0;
 543     for(int j=0;j<davcell[i].nb_num;j++)
 544       pCell->area += davcell[i].elen[j]*davcell[i].ilen[j]/4.0;
 545     field_area_vcell += pCell->area;
 546   }
 547   return 0;
 548
 549 }