wzsxgrid.cxx

Delaunay三角形的网格剖分程序

	PT2 *patch;
	DListIter<PT2> iter(m->castToMESH2()->points());
	while ((patch = iter.all())) 
	  {
	    if (patch->x< g->xmin) g->xmin=patch->x;
	    if (patch->x> g->xmax) g->xmax=patch->x;
	    if (patch->y< g->ymin) g->ymin=patch->y;
	    if (patch->y> g->ymax) g->ymax=patch->y;
	  }
	g->zmin=0.0;
	g->zmax=0.0;
      }
      break;
    case 3:
      {
	PT3 *patch;
	DListIter<PT3> iter(m->castToMESH3()->points());
	while ((patch =iter.all())) 
	  {
	    if (patch->x< g->xmin) g->xmin=patch->x;
	    if (patch->x> g->xmax) g->xmax=patch->x;
	    if (patch->y< g->ymin) g->ymin=patch->y;
	    if (patch->y> g->ymax) g->ymax=patch->y;
	    if (patch->z< g->zmin) g->zmin=patch->z;
	    if (patch->z> g->zmax) g->zmax=patch->z;
	  }
      }
      break;
      
    }
  
  ///////////////////////////////////////
  ///// parse cell materials

  PATCH *patch;
  i=0;
  g->cell_mat_num=0;
  m->resetElemIter();
  while ((patch = m->elemIterAll())) 
    { 
       if (patch->Class()>g->cell_mat_num) g->cell_mat_num=patch->Class();
    }
   
  g->cell_mat_cell_num=(int*)calloc(g->cell_mat_num+1,sizeof(int));
  
   /////////////////////////////////////////////////////
  ///// parse face materials (boundary condition types)
  

  g->face_node_num=g->spacedim;
  g->face_mat_num=0;
  
  switch(spacedim)
    { 
    case 1:
      {
      PT1* patch;
      DListIter<PT1>  iter(m->castToMESH1()->points());
      while ((patch = iter.all()))     
	if (patch->onBoundary()&&patch->Class())
	  {
	    if (patch->Class()>g->face_mat_num) 
	      g->face_mat_num=patch->Class();
	  }
      }
      break;
    case 2:
      {
      EDG2* patch;
      DListIter<EDG2>  iter(m->castToMESH2()->lines());
      i=0;
      while ((patch = iter.all()))
	if (patch->onBoundary()&&patch->Class())
	  {
	    if (patch->Class()>g->face_mat_num) 
	      g->face_mat_num=patch->Class();
	  }
      }
      break;
    case 3:
      {
      TR3* patch;
      DListIter<TR3>  iter(m->castToMESH3()->triangles());;
      while ((patch = iter.all()))
	if (patch->onBoundary()&&patch->Class())
	  {
	    if (patch->Class()>g->face_mat_num) 
	      g->face_mat_num=patch->Class();
	  }
      break;
      }
    }
  g->face_mat_face_num=(int*) calloc(g->face_mat_num+1,sizeof(int));
  
  number_stuff(g);
  
  sxgPrintData(g);  
  dprintf(inf)("sxGrid created from  mesh\n",meshname);
  return g;
}

void sxDestroy(sxGrid g)
{
  dprintf(inf)("Deleting sxgrid...\n");
  free(g->cell_mat_cell_num);
  free(g->face_mat_face_num);
  delete g->mesh;
  delete g;
  dprintf(inf)("Deletion done.\n");
  
}


/////////////////////////////////////////
// Interface to mesh refinement
/////////////////////////////////////////


// methods on simplices

int sxDepth(Simplex sx)
{
// My depth is relative to the finest level!!! 
  return sx->grid->mesh->maxDepth-sx->priv->patch->Depth();
}

void sxMark(Simplex sx)
{
 sx->priv->patch->setMark();
}

void sxUnMark(Simplex sx)
{
 sx->priv->patch->unMark();
}

int sxMarked(Simplex sx)
{
  if (sx->priv->patch->marked()) return 1; else return 0;
}


//////// methods on mesh
int sxMaxDepth(sxGrid g)
{
  return g->mesh->maxDepth;
}

void sxRefine(sxGrid g)
{
  MESH *m=g->mesh;
  int i;
  dprintf(inf)("Refining sxGrid...\n");
  m->Refine();
  g->prev_node_num=g->node_num;
  
  number_stuff(g);
  if (debug(inf)) sxgPrintData(g);
  dprintf(inf)("Refinement done\n");
}


void sxFlatten(sxGrid g)
{
  g->mesh->Flatten();
}

//// Interpolation of vectors 


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

static double gshape[30];
static int sx_nodes[5]={1,2,3,4};

#define COORD(DIM,NODE,DIRm1) coord[1+DIRm1+(cellnodes[NODE]-1)*DIM]

void sxFVMFactors(Simplex sx)
{
 sx->priv->FVMFactors(
		      sx->priv->coord,
		      sx_nodes,
		      &sx->volume,
		      sx->nodefac+1,
		      sx->linefac+1
		      );
}

void sxFVMMatrices(Simplex sx)
{
  register   int i,j;
  sx->priv->FVMFactors(
		       sx->priv->coord,
		       sx_nodes,
		       &sx->volume,
		       sx->nodefac+1,
		       sx->linefac+1
		       );
 for(i=1;i<=sx->nnodes;i++)
   for(j=1;j<=sx->nnodes;j++)
     sx->stiffmat[i][j]=sx->massmat[i][j]=0.0;

  for(i=1;i<=sx->nlines;i++) 
    {
      sx->stiffmat[sx->linenodes[i][1]][sx->linenodes[i][2]]=-sx->linefac[i];
      sx->stiffmat[sx->linenodes[i][2]][sx->linenodes[i][1]]=-sx->linefac[i];
      sx->stiffmat[sx->linenodes[i][1]][sx->linenodes[i][1]]+=sx->linefac[i];
      sx->stiffmat[sx->linenodes[i][2]][sx->linenodes[i][2]]+=sx->linefac[i];
    }
  for(i=1;i<=sx->nnodes;i++)
    sx->massmat[i][i]=sx->nodefac[i];
 
}


void sxFEMFactors(Simplex sx)
{
  register int id,in,ie;
  sx->priv->FEMFactors(
		       sx->priv->coord,
		       sx_nodes,
		       &sx->volume,
		       sx->nodefac+1,
		       sx->linefac+1,
		       gshape
		       );
  for(ie=1;ie<=sx->nlines;ie++) 
    sx->linefac[ie]= -sx->linefac[ie];
  for(in=1;in<=sx->nnodes;in++) 
    sx->nodefac[in]=sx->volume/((double)sx->nnodes);


  for(in=1;in<=sx->nnodes;in++)
    for(id=1;id<=sx->spacedim;id++)
	sx->dgrad[in][id]=gshape[(in-1)*sx->spacedim+(id-1)]; 

}

static double xn[4]=
{
  0,
  0.1666666666,
  0.0833333333,
  0.05
};


void sxFEMMatrices(Simplex sx)
{
  register   int i,j,in,id;
  double xv;
  sx->priv->FVMFactors(
		       sx->priv->coord,
		       sx_nodes,
		       &sx->volume,
		       sx->nodefac+1,
		       sx->linefac+1
		       );

   for(i=1;i<=sx->nnodes;i++)
     sx->stiffmat[i][i]=0.0;
 
  for(i=1;i<=sx->nlines;i++) 
    {
      sx->stiffmat[sx->linenodes[i][1]][sx->linenodes[i][2]]=-sx->linefac[i];
      sx->stiffmat[sx->linenodes[i][2]][sx->linenodes[i][1]]=-sx->linefac[i];
      sx->stiffmat[sx->linenodes[i][1]][sx->linenodes[i][1]]+=sx->linefac[i];
      sx->stiffmat[sx->linenodes[i][2]][sx->linenodes[i][2]]+=sx->linefac[i];
    }
  /* npar is wrong in celltye!!!*/
  
  
  xv=xn[sx->dim]*sx->volume;
  for(i=1;i<=sx->nnodes;i++)
    for(j=1;j<=sx->nnodes;j++)
      sx->massmat[i][j]=xv;
  
  for(i=1;i<=sx->nnodes;i++)
    sx->massmat[i][i]+=xv;
 
  for(in=1;in<=sx->nnodes;in++)
    for(id=1;id<=sx->dim;id++)
	sx->dgrad[in][id]=gshape[(in-1)*sx->spacedim+(id-1)]; 
}


void sxFVEMFactors(Simplex sx)
{
  register int id,in;
  
  sx->priv->FEMFactors(
		       sx->priv->coord,
		       sx_nodes,
		       &sx->volume,
		       sx->nodefac+1,
		       sx->linefac+1,
		       gshape
		       );


  for(in=1;in<=sx->nnodes;in++)
    for(id=1;id<=sx->spacedim;id++)
	sx->dgrad[in][id]=gshape[(in-1)*sx->spacedim+(id-1)]; 

  sx->priv->FVMFactors(
		       sx->priv->coord,
		       sx_nodes,
		       &sx->volume,
		       sx->nodefac+1,
		       sx->linefac+1
		       );
}




double* sxCenter(Simplex sx, double *ccrd)
{
  register int i,idim;
  for (idim=1;idim<=sx->spacedim;idim++) ccrd[idim]=0.0;
  for(i=1;i<=sx->nnodes;i++)
    {
      for (idim=1;idim<=sx->spacedim;idim++)
	ccrd[idim]+=sx->coord[i][idim];
    }
  
  for (idim=1;idim<=sx->spacedim;idim++)
    ccrd[idim]/=(double)(sx->nnodes);
  return ccrd;
}





////////////////////////////////////////////////////////////////
sxFunction sxCreateFunction(sxGrid g, char *name,int num, int comp)
{
  double *d=NULL;
  switch (comp)
    {
    case SX_NODES: 
      d=new double[num*g->node_num+1]; return d;
    case SX_CELLS:
      d=new double[num*g->cell_num+1]; return d;
    case SX_EDGES:
      d=new double[num*g->line_num+1]; return d;
      return NULL;
    case SX_ALL_FACES: 
      d=new double[num*g->all_face_num+1]; return d;
      return NULL;
    default:
      dprintf(err)("Function type %d not supported",comp);
    }
  return NULL;
}

void sxDestroyFunction(sxGrid g, sxFunction f)
{
  delete f;
}

// this interpolates _only_ from old level to new level
sxFunction sxInterpolateFunctionFromPrevLevel
(sxGrid g, sxFunction f_old, int num, int comp)
{
  double *f_new=NULL;
  MESH *mesh=g->mesh;

  if (num>1) 
    {
      dprintf(err)("sorry, interpolation for multiple functions currently not implemented\n");
      return NULL;
    }
  

  switch (comp)
    {
    case SX_NODES: 
      int i;

      // inject coarse mesh values
      // assume old values come first!? 
      f_new=new double[num*g->node_num+1];
      for (i=1;i<=g->prev_node_num;i++) f_new[i]=f_old[i];
      for (;i<=g->node_num;i++) f_new[i]=0.0;
// MLNode seems to be  undefined unless
// LinElemInt:: setMLNodeNumbers has been called!
//      {
//	int i, node_old,node_new;
//	PATCH* p;
//	f_new=new double[num*g->node_num+1];
//	for (i=1;i<=g->node_num;i++) f_new[i]=0.0;
//	
//	mesh->resetPtIter();
//	while (p=mesh->ptIterAll()) 
//	  if (p->Depth()<mesh->maxDepth) 
//	    {
//	      node_new  = p->getNode();
//	      node_old = p->MLNode(1);
//	      f_new[node_new] = f_old[node_old]; 
//	    }
//     }

      //interplolate new values
      {
	
	
	PATCH* ed;
        PATCH *pm;
	int    node1, node2, pmNode;
	int iint=0;
	mesh->resetEdgIter();
	while ((ed=mesh->edgIterAllOfHistory())) 
	  {
	    if (ed->refined())	
	      {
		pm = ed->getMidPoint();
		
		pmNode = pm->getNode();
		if (pmNode>g->prev_node_num)
		  {
                    iint++;
		    const Vector<PT*>& p=ed->getPoints();
		    node1 = p[1]->getNode();
		    node2 = p[2]->getNode();
		    f_new[pmNode] = 0.5*(f_new[node1] + f_new[node2]);
		  }
	      }
	  }

      }
      
      delete f_old;
      return f_new;
      break;
    case SX_CELLS:
      dprintf(err)("sorry, interpolation of cell functions currently not implemented\n");
      return NULL;
      break;
    case SX_EDGES:
      dprintf(err)("sorry, interpolation of line functions currently not implemented\n");
      return NULL;
      break;
    case SX_ALL_FACES: 
    case SX_FACES: 
      dprintf(err)("sorry, interpolation of face functions currently not implemented\n");
      return NULL;
    }
  return NULL;
}

int sxAdaptive(sxGrid g)
{
  return 1;
}

//get data of cell with x in it. 
//void sxFindCell(SimplexGrid g,
//		double coord[],
//		double baryc_coord[],
//		SimplexData *sx)
//{
//
//  p=findPatch(coord,baryc_coord);
//  SetCellData(sx,p);
//  baryc_coord[2]=1.0-baryc_coord[1]-baryc_coord[0];
//  return;
//
//}


////////////////////////////////////////////////////////////////
enum pauseMode  Pause()
{
  static Bool Continue = False;
  
  if (Continue == True) return noPicture;
  if (!Cmd.isTrue("pause")) { Continue = True; return noPicture; }
  
  char s[5]; cout << " <CR>"; cout.flush(); gets(s); strToLower(s);
  if (strchr(s,'q') || strchr(s,'e')) { cout << "\nEXIT FORCED\n"; exit(1); }
  if (strchr(s,'c') || strchr(s,'g')) Continue = True;
  return noPicture;
}


void sxLocatePoint(sxGrid g, double *coord, Simplex sx, double *bcoord)
{
  Vector<Real>x(4); Vector<Real>xUnit(4);
  
  for(int i=1;i<=g->spacedim;i++)
    x[i]=coord[i];
  
  switch (g->spacedim)
    {
    case 1:
      EDG1 * cell1;
      cell1 = (EDG1*) g->mesh->findPatch(x,xUnit,cell1);
      if (cell1==0) {sx->index=0; return;}
      setCell1Data(sx,cell1);
      bcoord[1]=1.0-xUnit[1];
      bcoord[2]=xUnit[1];
      break;
    case 2:
      TR2 * cell2;
      cell2=(TR2*)g->mesh->findPatch(x,xUnit,cell2);
      if (cell2==0) {sx->index=0; return;}
      // this order is somehow hard coded, but how
      // at least we can look it up in elements2.cc!

      bcoord[1]=1.0-xUnit[1]-xUnit[2];
      bcoord[2]=xUnit[1];
      bcoord[3]=xUnit[2];
      setCell2Data(sx,cell2);
      break;
   case 3:
      TET3 * cell3;
      cell3=(TET3*)g->mesh->findPatch(x,xUnit,cell3);
      if (cell3==0) {sx->index=0; return;}
      setCell3Data(sx,cell3);
      bcoord[1]=1.0-xUnit[1]-xUnit[2]-xUnit[3];
      bcoord[2]=xUnit[1];
      bcoord[3]=xUnit[2];
      bcoord[4]=xUnit[3];
      break;

    }
}

sxFunction sxInterpolateFunctionFromPrevGrid
(sxGrid g, sxGrid g_old, sxFunction f_old, int num, int comp)
{
  double *f_new=NULL;
  MESH *mesh=g->mesh;

  if (num>1) 
    {
      dprintf(err)("sorry, interpolation for multiple functions currently not implemented\n");
      return NULL;
    }
  

  switch (comp)
    {
    case SX_NODES: 
      int i;
      SimplexData sx,sx_old;
      sx_old.priv=new PrivateSimplexData;
      for(i=1;i<=g->spacedim+1;i++)
	sx_old.coord[i]=sx_old.priv->coord+(i-1)*g->spacedim;
      sx_old.node= new int[5];
      sx_old.line= new int[5];
      
      double bcoord[5];
      f_new=new double[num*g->node_num+1];
      
      SX_LOOP(g,sx,SX_NODES)
	{
	  sxLocatePoint(g_old,sx.coord[1],&sx_old,bcoord);
	  double f=0.0;
	  for (int in=1;in<=g->spacedim+1;in++)
	    //  	    cout<< bcoord[in]<< ' ' << f_old[sx_old.node[in]]<<endl<<flush,
	    f+=bcoord[in]*f_old[sx_old.node[in]];
	  // cout<<endl;
	  f_new[sx.index]=f;
	}
      delete[] sx_old.node;
      delete[] sx_old.line;
      delete sx_old.priv;
      return f_new;
      break;
    case SX_CELLS:
      dprintf(err)("sorry, interpolation of cell functions currently not implemented\n");
      return NULL;
      break;
    case SX_EDGES:
      dprintf(err)("sorry, interpolation of line functions currently not implemented\n");
      return NULL;
      break;
    case SX_ALL_FACES: 
    case SX_FACES: 
      dprintf(err)("sorry, interpolation of face functions currently not implemented\n");
      return NULL;
    }
  return NULL;
}


sxGrid global_sxg=NULL;