wzsxgrid.cxx

Delaunay三角形的网格剖分程序

	  sx->nlines=1;
	  sx->priv->FVMFactors=FVM2Edge; 
	  sx->priv->FEMFactors=FEM2Edge; 
          sx->priv->sxNext=sxNextFace2;
	  sx->priv->faceiter2= new DListIter<EDG2>(sx->grid->mesh->castToMESH2()->lines());
	  break;
  	case 3: 
	  sx->nlines=3;
	  sx->priv->FVMFactors=FVM3Triangle; 
	  sx->priv->FEMFactors=FEM3Triangle; 
          sx->priv->sxNext=sxNextFace3;
	  sx->priv->faceiter3= new DListIter<TR3>(sx->grid->mesh->castToMESH3()->triangles());
	  break;
	}
      sx->dim=g->spacedim-1;
      sx->codim=1;
      break;
      

    case SX_EDGES:
      sx->nnodes=2;
      sx->nlines=1;
      sx->linenodes=0;

      switch(sx->spacedim)
	{
	case 1: 
	  sx->priv->FVMFactors=0; 
	  sx->priv->FEMFactors=0;
	  sx->priv->sxNext=sxNextLine1;
	  sx->priv->lineiter1= new DListIter<EDG1>(sx->grid->mesh->castToMESH1()->lines());
	  break;
	case 2: 
	  sx->priv->FVMFactors=0;
	  sx->priv->FEMFactors=0; 
	  sx->priv->sxNext=sxNextLine2;
	  sx->priv->lineiter2= new DListIter<EDG2>(sx->grid->mesh->castToMESH2()->lines());
	  break;
  	case 3: 
	  sx->priv->FVMFactors=0; 
	  sx->priv->FEMFactors=0;
	  sx->priv->sxNext=sxNextLine3;
	  sx->priv->lineiter3= new DListIter<EDG3>(sx->grid->mesh->castToMESH3()->lines());
	  break;
	}
      sx->dim=1;
      sx->codim=g->spacedim-1;
      break;
      
    case SX_NODES:
      sx->nnodes=1; /* hier koennten Nachbarknoten/kanten stehen */
      sx->nlines=0;
      sx->linenodes=NULL;
      sx->node=node_node_ptr;
      
      sx->dim=0;
      sx->codim=g->spacedim;
      switch(sx->spacedim)
	{
	case 1: 
          sx->priv->sxNext=sxNextNode1;
	  sx->priv->nodeiter1= new DListIter<PT1>(sx->grid->mesh->castToMESH1()->points());
	  break;
	case 2: 
          sx->priv->sxNext=sxNextNode2;
	  sx->priv->nodeiter2= new DListIter<PT2>(sx->grid->mesh->castToMESH2()->points());
	  break;
  	case 3: 
          sx->priv->sxNext=sxNextNode3;
	  sx->priv->nodeiter3= new DListIter<PT3>(sx->grid->mesh->castToMESH3()->points());
	  break;
	}
      sx->priv->FVMFactors=FVM0Vertex;
      sx->priv->FEMFactors=FEM0Vertex;
      break;

    default:
      dprintf(err)("sorry, loop number %d not implemented\n",WhatToTraverse);
    }

   {
    int i,j;
    sx->stiffmat[0]=0;
    sx->massmat[0]=0;
    for (i=1,j=1;i<=sx->nnodes;i++,j+=sx->nnodes)
      sx->stiffmat[i]=&sx->priv->stiffmat[j]-1,
	sx->massmat[i]=&sx->priv->massmat[j]-1;
  }
 for(in=1;in<=sx->nnodes;in++)
    sx->coord[in]=sx->priv->coord+(in-1)*sx->spacedim;
  
  sx->priv->sxNext(sx);
}

int sxStop(Simplex sx)
{
  return (sx->priv==0);
}




////////////////////////////////////////////////////
// various data

void* sxMesh(sxGrid g)  { return g->mesh;     }
int sxDimension(sxGrid g) { return g->spacedim; }
int sxNumber(sxGrid g, int comp) 
{ 
  switch (comp)
  {
  case SX_NODES: return g->node_num; 
  case SX_CELLS: return g->cell_num; 
  case SX_FACES: return g->face_num; 
  case SX_EDGES: return g->line_num; 
  case SX_ALL_FACES: return g->all_face_num; 
  default: return 0;
  }
}

int sxNumberOfMaterials(sxGrid g, int comp) 
{ 
  switch (comp)
  {
  case SX_NODES: return g->node_mat_num; 
  case SX_CELLS: return g->cell_mat_num; 
  case SX_FACES: return g->face_mat_num; 
  case SX_EDGES: return 0; 
  default: return 0;
  }
}


double sxXMin(sxGrid g) { return g->xmin; }
double sxXMax(sxGrid g) { return g->xmax; }
double sxYMin(sxGrid g) { return g->ymin; }
double sxYMax(sxGrid g) { return g->ymax; }
double sxZMin(sxGrid g) { return g->zmin; }
double sxZMax(sxGrid g) { return g->zmax; }


CmdPars  Cmd;				// command Parameters
//extern - on SGI initialization seems to work only in main prog! 

    
////////////////////////////////////////////////////
// Create/Destroy
///////////////////////////////////////////////////


void sxgPrintData(sxGrid g) // same for creation and refinement!
{
  int i;
  int bad;
  double *cell_mat_volume;
  SimplexData sx;
  int bad_cell_count_zero_vol;
  int bad_cell_count_neg_line;
  double bad_cell_volume;
  
  bad_cell_count_zero_vol=0;
  bad_cell_count_neg_line=0;
  bad_cell_volume=0;

  cell_mat_volume=(double*)calloc(g->cell_mat_num+1,sizeof(double));
  SX_LOOP(g,sx,SX_CELLS)
    {
      bad=0;
      sxFVMFactors(&sx);
      cell_mat_volume[0]+=sx.volume;
      cell_mat_volume[sx.material]+=sx.volume;
      if (sx.volume==0.0) 
	{
	  bad_cell_count_zero_vol++;
	  bad=1;
	}
      for (i=1;i<=sx.nlines;i++)
	{
	  if (sx.linefac[i]<0.0) 
	    {
	      bad_cell_count_neg_line++;
	      bad=1;
	      break;
	    }
	}
      if (bad) bad_cell_volume+=sx.volume;
    }
  if (debug(inf))
  printf("-----------------------------------------------------------------\n");
  dprintf(inf)("cell data:\n");
  dprintf(inf)("grid:       %4d cells, vol=%g\n",g->cell_num,cell_mat_volume[0]);
  for(i=1;i<=g->cell_mat_num;i++)
    if(g->cell_mat_cell_num[i])
      dprintf(inf)("material %d: %4d cells, vol=%g\n",
		   i,
		   g->cell_mat_cell_num[i],cell_mat_volume[i]);
      dprintf(inf)("bad cells:   zero vol: %d,  neg line: %d, vol=%g\n",
		   bad_cell_count_zero_vol,
		   bad_cell_count_neg_line,
		   bad_cell_volume);
  free(cell_mat_volume);

  if (debug(inf))
  printf("-----------------------------------------------------------------\n");
  dprintf(inf)("face data:\n");
  for(i=1;i<=g->face_mat_num;i++)
    if(g->face_mat_face_num[i])
	  dprintf(inf)("material %d: %d faces\n",i,g->face_mat_face_num[i]);
  dprintf(inf)("boundary faces: %d\n",g->face_num);
  dprintf(inf)("all faces: %d\n",g->all_face_num);
  
  
  if (debug(inf))
  printf("-----------------------------------------------------------------\n");
  dprintf(inf)("node data:\n");
  dprintf(inf)("nodes:    %d\n",g->node_num);
  dprintf(inf)("spacedim: %d\n",g->spacedim);
  dprintf(inf)("xcoord: [ %g , %g ]\n",g->xmin,g->xmax);
  dprintf(inf)("ycoord: [ %g , %g ]\n",g->ymin,g->ymax);
  if (g->spacedim>2)
    dprintf(inf)("zcoord: [ %g , %g ]\n",g->zmin,g->zmax);
  
  if (debug(inf))
  printf("-----------------------------------------------------------------\n");
  
}


static int bisect=0;
void sxSetBisection(void)
{
  bisect=1;
}


void number_stuff(sxGrid g)
{

  MESH *m=g->mesh;
  register int i;
// sorry for this!
  switch (g->spacedim)
    {
    case 1:
      {
	PT1 *patch;
	DListIter<PT1> iter(m->castToMESH1()->points());
	i=0;
	while ((patch = iter.all())) 
	  {
	    i++;
	    patch->setNode(i);
	  }
	g->node_num=i;
      }
      break;
    case 2:
      {
	PT2 *patch;
	DListIter<PT2> iter(m->castToMESH2()->points());
	i=0;
	while ((patch = iter.all())) 
	  {
	    i++;
	    patch->setNode(i);
	  }
	g->node_num=i;
      }
      break;
    case 3:
      {
	PT3 *patch;
	DListIter<PT3> iter(m->castToMESH3()->points());
	i=0;
	while ((patch =iter.all())) 
	  {
	    i++;
	    patch->setNode(i++);
	  }
	g->node_num=i;
      }
      break;
      
    }
    
  PATCH *patch;
  m->resetElemIter();
  for(i=1;i<=g->cell_mat_num;i++) g->cell_mat_cell_num[i]=0;
  i=0;
  while ((patch = m->elemIterAll())) 
    { 
      i++;
      patch->setNode(i);
      g->cell_mat_cell_num[patch->Class()]++;
    }

  
  
  switch(g->spacedim)
    { 
    case 1:
      {
	EDG1* patch;
	DListIter<EDG1>  iter(m->castToMESH1()->lines());
	i=0;
	while ((patch = iter.all()))     
	  {
	    i++;
	    patch->setNode(i);
	  }
	g->line_num=i;
      }
      break;
    case 2:
      {
	EDG2* patch;
	DListIter<EDG2>  iter(m->castToMESH2()->lines());
	i=0;
	while ((patch = iter.all()))     
	  {
	    i++;
	    patch->setNode(i);
	  }
	g->line_num=i;
      }
      break;
     case 3:
       {
	EDG3* patch;
	DListIter<EDG3>  iter(m->castToMESH3()->lines());
	i=0;
	while ((patch = iter.all()))     
	  {
	    i++;
	    patch->setNode(i);
	  }
	g->line_num=i;
      }
      break;
    }
  
  
  for(i=1;i<=g->face_mat_num;i++) g->face_mat_face_num[i]=0;
  switch(g->spacedim)
    { 
    case 1:
      {
	PT1* patch;
	DListIter<PT1>  iter(m->castToMESH1()->points());
	while ((patch = iter.all()))     
	  if (patch->onBoundary()&&patch->Class())
	    g->face_mat_face_num[patch->Class()]++;

	g->all_face_num=g->line_num;

      }
      break;
    case 2:
      {
	EDG2* patch;
	DListIter<EDG2>  iter(m->castToMESH2()->lines());
	while ((patch = iter.all()))     
	  if (patch->onBoundary()&&patch->Class())
	    g->face_mat_face_num[patch->Class()]++;
	g->all_face_num=g->line_num;
      }
      break;
    case 3:
      {
	TR3* patch;
	DListIter<TR3>  iter(m->castToMESH3()->triangles());;
	iter.reset();
	int i=0;
	while ((patch = iter.all()))     
	{
	  i++;
	  patch->setNode(i);
	  if (patch->onBoundary()&&patch->Class())
	    g->face_mat_face_num[patch->Class()]++;
	}
	g->all_face_num=i;
      }
      break;
    }
  g->face_num=0;
  for(i=1;i<=g->face_mat_num;i++) g->face_num+=g->face_mat_face_num[i];
  g->cell_num=g->mesh->noOfElements();
  /*
  {
  
    void find_face_in_list(n1,n2,n3)
      {
	sort(n1,n2,n3);
	
      }


  }

  */
}


sxGrid sxCreate(char *meshname, void *user_data)
{

  register int i;
  sxGrid g;
  char fname[64];
  MESH *m=NULL;
  int spacedim=0;
 
  init();
  Cmd.ReadCmdFile(meshname);  
  
///////////////////////////////////////
///// MESH creation
  if (Cmd.isTrue("printParameters")) cout << Cmd;
  
   if (spacedim==0)
    { 
      if (!Cmd.get("spaceDim",&spacedim)) 
	{
	  dprintf(wrn)("Spacedim said to be 0, assuming 2!\n");
	  spacedim=2;
	}
    }
  
  sprintf(fname,"%s.geo",meshname);
  switch (spacedim)
    {
    case 1: m=new MESH1(fname);
	    break;
    case 2: 
	    if (bisect)
	      {
		m=new MESH2BISECT(fname);
	      }
	    else 
	      m=new MESH2(fname); 
	    break;
    case 3: m=new MESH3(fname);
	    break;
    }
  
  assert(m!=NULL);



  g=new sxGridData;
  assert(g!=NULL);
  
  g->mesh=m;
  g->spacedim=spacedim;

  if (bisect) 
    {
      g->bisect=bisect;
      bisect=0;
    }
///////////////////////////////////////
///// Set Simplex format data

  switch (g->spacedim)
    {
    case 1:
      g->cell_line_num=1;
      g->cell_face_num=2;
      g->cell_line_node_ptr[1][1]=1;
      g->cell_line_node_ptr[1][2]=2;
   
      g->cell_face_node_ptr[1][1]=1;
      g->cell_face_node_ptr[1][2]=2;
 
      g->cell_face_node_ptr[2][1]=2;
      g->cell_face_node_ptr[2][2]=1;
      break;
    case 2:
      g->cell_line_num=3;
      g->cell_face_num=3;
      g->cell_line_node_ptr[1][1]=2;
      g->cell_line_node_ptr[1][2]=3;
      g->cell_line_node_ptr[1][3]=1;

      g->cell_line_node_ptr[2][1]=1;
      g->cell_line_node_ptr[2][2]=3;
      g->cell_line_node_ptr[2][3]=2;

      g->cell_line_node_ptr[3][1]=1;
      g->cell_line_node_ptr[3][2]=2;
      g->cell_line_node_ptr[3][3]=3;
 
      g->cell_face_node_ptr[1][1]=2;
      g->cell_face_node_ptr[1][2]=3;
      g->cell_face_node_ptr[1][3]=1;

      g->cell_face_node_ptr[2][1]=1;
      g->cell_face_node_ptr[2][2]=3;
      g->cell_face_node_ptr[2][3]=2;

      g->cell_face_node_ptr[3][1]=1;
      g->cell_face_node_ptr[3][2]=2;
      g->cell_face_node_ptr[3][3]=3;
      break;
    case 3:
      g->cell_line_num=6;
      g->cell_face_num=4;
  /* die ersten beiden muessen mit celltype konform gehen !*/
      g->cell_line_node_ptr[1][1]=3;
      g->cell_line_node_ptr[1][2]=4;
      g->cell_line_node_ptr[1][3]=1;
      g->cell_line_node_ptr[1][4]=2;


      g->cell_line_node_ptr[2][1]=4;
      g->cell_line_node_ptr[2][2]=2;
      g->cell_line_node_ptr[2][3]=1;
      g->cell_line_node_ptr[2][4]=3;


      g->cell_line_node_ptr[3][1]=2;
      g->cell_line_node_ptr[3][2]=3;
      g->cell_line_node_ptr[3][3]=4;
      g->cell_line_node_ptr[3][4]=1;

      g->cell_line_node_ptr[4][1]=1;
      g->cell_line_node_ptr[4][2]=2;
      g->cell_line_node_ptr[4][3]=3;
      g->cell_line_node_ptr[4][4]=4;

      g->cell_line_node_ptr[5][1]=1;
      g->cell_line_node_ptr[5][2]=3;
      g->cell_line_node_ptr[5][3]=2;
      g->cell_line_node_ptr[5][4]=4;

      g->cell_line_node_ptr[6][1]=1;
      g->cell_line_node_ptr[6][2]=4;
      g->cell_line_node_ptr[6][3]=2;
      g->cell_line_node_ptr[6][4]=3;



      g->cell_face_node_ptr[1][1]=2;
      g->cell_face_node_ptr[1][2]=3;
      g->cell_face_node_ptr[1][3]=4;
      g->cell_face_node_ptr[1][4]=1;

      g->cell_face_node_ptr[2][1]=1;
      g->cell_face_node_ptr[2][2]=3;
      g->cell_face_node_ptr[2][3]=4;
      g->cell_face_node_ptr[2][4]=2;

      g->cell_face_node_ptr[3][1]=1;
      g->cell_face_node_ptr[3][2]=2;
      g->cell_face_node_ptr[3][3]=4;
      g->cell_face_node_ptr[3][4]=3;

      g->cell_face_node_ptr[4][1]=1;
      g->cell_face_node_ptr[4][2]=2;
      g->cell_face_node_ptr[4][3]=3;
      g->cell_face_node_ptr[4][4]=4;


      break;
    }
  g->cell_node_num= g->spacedim+1;
  
///////////////////////////////////////
///// Calculate min/max and node numbers ?!
///  MESH*::noOfPoints is   "private (at least in principle)" !!
  g->prev_node_num=0;

  g->xmin=g->ymin=g->zmin=  1.0e30;
  g->xmax=g->ymax=g->zmax= -1.0e30;
  
  switch (spacedim)
    {
    case 1:
      {
	PT1 *patch;
	DListIter<PT1> iter(m->castToMESH1()->points());
	while ((patch = iter.all())) 
	  {
	    if (patch->x< g->xmin) g->xmin=patch->x;
	    if (patch->x> g->xmax) g->xmax=patch->x;
	  }
	g->ymin=g->zmin=0.0;
	g->ymax=g->zmax=0.0;
      }
      break;
    case 2:
      {