sxgridint.c

Delaunay三角形的网格剖分程序

#include <stdio.h>
#include <stdlib.h>
#define wzgrid_hxx
#include "sxgridint.h"
#include "memory.h"

static void sxCountMaterials(sxGrid g, sxIteratorType t)
{
  sxSimplexData sx;
  int nmat=0;
  if (t==sxEdges) return;
  if (t==sxFaces) return;
  sxLoop(g,&sx,t)
    {
      if (sx.material>nmat) nmat=sx.material;
    }
  g->numOfMaterials[t]=nmat;
}

static void sxCountSimplices(sxGrid g, sxIteratorType t)
{
  sxSimplexData sx;
  int nsimp=0;
  if (t==sxEdges) return;
  if (t==sxFaces) return;

  sxLoop(g,&sx,t)
    { 
      nsimp++;
    }
  g->numOfSimplices[t]=nsimp;
}

static void sxDefineBox(sxGrid g)
{
  sxSimplexData sx;
  g->mincoord[1]=g->mincoord[2]=g->mincoord[3]=  1.0e30;
  g->maxcoord[1]=g->maxcoord[2]=g->maxcoord[3]= -1.0e30;
  switch(g->spaceDim)
    {
    case 1:
      sxLoop(g,&sx,sxNodes)
	{
	  if(sx.coordinates[1][1]<g->mincoord[1]) g->mincoord[1]=sx.coordinates[1][1];
	  else if(sx.coordinates[1][1]>g->maxcoord[1]) g->maxcoord[1]=sx.coordinates[1][1];
	}
      
    case 2:
      sxLoop(g,&sx,sxNodes)
	{
	  if(sx.coordinates[1][1]<g->mincoord[1]) g->mincoord[1]=sx.coordinates[1][1];
	  else if(sx.coordinates[1][1]>g->maxcoord[1]) g->maxcoord[1]=sx.coordinates[1][1];

	  if(sx.coordinates[1][2]<g->mincoord[2]) g->mincoord[2]=sx.coordinates[1][2];
	  else if(sx.coordinates[1][2]>g->maxcoord[2]) g->maxcoord[2]=sx.coordinates[1][2];
	}
      
    case 3:
      sxLoop(g,&sx,sxNodes)
	{
	  if(sx.coordinates[1][1]<g->mincoord[1]) g->mincoord[1]=sx.coordinates[1][1];
	  else if(sx.coordinates[1][1]>g->maxcoord[1]) g->maxcoord[1]=sx.coordinates[1][1];

	  if(sx.coordinates[1][2]<g->mincoord[2]) g->mincoord[2]=sx.coordinates[1][2];
	  else if(sx.coordinates[1][2]>g->maxcoord[2]) g->maxcoord[2]=sx.coordinates[1][2];

	  if(sx.coordinates[1][3]<g->mincoord[3]) g->mincoord[3]=sx.coordinates[1][3];
	  else if(sx.coordinates[1][3]>g->maxcoord[3]) g->maxcoord[3]=sx.coordinates[1][3];
	}
    }
  g->boxNotDefined=0;
}

double sxXMin(sxGrid g)
{
  if(g->boxNotDefined) sxDefineBox(g); 
  return g->mincoord[1];
}

double sxXMax(sxGrid g)
{
  if(g->boxNotDefined) sxDefineBox(g); 
  return g->maxcoord[1];
}

double sxYMin(sxGrid g)
{
  if(g->boxNotDefined) sxDefineBox(g); 
  return g->mincoord[2];
}

double sxYMax(sxGrid g)
{
  if(g->boxNotDefined) sxDefineBox(g); 
  return g->maxcoord[2];
}

double sxZMin(sxGrid g)
{
  if(g->boxNotDefined) sxDefineBox(g); 
  return g->mincoord[3];
}

double sxZMax(sxGrid g)
{
  if(g->boxNotDefined) sxDefineBox(g); 
  return g->maxcoord[3];
}

int sxDimension(sxGrid g)
{
  return g->gridDim;
}

int sxNumber(sxGrid g, sxIteratorType t)
{
  if(!g->numOfSimplices[t]) sxCountSimplices(g,t); 
  return g->numOfSimplices[t];
}

int sxNumberOfMaterials(sxGrid g, sxIteratorType t)
{
  if(!g->numOfMaterials[t]) sxCountMaterials(g,t); 
  return g->numOfMaterials[t];
}


void sxDefaultDestroyGrid(sxGrid g)
{free(g);}

void sxDefaultInitLoop(sxGrid g, sxSimplex sx)
{
  int i,j;

  sx->priv.grid = g;
  sx->priv.internal = NULL;
  sx->node = &(sx->priv.node_space[0]);
 
  
  sx->spaceDimension = g->spaceDim;
  switch(sx->iterator_type)
    {
    case sxNodes:	sx->dimension = 0; break;
    case sxEdges:	sx->dimension = 1; break;
    case sxCells:	sx->dimension = g->gridDim; break;
    case sxFaces:	sx->dimension = g->gridDim - 1; break;
    case sxBoundaryFaces:	sx->dimension = g->gridDim - 1; break;
    default: break;
    }
  sx->codimension = sx->spaceDimension-sx->dimension;


  switch(sx->dimension)
    {
    case 3:/* this is a cell */
      sx->numberOfNodes = 4;
      sx->numberOfEdges = 6;
      sx->numberOfFaces = 4;
      sx->numberOfEdgeNodes=2;
      sx->numberOfFaceNodes=3;

      sx->edgeNodes[1][1]=3;
      sx->edgeNodes[1][2]=4;
      sx->edgeNodes[1][3]=1;
      sx->edgeNodes[1][4]=2;


      sx->edgeNodes[2][1]=4;
      sx->edgeNodes[2][2]=2;
      sx->edgeNodes[2][3]=1;
      sx->edgeNodes[2][4]=3;


      sx->edgeNodes[3][1]=2;
      sx->edgeNodes[3][2]=3;
      sx->edgeNodes[3][3]=4;
      sx->edgeNodes[3][4]=1;

      sx->edgeNodes[4][1]=1;
      sx->edgeNodes[4][2]=2;
      sx->edgeNodes[4][3]=3;
      sx->edgeNodes[4][4]=4;

      sx->edgeNodes[5][1]=1;
      sx->edgeNodes[5][2]=3;
      sx->edgeNodes[5][3]=2;
      sx->edgeNodes[5][4]=4;

      sx->edgeNodes[6][1]=1;
      sx->edgeNodes[6][2]=4;
      sx->edgeNodes[6][3]=2;
      sx->edgeNodes[6][4]=3;

      sx->faceNodes[1][1]=2;
      sx->faceNodes[1][2]=3;
      sx->faceNodes[1][3]=4;
      sx->faceNodes[1][4]=1;


      sx->faceNodes[2][1]=1;
      sx->faceNodes[2][2]=3;
      sx->faceNodes[2][3]=4;
      sx->faceNodes[2][4]=2;


      sx->faceNodes[3][1]=1;
      sx->faceNodes[3][2]=2;
      sx->faceNodes[3][3]=4;
      sx->faceNodes[3][4]=3;

      sx->faceNodes[4][1]=1;
      sx->faceNodes[4][2]=2;
      sx->faceNodes[4][3]=3;
      sx->faceNodes[4][4]=4;

      break;
    case 2:
      switch (sx->spaceDimension)
	{
	case 2: /* this is cell */
	  sx->numberOfNodes = 3;
	  sx->numberOfEdges = 3;
	  sx->numberOfFaces = 3;
	  sx->numberOfFaceNodes = 2;
	  sx->numberOfEdgeNodes = 2;
	  
	  sx->edgeNodes[1][1]=2;
	  sx->edgeNodes[1][2]=3;
	  sx->edgeNodes[1][3]=1;
	  
	  sx->edgeNodes[2][1]=1;
	  sx->edgeNodes[2][2]=3;
	  sx->edgeNodes[2][3]=2;
	  
	  sx->edgeNodes[3][1]=1;
	  sx->edgeNodes[3][2]=2;
	  sx->edgeNodes[3][3]=3;
	  
	  sx->faceNodes[1][1]=2;
	  sx->faceNodes[1][2]=3;
	  sx->faceNodes[1][3]=1;
	  
	  sx->faceNodes[2][1]=1;
	  sx->faceNodes[2][2]=3;
	  sx->faceNodes[2][3]=2;
	  
	  sx->faceNodes[3][1]=1;
	  sx->faceNodes[3][2]=2;
	  sx->faceNodes[3][3]=3;
	  break;
	case 3:/* this is a boundary face */
	  sx->numberOfNodes = 3;
	  sx->numberOfEdges = 3;
	  sx->numberOfFaces = 1;
	  sx->numberOfFaceNodes = 3;
	  sx->numberOfEdgeNodes = 2;
	  
	  sx->edgeNodes[1][1]=2;
	  sx->edgeNodes[1][2]=3;
	  sx->edgeNodes[1][3]=1;
	  
	  sx->edgeNodes[2][1]=1;
	  sx->edgeNodes[2][2]=3;
	  sx->edgeNodes[2][3]=2;
	  
	  sx->edgeNodes[3][1]=1;
	  sx->edgeNodes[3][2]=2;
	  sx->edgeNodes[3][3]=3;
	  
	  sx->faceNodes[1][1]=1;
	  sx->faceNodes[1][2]=2;
	  sx->faceNodes[1][3]=3;
	  
	  break;
	}

      break;
    case 1:
      switch(sx->spaceDimension)
	{
	case 1:/* this is a cell */
	  sx->numberOfNodes = 2;
	  sx->numberOfEdges = 1;
	  sx->numberOfFaces = 2;
	  sx->numberOfFaceNodes = 1;
	  sx->numberOfEdgeNodes = 2;
	  
	  sx->edgeNodes[1][1]=1;
	  sx->edgeNodes[1][2]=2;

	  sx->faceNodes[1][1]=1;
	  sx->faceNodes[2][1]=2;
	  break;
	case 2:/* this is a boundary face */
	  sx->numberOfNodes = 2;
	  sx->numberOfEdges = 1;
	  sx->numberOfFaces = 1;
	  sx->numberOfFaceNodes = 2;
	  sx->numberOfEdgeNodes = 2;
	  
	  sx->edgeNodes[1][1]=1;
	  sx->edgeNodes[1][2]=2;
	  sx->faceNodes[1][1]=1;
	  sx->faceNodes[1][2]=2;
	  break;
	}
      break;
	case 0:
	  sx->numberOfNodes = 1;
	  sx->numberOfEdges = 0;
	  sx->numberOfFaces = 0;
	  sx->numberOfFaceNodes = 0;
	  sx->numberOfEdgeNodes = 0;
      break;
    }

  /* 
  switch (g->spaceDim)
    {
    case 1:
      sx->numberOfEdgeNodes=2;
      sx->edgeNodes[1][1]=1;
      sx->edgeNodes[1][2]=2;
      sx->numberOfFaceNodes=1;
      sx->faceNodes[1][1]=1;
      sx->faceNodes[2][1]=2;
      break;
    case 2:
      sx->numberOfEdgeNodes=3;
      sx->edgeNodes[1][1]=2;
      sx->edgeNodes[1][2]=3;
      sx->edgeNodes[1][3]=1;

      sx->edgeNodes[2][1]=1;
      sx->edgeNodes[2][2]=3;
      sx->edgeNodes[2][3]=2;

      sx->edgeNodes[3][1]=1;
      sx->edgeNodes[3][2]=2;
      sx->edgeNodes[3][3]=3;

      sx->numberOfFaceNodes=3;
      sx->faceNodes[1][1]=2;
      sx->faceNodes[1][2]=3;
      sx->faceNodes[1][3]=1;

      sx->faceNodes[2][1]=1;
      sx->faceNodes[2][2]=3;
      sx->faceNodes[2][3]=2;

      sx->faceNodes[3][1]=1;
      sx->faceNodes[3][2]=2;
      sx->faceNodes[3][3]=3;
      break;
    case 3:
      sx->numberOfEdgeNodes=4;
      sx->edgeNodes[1][1]=3;
      sx->edgeNodes[1][2]=4;
      sx->edgeNodes[1][3]=1;
      sx->edgeNodes[1][4]=2;


      sx->edgeNodes[2][1]=4;
      sx->edgeNodes[2][2]=2;
      sx->edgeNodes[2][3]=1;
      sx->edgeNodes[2][4]=3;


      sx->edgeNodes[3][1]=2;
      sx->edgeNodes[3][2]=3;
      sx->edgeNodes[3][3]=4;
      sx->edgeNodes[3][4]=1;

      sx->edgeNodes[4][1]=1;
      sx->edgeNodes[4][2]=2;
      sx->edgeNodes[4][3]=3;
      sx->edgeNodes[4][4]=4;

      sx->edgeNodes[5][1]=1;
      sx->edgeNodes[5][2]=3;
      sx->edgeNodes[5][3]=2;
      sx->edgeNodes[5][4]=4;

      sx->edgeNodes[6][1]=1;
      sx->edgeNodes[6][2]=4;
      sx->edgeNodes[6][3]=2;
      sx->edgeNodes[6][4]=3;

      sx->numberOfFaceNodes=4;
      sx->faceNodes[1][1]=2;
      sx->faceNodes[1][2]=3;
      sx->faceNodes[1][3]=4;
      sx->faceNodes[1][4]=1;


      sx->faceNodes[2][1]=1;
      sx->faceNodes[2][2]=3;
      sx->faceNodes[2][3]=4;
      sx->faceNodes[2][4]=2;


      sx->faceNodes[3][1]=1;
      sx->faceNodes[3][2]=2;
      sx->faceNodes[3][3]=4;
      sx->faceNodes[3][4]=3;

      sx->faceNodes[4][1]=1;
      sx->faceNodes[4][2]=2;
      sx->faceNodes[4][3]=3;
      sx->faceNodes[4][4]=4;

      break;
    }
    */

  for(i=1;i<=sx->numberOfNodes;i++)
    sx->coordinates[i]=sx->priv.coord_space+(i-1)*sx->spaceDimension;
  sx->priv.coord_nodes=sx->priv.coord_node_space;
  
  sx->priv.coord_ptr=sx->priv.coord_space;
  sx->priv.coord_nodes[1]=1;
  sx->priv.coord_nodes[2]=2;
  sx->priv.coord_nodes[3]=3;
  sx->priv.coord_nodes[4]=4;

  if (g->setAdditionalData!=NULL)
    g->setAdditionalData(g,sx);

}