src/solver/refine.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: Dec 19, 2005                                                */
  14 /*                                                                           */
  15 /*  Gong Ding                                                                */
  16 /*  gdiso@ustc.edu                                                           */
  17 /*  NINT, No.69 P.O.Box, Xi'an City, China                                   */
  18 /*                                                                           */
  19 /*****************************************************************************/
  20
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <string.h>
  24 #include <cgnslib.h>
  25 #include <math.h>
  26 #include <iostream>
  27 #include "bsolver.h"
  28 #include "geom.h"
  29 #include "typedef.h"
  30 #define ANSI_DECLARATORS
  31
  32 extern "C"
  33 {
  34 #include "triangle.h"
  35 }
  36
  37 inline double FunLinear(PetscScalar value)
  38 {
  39   return double(value);
  40 }
  41
  42 inline double FunSignedLog(PetscScalar value)
  43 {
  44   return double(dsign(value)*log(1.0+fabs(value)));
  45 }
  46
  47 /* ----------------------------------------------------------------------------
  48  * BSolver::set_scale  This function compute the refine scale
  49  */
  50 double  BSolver::set_refine_scale(Tri & tri, RefineDefine & ref)
  51 {
  52   double M1=0,M2=0,M3=0,M,Scale=1.1;
  53   double (*Measure) (PetscScalar);
  54   int z = tri.zone_index;
  55
  56   if(zonedata[z]->material_type == Semiconductor)
  57   {
  58     if(ref.Measure==Linear)
  59       Measure = FunLinear;
  60     if(ref.Measure==SignedLog)
  61       Measure = FunSignedLog;
  62     SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
  63
  64     int p1 = tri.node[0];
  65     int p2 = tri.node[1];
  66     int p3 = tri.node[2];
  67
  68     if(ref.Variable == Potential)
  69     {
  70       M1 = Measure(pzonedata->fs[p1].P);
  71       M2 = Measure(pzonedata->fs[p2].P);
  72       M3 = Measure(pzonedata->fs[p3].P);
  73       M =  dmax(dmax(fabs(M1-M2),fabs(M2-M3)),fabs(M3-M1))+1e-6;
  74       if(M>ref.Dispersion)
  75         Scale = ref.Dispersion/M;
  76       if(Scale<ref.DivisionRatio)        Scale = ref.DivisionRatio;
  77       return Scale;
  78     }
  79     if(ref.Variable == Doping)
  80     {
  81       // the mesh size must small than Debye length.
  82       PetscScalar doping = fabs(pzonedata->aux[p1].Nd - pzonedata->aux[p1].Na
  83                                +pzonedata->aux[p2].Nd - pzonedata->aux[p2].Na
  84                                +pzonedata->aux[p3].Nd - pzonedata->aux[p3].Na)/3.0
  85                            +1.0*pow(scale_unit.s_centimeter,-3);
  86       PetscScalar  Ld = sqrt(pzonedata->aux[p1].eps*pzonedata->mt->kb*pzonedata->fs[p1].T/(pzonedata->mt->e*pzonedata->mt->e*doping));
  87
  88       M1 = Measure(pzonedata->aux[p1].Nd*pow(scale_unit.s_centimeter,3)-
  89                    pzonedata->aux[p1].Na*pow(scale_unit.s_centimeter,3));
  90       M2 = Measure(pzonedata->aux[p2].Nd*pow(scale_unit.s_centimeter,3)-
  91                    pzonedata->aux[p2].Na*pow(scale_unit.s_centimeter,3));
  92       M3 = Measure(pzonedata->aux[p3].Nd*pow(scale_unit.s_centimeter,3)-
  93                    pzonedata->aux[p3].Na*pow(scale_unit.s_centimeter,3));
  94       M =  dmax(dmax(fabs(M1-M2),fabs(M2-M3)),fabs(M3-M1))+1e-6;
  95       if(M>ref.Dispersion)
  96         Scale = Ld*Ld/(2*tri.area) < ref.Dispersion/M ? Ld*Ld/(2*tri.area):ref.Dispersion/M;
  97       //if(Scale>Ld/sqrt(2*tri.area))   Scale = Ld*Ld/(2*tri.area);
  98       if(Scale<ref.DivisionRatio)       Scale = ref.DivisionRatio;
  99       return Scale;
 100     }
 101
 102   }
 103   else if(zonedata[z]->material_type == Insulator)
 104   {
 105     return      1.1;
 106   }
 107   else if(zonedata[z]->material_type == Conductor)
 108   {
 109     return      1.1;
 110   }
 111   else if(zonedata[z]->material_type == Vacuum)
 112   {
 113     return      1.1;
 114   }
 115   else if(zonedata[z]->material_type == PML)
 116   {
 117     return      1.1;
 118   }
 119   return 0.0; //prevent warning of complier
 120 }
 121
 122 /* ----------------------------------------------------------------------------
 123  * BSolver::refine  This function call triangulate to refine the mesh
 124  *
 125  */
 126
 127 int BSolver::refine(RefineDefine & ref)
 128 {
 129   struct triangulateio in, out;
 130
 131   //do necessarily prepare for call triangulate
 132   in.numberofpoints = gnode.size();
 133   in.numberofpointattributes = 6; //P, n, p, Nd, Na and mole_x
 134
 135   in.pointattributelist = (double *)malloc(in.numberofpoints * in.numberofpointattributes * sizeof(double));
 136   in.pointlist = (double *) malloc(in.numberofpoints * 2 * sizeof(double));
 137   in.pointmarkerlist = (int *) malloc(in.numberofpoints * sizeof(int));
 138   for(int i=0;i<in.numberofpoints;i++)
 139   {
 140     in.pointlist[2*i] = gnode[i].x;
 141     in.pointlist[2*i+1] = gnode[i].y;
 142     in.pointmarkerlist[i] = gnode[i].bc_index;
 143   }
 144   double *ppointattributelist = in.pointattributelist;
 145 /*
 146   // set all the interface node in global node array as semiconductor node
 147   for(int z=0;z<zone_num;z++)
 148     if(zonedata[z]->material_type == Semiconductor)
 149     {
 150       for(int i=0;i<zone[z].danode.size();i++)
 151       {
 152         gnode[zone[z].danode[i].g_index].local_index = i;
 153         gnode[zone[z].danode[i].g_index].zone_index = z;
 154       }
 155     }
 156 */
 157   for(int i=0;i<in.numberofpoints;i++)
 158   {
 159     int z = gnode[i].zone_index;
 160     if(zonedata[z]->material_type == Semiconductor)
 161     {
 162       SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
 163       int n = gnode[i].local_index;
 164       //*ppointattributelist++ = log(1.0+pzonedata->fs[n].Nd*pow(scale_unit.s_centimeter,3));
 165       //*ppointattributelist++ = log(1.0+pzonedata->fs[n].Na*pow(scale_unit.s_centimeter,3));
 166       *ppointattributelist++ = double(pzonedata->fs[n].P);
 167       *ppointattributelist++ = double(pzonedata->fs[n].n*pow(scale_unit.s_centimeter,3));
 168       *ppointattributelist++ = double(pzonedata->fs[n].p*pow(scale_unit.s_centimeter,3));
 169       *ppointattributelist++ = double(pzonedata->aux[n].Nd*pow(scale_unit.s_centimeter,3));
 170       *ppointattributelist++ = double(pzonedata->aux[n].Na*pow(scale_unit.s_centimeter,3));
 171       *ppointattributelist++ = double(pzonedata->aux[n].mole_x);
 172     }
 173     else if(zonedata[z]->material_type == Insulator)
 174     {
 175       ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
 176       int n = gnode[i].local_index;
 177       *ppointattributelist++ = double(pzonedata->fs[n].P);
 178       *ppointattributelist++ = 0;
 179       *ppointattributelist++ = 0;
 180       *ppointattributelist++ = 0;
 181       *ppointattributelist++ = 0;
 182       *ppointattributelist++ = 0;
 183     }
 184     else
 185     {
 186       *ppointattributelist++ = 0;
 187       *ppointattributelist++ = 0;
 188       *ppointattributelist++ = 0;
 189       *ppointattributelist++ = 0;
 190       *ppointattributelist++ = 0;
 191       *ppointattributelist++ = 0;
 192     }
 193   }
 194
 195   in.numberoftriangles = gtri.size();
 196   in.numberofcorners = 3;
 197   in.numberoftriangleattributes = 1;
 198   in.trianglelist =  (int *) malloc(in.numberoftriangles * 3 * sizeof(int));
 199   in.trianglearealist = (double *) malloc(in.numberoftriangles * sizeof(double));
 200   in.triangleattributelist = (double *) malloc(in.numberoftriangles * sizeof(double));
 201   int *ptrianglelist = in.trianglelist;
 202   double *ptrianglearealist = in.trianglearealist;
 203   double *ptriangleattributelist = in.triangleattributelist;
 204   for(int i=0;i<gtri.size();i++)
 205   {
 206      *ptrianglelist++ = gtri[i].g_node[0];
 207      *ptrianglelist++ = gtri[i].g_node[1];
 208      *ptrianglelist++ = gtri[i].g_node[2];
 209      *ptriangleattributelist++ = gtri[i].zone_index;
 210   }
 211   for(int i=0;i<gtri.size();i++)
 212   {
 213     double Scale=0.5;
 214     Scale = set_refine_scale(gtri[i],ref);
 215     *ptrianglearealist++ = gtri[i].area*Scale;
 216   }
 217   in.numberofsegments =  0;
 218   for(int i=0;i<gsegment.size();i++)
 219     in.numberofsegments += gsegment[i].edge_num;
 220   in.segmentlist =  (int *) malloc(in.numberofsegments * 2 * sizeof(int));
 221   in.segmentmarkerlist = (int *) malloc(in.numberofsegments * sizeof(int));
 222   int *psegmentlist =  in.segmentlist;
 223   int *psegmentmarkerlist = in.segmentmarkerlist;
 224   for(int i=0;i<gsegment.size();i++)
 225     for(int j=0;j<gsegment[i].edge_num;j++)
 226     {
 227       *psegmentlist++ = gsegment[i].edge_array[j].p1;
 228       *psegmentlist++ = gsegment[i].edge_array[j].p2;
 229       *psegmentmarkerlist++ = gsegment[i].bc_index;
 230     }
 231
 232   in.numberofholes = 0;
 233   in.numberofregions = 0;
 234
 235
 236   // Make necessary initializations so that Triangle can return a
 237   // triangulation in `out'.
 238
 239   out.pointlist = (double *) NULL;
 240   out.pointattributelist = (double *) NULL;
 241   out.pointmarkerlist = (int *) NULL;
 242   out.trianglelist = (int *) NULL;
 243   out.triangleattributelist = (double *) NULL;
 244   out.segmentlist = (int *) NULL;
 245   out.segmentmarkerlist = (int *) NULL;
 246
 247   // Refine the triangulation according to the attached
 248   // triangle area constraints.
 249
 250   triangulate(ref.tri_cmd, &in, &out, (struct triangulateio *) NULL);
 251
 252   gnode.resize(out.numberofpoints);
 253   gtri.resize(out.numberoftriangles);
 254   for(int i=0; i<gsegment.size(); i++)
 255   {
 256     gsegment[i].edge_array.clear();
 257     gsegment[i].edge_num = 0;
 258   }
 259
 260   ppointattributelist = out.pointattributelist;
 261   for(int i=0;i<out.numberofpoints;i++)
 262   {
 263     gnode[i].x = out.pointlist[2*i];
 264     gnode[i].y = out.pointlist[2*i+1];
 265     gnode[i].bc_index = out.pointmarkerlist[i];
 266     gnode[i].g_index = i;
 267     gnode[i].zone_index = -1;
 268     gnode[i].local_index = -1;
 269   }
 270
 271   ptrianglelist = out.trianglelist;
 272   for(int i=0;i<out.numberoftriangles;i++)
 273   {
 274     gtri[i].zone_index = static_cast<int>( out.triangleattributelist[i]+0.5 );
 275     gtri[i].g_node[0] = *ptrianglelist++;
 276     gtri[i].g_node[1] = *ptrianglelist++;
 277     gtri[i].g_node[2] = *ptrianglelist++;
 278     gtri[i].bc[0] = 0;
 279     gtri[i].bc[1] = 0;
 280     gtri[i].bc[2] = 0;
 281   }
 282
 283   for(int j=0;j<gsegment.size();j++)
 284     for(int i=0;i<out.numberofsegments;i++)
 285       if(out.segmentmarkerlist[i] == gsegment[j].bc_index)
 286       {
 287         Edge tmp_edge;
 288         tmp_edge.p1 = out.segmentlist[2*i];
 289         tmp_edge.p2 = out.segmentlist[2*i+1];
 290         tmp_edge.bc_index = out.segmentmarkerlist[i];
 291         gsegment[j].edge_array.push_back(tmp_edge);
 292         gsegment[j].edge_num++;
 293       }
 294
 295   field_to_zone();
 296   build_zone();
 297
 298
 299   clear_data();
 300   bc.clear();
 301   if(setup_bc())   return 1;
 302   if(build_zonedata())    return 1;
 303
 304   //re-arrange data
 305   if(doping_func.size())
 306         setup_doping();
 307   else
 308   {
 309     ppointattributelist = out.pointattributelist;
 310     for(int z=0;z<zone_num;z++)
 311       if(zonedata[z]->material_type == Semiconductor)
 312     {
 313       SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
 314       for(int i=0;i<zone[z].danode.size();i++)
 315       {
 316         int index = zone[z].danode[i].g_index;
 317         pzonedata->aux[i].Nd = fabs((ppointattributelist[6*index+3]))/pow(scale_unit.s_centimeter,3);
 318         pzonedata->aux[i].Na = fabs((ppointattributelist[6*index+4]))/pow(scale_unit.s_centimeter,3);
 319       }
 320     }
 321   }
 322
 323   //mole function
 324   for(int z=0;z<zone_num;z++)
 325     if(IsSingleCompSemiconductor(zonedata[z]->pzone->zonelabel))
 326     {
 327       SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
 328       for(int i=0;i<zone[z].danode.size();i++)
 329       {
 330         int index = zone[z].danode[i].g_index;
 331         pzonedata->aux[i].mole_x = fabs(ppointattributelist[6*index+5]);
 332         if(pzonedata->aux[i].mole_x > 1.0) pzonedata->aux[i].mole_x = 1.0;
 333       }
 334     }
 335   setup_init_data();
 336 /*
 337   //The re-computing of physical values on new mesh even make things worse.
 338   //just use initial value.
 339   for(int z=0;z<zone_num;z++)
 340   {
 341     if(zonedata[z]->material_type == Semiconductor)
 342     {
 343       SMCZone *pzonedata = dynamic_cast< SMCZone * >(zonedata[z]);
 344       for(int i=0;i<zone[z].danode.size();i++)
 345       {
 346         int index = zone[z].danode[i].g_index;
 347         pzonedata->fs[i].P = ppointattributelist[6*index+0];
 348         pzonedata->fs[i].n = fabs((ppointattributelist[6*index+1]))/pow(scale_unit.s_centimeter,3)+1e-12;
 349         pzonedata->fs[i].p = fabs((ppointattributelist[6*index+2]))/pow(scale_unit.s_centimeter,3)+1e-12;
 350         //printf("%e %e %e\n",pzonedata->fs[i].P,pzonedata->fs[i].n,pzonedata->fs[i].p);
 351       }
 352     }
 353     else if(zonedata[z]->material_type == Insulator)
 354     {
 355       ISZone *pzonedata = dynamic_cast< ISZone * >(zonedata[z]);
 356       for(int i=0;i<zone[z].danode.size();i++)
 357       {
 358         int index = zone[z].danode[i].g_index;
 359         pzonedata->fs[i].P = fabs((ppointattributelist[6*index+0]));
 360       }
 361     }
 362   }
 363 */
 364   // Free all allocated arrays, including those allocated by Triangle.
 365
 366   free(in.pointlist);
 367   free(in.pointmarkerlist);
 368   free(in.pointattributelist);
 369   free(in.trianglelist);
 370   free(in.triangleattributelist);
 371   free(in.trianglearealist);
 372   free(in.segmentlist);
 373   free(in.segmentmarkerlist);
 374
 375   free(out.pointlist);
 376   free(out.pointmarkerlist);
 377   free(out.pointattributelist);
 378   free(out.trianglelist);
 379   free(out.triangleattributelist);
 380   free(out.segmentlist);
 381   free(out.segmentmarkerlist);
 382
 383   return 0;
 384 }