wzsxgrid.cxx

Delaunay三角形的网格剖分程序

//static char *rcsid="$Id: wzsxgrid.cxx,v 1.2 2000/07/28 06:37:05 schmelze Exp $";
//
// $Log: wzsxgrid.cxx,v $
// Revision 1.2  2000/07/28 06:37:05  schmelze
// COG 2.1 alpha
//
// Revision 1.1  1997/11/05  19:14:36  schmelze
// Version mit neuer Dokumentation, erste Version des neuen Gittergenerators
// noch ohne Randlinienbehandlung
//
// Revision 1.5  1997/10/15  16:58:07  fuhrmann
// face access for interfaces in 1d/2d corrected
//
// Revision 1.4  1997/10/10  14:16:53  fuhrmann
// femmatrices in kaskade, ifm stuff checkin
//
// Revision 1.3  1997/10/08  12:31:30  fuhrmann
// 1d,2d face cycles, verbose test
//
// Revision 1.2  1997/10/01  18:27:32  fuhrmann
// Merge in of kaskade-3.2
//
// Revision 1.1.1.1  1997/09/26  18:36:17  darcy
// preliminary checkin
//
// Revision 1.3  1995/11/10  20:59:24  fuhrmann
// Refinement/Interpolation interface built in.
//
//

extern "C"
{
#include "stdsys/dprintf.h"
#include "stdsys/memory.h"
#include "sxgrid/sxgrid.h"
#include "sxgrid/celltype.h"


// These should come from celltype

  extern int cellInit();
  extern void FVM1Edge(double*,int*,double*,double *,double*),
  FVM2Triangle(double*,int*,double*,double *,double*),
  FVM3Tetrahedron(double*,int*,double*,double *,double*);
  extern void FEM1Edge(double*,int*,double*,double *,double*,double*),
  FEM2Triangle(double*,int*,double*,double *,double*,double*),
  FEM3Tetrahedron(double*,int*,double*,double *,double*,double*);
  extern void FVM0Vertex(double*,int*,double*,double *,double*),
  FVM2Edge(double*,int*,double*,double *,double*),
  FVM3Triangle(double*,int*,double*,double *,double*);
  extern void FEM0Vertex(double*,int*,double*,double *,double*,double*),
  FEM2Edge(double*,int*,double*,double *,double*,double*),
  FEM3Triangle(double*,int*,double*,double *,double*,double*);
  
}

#include <iostream.h>
#include "general.h"
#include "utils.h"
#include "cmdpars.h"
#include "dimension.h"
#include "triang.h"
#include "triang1.h"
#include "triang2.h"
#include "triang3.h"

#include "vector.h"
#ifndef KASKADE_LIB
#include "feplot.h"
#endif


//  !aufteilen in Basis & Ableitung
struct PrivateSimplexDataStruct
{
  void (*FVMFactors)(double*,int*,double*,double*,double*);
  void (*FEMFactors)(double*,int*,double*,double*,double*,double*);
  void (*FVEMFactors)(double*,int*,double*,double*,double*,double*);

  void (*sxNext)(Simplex sx);

  EDG1 *cell1;
  TR2 *cell2;
  TET3 *cell3;
  DListIter<EDG1> *celliter1;
  DListIter<TR2> *celliter2;
  DListIter<TET3> *celliter3;
  

  PT1 *face1;
  EDG2 *face2;
  TR3  *face3;
  DListIter<PT1> *faceiter1;
  DListIter<EDG2> *faceiter2;
  DListIter<TR3> *faceiter3;
  
  PT1 *node1;
  PT2 *node2;
  PT3  *node3;
  DListIter<PT1> *nodeiter1;
  DListIter<PT2> *nodeiter2;
  DListIter<PT3> *nodeiter3;

  EDG1 *line1;
  EDG2 *line2;
  EDG3 *line3;
  DListIter<EDG1> *lineiter1;
  DListIter<EDG2> *lineiter2;
  DListIter<EDG3> *lineiter3;


  PATCH* patch;
  double coord[15];

  double stiffmat[30];
  double massmat[30];
};


struct sxGridStruct
{
  MESH *mesh;
  int spacedim;
  
  int node_num;
  int cell_num;
  int face_num;
  int all_face_num;
  int line_num;
  int cell_node_num;
  int face_node_num;
  
  int  cell_line_num;
  int  cell_face_num;
  int  cell_line_node_ptr[7][7];
  int  cell_face_node_ptr[7][7];
  
  
  int *node_mat_node_num;
  int *cell_mat_cell_num;
  int *face_mat_face_num;
  
  int node_mat_num;
  int cell_mat_num;
  int face_mat_num;
  
  
  double xmin,xmax;
  double ymin,ymax;
  double zmin,zmax;

  int sx_node[5];
  int sx_line[8];
  int sx_face[5];
  
  int prev_node_num; // for interpolation
  int bisect;
};


typedef struct sxGridStruct sxGridData;


// debug flags, siehe dprintf.h
static int dbg=1;
static int inf=1;
static int wrn=1;
static int err=1;

static int initflag=0;


static void init(void)
{
  if (initflag) return;

  cout.sync_with_stdio();

  dbgRegister("sxgrid-dbg",&dbg);
  dbgRegister("sxgrid",&inf);
  dbgRegister("sxgrid-err",&err);
  dbgRegister("sxgrid-wrn",&wrn);

  cellInit();

  initflag=1;
}



//////////////////////////////////
int node_node_ptr[2]={ 0 ,1 };


//////////////////////////////////
// CELL Loops
//////////////////////////////////


inline void setCell1Data(Simplex sx, EDG1 *cell1)
{
  sx->priv->patch=cell1;
  sx->material=cell1->Class();
  sx->index=cell1->getNode();
  sx->coord[1][1]=cell1->p1->x;
  sx->coord[2][1]=cell1->p2->x;
  sx->node[1]=cell1->p1->getNode();
  sx->node[2]=cell1->p2->getNode();
  sx->line[1]=sx->index;
}

void sxNextCell1(Simplex sx)
{
  sx->priv->cell1=sx->priv->celliter1->all();
  if (!sx->priv->cell1)
    {
      delete sx->priv->celliter1;
      delete sx->priv;
      sx->priv=NULL;
      return;
    }
  setCell1Data(sx,sx->priv->cell1);
} 

void setCell2Data(Simplex sx, TR2 *cell2)
{
  sx->priv->patch=cell2;
  sx->material=cell2->Class();
  sx->index=cell2->getNode();
  sx->coord[1][1]=cell2->p1->x;
  sx->coord[1][2]=cell2->p1->y;
  sx->coord[2][1]=cell2->p2->x;
  sx->coord[2][2]=cell2->p2->y;
  sx->coord[3][1]=cell2->p3->x;
  sx->coord[3][2]=cell2->p3->y;
  sx->node[1]=cell2->p1->getNode();
  sx->node[2]=cell2->p2->getNode();
  sx->node[3]=cell2->p3->getNode();

  sx->line[1]=cell2->e1->getNode();
  sx->line[2]=cell2->e2->getNode();
  sx->line[3]=cell2->e3->getNode();

}


void sxNextCell2(Simplex sx)
{
  sx->priv->cell2=sx->priv->celliter2->all();
  if (!sx->priv->cell2)
    {
      delete sx->priv->celliter2;
      delete sx->priv;
      sx->priv=0;
      return;
    }
  setCell2Data(sx,sx->priv->cell2);

}

inline void setCell3Data(Simplex sx, TET3 *cell3)
{
  sx->priv->patch=cell3;
  sx->material=cell3->Class();
  sx->index=cell3->getNode();
  sx->coord[1][1]=cell3->p1->x;
  sx->coord[1][2]=cell3->p1->y;
  sx->coord[1][3]=cell3->p1->z;
  sx->coord[2][1]=cell3->p2->x;
  sx->coord[2][2]=cell3->p2->y;
  sx->coord[2][3]=cell3->p2->z;
  sx->coord[3][1]=cell3->p4->x;
  sx->coord[3][2]=cell3->p4->y;
  sx->coord[3][3]=cell3->p4->z;
  sx->coord[4][1]=cell3->p3->x;
  sx->coord[4][2]=cell3->p3->y;
  sx->coord[4][3]=cell3->p3->z;
// Orientierung !!! ??? 

  sx->node[1]=cell3->p1->getNode();
  sx->node[2]=cell3->p2->getNode();
  sx->node[3]=cell3->p4->getNode();
  sx->node[4]=cell3->p3->getNode();


  sx->line[1]=cell3->e1->getNode();
  sx->line[2]=cell3->e2->getNode();
  sx->line[3]=cell3->e3->getNode();
  sx->line[4]=cell3->e4->getNode();
  sx->line[5]=cell3->e5->getNode();
  sx->line[6]=cell3->e6->getNode();

  sx->face[1]=cell3->n1->getNode();
  sx->face[2]=cell3->n2->getNode();
  sx->face[3]=cell3->n3->getNode();
  sx->face[4]=cell3->n4->getNode();

}


void sxNextCell3(Simplex sx)
{
  sx->priv->cell3=sx->priv->celliter3->all();
  if (!sx->priv->cell3)
    {
      delete sx->priv->celliter3;
      delete sx->priv;
      sx->priv=0;
      return;
    }

  setCell3Data(sx,sx->priv->cell3);
}

//////////////////////////////////
// FACE Loops
//////////////////////////////////




void sxNextFace1(Simplex sx)
{
  while ((sx->priv->face1=sx->priv->faceiter1->all())&&!(sx->priv->face1->isDirichlet()));
  if (!sx->priv->face1)
    {
      delete sx->priv->faceiter1;
      delete sx->priv;
      sx->priv=0;
      return;
    }
  sx->priv->patch=sx->priv->face1;
  sx->material=sx->priv->face1->Class();
  sx->index=sx->priv->face1->getNode();
  sx->coord[1][1]=sx->priv->face1->x;
  sx->node[1]=sx->priv->face1->getNode();
}


void sxNextFace2(Simplex sx)
{
  while ((sx->priv->face2=sx->priv->faceiter2->all())&&!sx->priv->face2->isDirichlet());
  if (!sx->priv->face2)
    {
      delete sx->priv->faceiter2;
      delete sx->priv;
      sx->priv=0;
      return;
    }

  sx->priv->patch=sx->priv->face2;
  sx->material=sx->priv->face2->Class();
  sx->index=sx->priv->face2->getNode();
  sx->coord[1][1]=sx->priv->face2->p1->x;
  sx->coord[1][2]=sx->priv->face2->p1->y;
  sx->coord[2][1]=sx->priv->face2->p2->x;
  sx->coord[2][2]=sx->priv->face2->p2->y;
  sx->node[1]=sx->priv->face2->p1->getNode();
  sx->node[2]=sx->priv->face2->p2->getNode();
  sx->line[1]=sx->index;
}

void sxNextFace3(Simplex sx)
{
  while ((sx->priv->face3=sx->priv->faceiter3->all())&&!sx->priv->face3->isDirichlet());
  if (!sx->priv->face3)
    {
      delete sx->priv->faceiter3;
      delete sx->priv;
      sx->priv=0;
      return;
    }
  sx->priv->patch=sx->priv->face3;
  sx->material=sx->priv->face3->Class();
  sx->index=sx->priv->face3->getNode();
  sx->coord[1][1]=sx->priv->face3->p1->x;
  sx->coord[1][2]=sx->priv->face3->p1->y;
  sx->coord[1][3]=sx->priv->face3->p1->z;
  sx->coord[2][1]=sx->priv->face3->p2->x;
  sx->coord[2][2]=sx->priv->face3->p2->y;
  sx->coord[2][3]=sx->priv->face3->p2->z;
  sx->coord[3][1]=sx->priv->face3->p3->x;
  sx->coord[3][2]=sx->priv->face3->p3->y;
  sx->coord[3][3]=sx->priv->face3->p3->z;
 
  sx->node[1]=sx->priv->face3->p1->getNode();
  sx->node[2]=sx->priv->face3->p2->getNode();
  sx->node[3]=sx->priv->face3->p3->getNode();

  sx->line[1]=sx->priv->face3->e1->getNode();
  sx->line[2]=sx->priv->face3->e2->getNode();
  sx->line[3]=sx->priv->face3->e3->getNode();

}

//////////////////////////////////
// EDGE Loops
//////////////////////////////////
void sxNextLine1(Simplex sx)
{
  sx->priv->line1=sx->priv->lineiter1->all();
  if (!sx->priv->line1)
    {
      delete sx->priv->lineiter1;
      delete sx->priv;
      sx->priv=0;
      return;
    }
  sx->priv->patch=sx->priv->line1;
  sx->index=sx->priv->line1->getNode();
  sx->node[1]=sx->priv->line1->p1->getNode();
  sx->node[2]=sx->priv->line1->p2->getNode();
}


void sxNextLine2(Simplex sx)
{
  sx->priv->line2=sx->priv->lineiter2->all();
  if (!sx->priv->line2)
    {
      delete sx->priv->lineiter2;
      delete sx->priv;
      sx->priv=0;
      return;
    }
  
  sx->priv->patch=sx->priv->line2;
  sx->index=sx->priv->line2->getNode();
  sx->node[1]=sx->priv->line2->p1->getNode();
  sx->node[2]=sx->priv->line2->p2->getNode();
}

void sxNextLine3(Simplex sx)
{
  sx->priv->line3=sx->priv->lineiter3->all();
  if (!sx->priv->line3)
    {
      delete sx->priv->lineiter3;
      delete sx->priv;
      sx->priv=0;
      return;
    }
  sx->priv->patch=sx->priv->line3;
  sx->index=sx->priv->line3->getNode();
  sx->node[1]=sx->priv->line3->p1->getNode();
  sx->node[2]=sx->priv->line3->p2->getNode();

}


//////////////////////////////////
// NODE Loops
//////////////////////////////////



void sxNextNode1(Simplex sx)
{
  sx->priv->node1=sx->priv->nodeiter1->all();
  if (!sx->priv->node1)
    {
      delete sx->priv->nodeiter1;
      delete sx->priv;
      sx->priv=0;
      return;
    }
  sx->priv->patch=sx->priv->node1;
  sx->material=sx->priv->node1->Class();
  sx->index=sx->priv->node1->getNode();
  sx->coord[1][1]=sx->priv->node1->x;
  sx->node[1]=sx->priv->node1->getNode();
}


void sxNextNode2(Simplex sx)
{
  sx->priv->node2=sx->priv->nodeiter2->all();
  if (!sx->priv->node2)
    {
      delete sx->priv->nodeiter2;
      delete sx->priv;
      sx->priv=0;
      return;
    }
  
  sx->priv->patch=sx->priv->node2;
  sx->material=sx->priv->node2->Class();
  sx->index=sx->priv->node2->getNode();
  sx->coord[1][1]=sx->priv->node2->x;
  sx->coord[1][2]=sx->priv->node2->y;
  sx->node[1]=sx->priv->node2->getNode();
}

void sxNextNode3(Simplex sx)
{
  sx->priv->node3=sx->priv->nodeiter3->all();
  if (!sx->priv->node3)
    {
      delete sx->priv->nodeiter3;
      delete sx->priv;
      sx->priv=0;
      return;
    }
  sx->priv->patch=sx->priv->node3;
  sx->material=sx->priv->node3->Class();
  sx->index=sx->priv->node3->getNode();
  sx->coord[1][1]=sx->priv->node3->x;
  sx->coord[1][2]=sx->priv->node3->y;
  sx->coord[1][3]=sx->priv->node3->z;
 
  sx->node[1]=sx->priv->node3->getNode();
}

//////////////////////////////////



void sxNext(Simplex sx)
{
  sx->priv->sxNext(sx);
}


void sxStart(sxGrid g, Simplex sx, int WhatToTraverse)
{
  register int in;
  sx->grid=g;
  sx->spacedim=g->spacedim;
  sx->priv=new PrivateSimplexData;
  for(in=0;in<15;in++) sx->priv->coord[in]=3.3333333e33;
  sx->node=sx->grid->sx_node;
  sx->line=sx->grid->sx_line;
  sx->face=sx->grid->sx_face;
  
  switch (WhatToTraverse)
    {
    case SX_CELLS:
      sx->nnodes=g->cell_node_num;
      sx->nlines=g->cell_line_num;
      sx->nfaces=g->cell_face_num;
      sx->linenodenum=2;
      sx->facenodenum=g->face_node_num;


      sx->linenodes=g->cell_line_node_ptr;
      sx->facenodes=g->cell_face_node_ptr;
      sx->dim=g->spacedim;
      sx->codim=0;

      switch(sx->spacedim)
	{
	case 1: 
	  sx->face=sx->node;
	  sx->priv->FVMFactors=FVM1Edge; 
	  sx->priv->FEMFactors=FEM1Edge; 
          sx->priv->sxNext=sxNextCell1;
	  sx->priv->celliter1= new DListIter<EDG1>(sx->grid->mesh->castToMESH1()->lines());
	  break;
	case 2: 
	  sx->face=sx->line;
	  sx->priv->FVMFactors=FVM2Triangle; 
	  sx->priv->FEMFactors=FEM2Triangle; 
          sx->priv->sxNext=sxNextCell2;
	  sx->priv->celliter2= new DListIter<TR2>(sx->grid->mesh->castToMESH2()->triangles());
	  break;
  	case 3: 
	  sx->priv->FVMFactors=FVM3Tetrahedron; 
	  sx->priv->FEMFactors=FEM3Tetrahedron; 
          sx->priv->sxNext=sxNextCell3;
	  sx->priv->celliter3= new DListIter<TET3>(sx->grid->mesh->castToMESH3()->tetras());
	  break;
	}

      break;

    case SX_FACES:
      sx->nnodes=g->face_node_num;
      sx->nfaces=1;
      sx->linenodes=NULL;

      switch(sx->spacedim)
	{
	case 1: 
	  sx->face=sx->node;
	  sx->nlines=0;
	  sx->priv->FVMFactors=FVM0Vertex; 
	  sx->priv->FEMFactors=FEM0Vertex; 
          sx->priv->sxNext=sxNextFace1;
	  sx->priv->faceiter1= new DListIter<PT1>(sx->grid->mesh->castToMESH1()->points());
	  break;
	case 2: 
	  sx->face=sx->line;