mesh_smoother_vsmoother.c

一个用来实现偏微分方程中网格的计算库

// $Id: mesh_smoother_vsmoother.C 2858 2008-06-13 18:59:07Z benkirk $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#include "libmesh_config.h"
#ifdef ENABLE_VSMOOTHER

// C++ includes
#include <algorithm> // for std::copy, std::sort


// Local includes
#include "mesh_smoother_vsmoother.h"
#include "mesh_tools.h"
#include "elem.h"
#include "unstructured_mesh.h"
#include "utility.h"

// Member functions for the Variational Smoother
double VariationalMeshSmoother::smooth(unsigned int)
{
  int n, me, gr, adp, miniter, miniterBC, maxiter, N, ncells, nedges, s, err=0;
  double theta;
  char grid[40], metr[40], grid_old[40], adap[40];
//  FILE *stream;
  LPLPLPDOUBLE H;
  LPLPDOUBLE R;
  LPLPINT cells;
  LPINT mask, edges, mcells, hnodes;
  int iter[4],i,j;
  clock_t ticks1, ticks2;

  FILE *sout;
  sout=fopen("smoother.out","wr");//stdout;

  n=_dim;
  me=_metric;
  gr=_generate_data?0:1;
  adp=_adaptive_func;
  theta=_theta;
  miniter=_miniter;
  maxiter=_maxiter;
  miniterBC=_miniterBC;

  if((gr==0)&&(me>1))
  {
    for(i=0;i<40;i++) grid_old[i]=grid[i];
      metr_data_gen(grid, metr, n, me, sout); //generate metric from initial mesh (me=2,3)
  }

  s=_dim;
  N=_mesh.n_nodes();
  ncells=_mesh.n_active_elem();

  //I wish I could do this in readgr... but the way she allocs
  //memory we need to do this here... or pass a lot of things around
  MeshTools::find_hanging_nodes_and_parents(_mesh,_hanging_nodes);

  nedges=_hanging_nodes.size();
  //nedges=0;
  if(s!=n){ fprintf(sout,"Error: dim in input file is inconsistent with dim in .cfg \n"); return _dist_norm; }

  mask=alloc_i_n1(N);
  edges=alloc_i_n1(2*nedges);
  mcells=alloc_i_n1(ncells);
  hnodes=alloc_i_n1(nedges);
  H=alloc_d_n1_n2_n3(ncells,n,n);
  R=alloc_d_n1_n2(N,n);
  cells=alloc_i_n1_n2(ncells,3*n+n%2);
  for(i=0;i<ncells;i++){
     cells[i]=alloc_i_n1(3*n+n%2);
	 if(me>1){
		 H[i]=alloc_d_n1_n2(n,n);
		 for(j=0;j<n;j++) H[i][j]=alloc_d_n1(n);
	 }
  }
  for(i=0;i<N;i++) R[i]=alloc_d_n1(n);

  /*----------initial grid---------*/
  err=readgr(n,N,R,mask,ncells,cells,mcells,nedges,edges,hnodes,sout);
  if(err<0) {fprintf(sout,"Error reading input mesh file\n"); return _dist_norm;}
  if(me>1) err=readmetr(metr,H,ncells,n,sout);
  if(err<0) {fprintf(sout,"Error reading metric file\n"); return _dist_norm;}

  iter[0]=miniter; iter[1]=maxiter; iter[2]=miniterBC;


  /*---------grid optimization--------*/
  fprintf(sout,"Starting Grid Optimization \n");
  ticks1=clock();
  full_smooth(n,N,R,mask,ncells,cells,mcells,nedges,edges,hnodes,theta,iter,me,H,adp,adap,gr,sout);
  ticks2=clock();
  fprintf(sout,"full_smooth took (%d-%d)/%ld = %ld seconds \n",
	  (int)ticks2,(int)ticks1,(long int)CLOCKS_PER_SEC,(long int)(ticks2-ticks1)/CLOCKS_PER_SEC);

  /*---------save result---------*/
  fprintf(sout,"Saving Result \n");
  writegr(n,N,R,mask,ncells,cells,mcells,nedges,edges,hnodes,"fin_grid.dat", 4-me+3*gr, grid_old, sout);

  for(i=0;i<N;i++) free(R[i]);
  for(i=0;i<ncells;i++){
	  if(me>1){
		  for(j=0;j<n;j++) free(H[i][j]);
		  free(H[i]);
	  }
     free(cells[i]);
  }

  free(cells);
  free(H);
  free(R);
  free(mask);
  free(edges);
  free(mcells);
  free(hnodes);
  //fclose(sout);
  libmesh_assert(_dist_norm > 0);

  return _dist_norm;
}

/**
 * save grid
 */
//int VariationalMeshSmoother::writegr(int n, int N, LPLPDOUBLE R, LPINT mask, int ncells, LPLPINT cells,
//            LPINT mcells, int nedges, LPINT edges, LPINT hnodes, char grid[], int me, char grid_old[], FILE *sout)
int VariationalMeshSmoother::writegr(int, int, LPLPDOUBLE R, LPINT, int, LPLPINT,
	    LPINT, int, LPINT, LPINT, const char [], int, const char [], FILE*)
{
  std::cout<<"Starting writegr"<<std::endl;
  int i;

  //Adjust nodal coordinates to new positions
  {
    MeshBase::const_node_iterator       it  = _mesh.nodes_begin();
    const MeshBase::const_node_iterator end = _mesh.nodes_end();

    libmesh_assert(_dist_norm == 0.0);
    _dist_norm=0;
    i=0;
    for (; it != end; ++it)
    {
      double total_dist=0.0;

      //For each node set its X Y [Z] coordinates
      for(unsigned int j=0;j<_dim;j++)
      {
	double distance = R[i][j]-(*(*it))(j);

	//Save the squares of the distance
	total_dist+=Utility::pow<2>(distance);

	(*(*it))(j)=(*(*it))(j)+(distance*_percent_to_move);
      }

      libmesh_assert(total_dist >= 0.0);

      //Add the distance this node moved to the global distance
      _dist_norm+=total_dist;

      i++;
    }

    //Relative "error"
    _dist_norm=sqrt(_dist_norm/_mesh.n_nodes());
  }
  /*
if(me>=3){
	fprintf(stream,"%d \n%d \n%d \n%d \n",n,N,ncells,nedges);

    for(i=0;i<N;i++){//node coordinates
	    for(unsigned int j=0;j<n;j++) fprintf(stream,"%e ",R[i][j]);
	fprintf(stream,"%d \n",mask[i]);
    }
    for(i=0;i<ncells;i++){//cell connectivity
	    int nvert=0;
	    while(cells[i][nvert]>=0){
	       fprintf(stream,"%d ",cells[i][nvert]);
	       nvert++;
	    }
	    for(unsigned int j=nvert;j<3*n+n%2;j++) fprintf(stream,"-1 ");
	    fprintf(stream,"%d \n",mcells[i]);
	}
    //hanging nodes on edges
    for(i=0;i<nedges;i++) fprintf(stream,"%d %d %d \n",edges[2*i],edges[2*i+1],hnodes[i]);
} else {//restoring original grid connectivity when me>1
	int lo, Ncells;
	double d1;
	stream1=fopen(grid_old,"r");
	fscanf(stream1,"%d \n%d \n%d \n%d \n",&lo,&i,&Ncells,&nedges);
	fprintf(stream,"%d \n%d \n%d \n%d \n",n,N,Ncells,nedges);

	for(i=0;i<N;i++){//node coordinates
		for(unsigned int j=0;j<n;j++) {fprintf(stream,"%e ",R[i][j]); fscanf(stream1,"%le ",&d1);}
	fprintf(stream,"%d \n",mask[i]); fscanf(stream1,"%d \n",&lo);
    }
	for(i=0;i<Ncells;i++){
		for(unsigned int j=0;j<=3*n+n%2;j++){
			fscanf(stream1,"%d ",&lo); fprintf(stream,"%d ",lo);
		}
		fscanf(stream1,"\n"); fprintf(stream,"\n");
	}
	for(i=0;i<nedges;i++){
		for(unsigned int j=0;j<3;j++){
	    fscanf(stream1,"%d ",&lo); fprintf(stream,"%d ",lo);
		}
	fscanf(stream1,"\n"); fprintf(stream,"\n");
	}
	fclose(stream1);
}

fclose(stream);
*/
  std::cout<<"Finished writegr"<<std::endl;
  return 0;
}


/**
 * reading grid from input file
 */
// int VariationalMeshSmoother::readgr(int n, int N, LPLPDOUBLE R, LPINT mask, int ncells, LPLPINT cells,
//            LPINT mcells, int nedges, LPINT edges, LPINT hnodes, FILE *sout)
int VariationalMeshSmoother::readgr(int n, int, LPLPDOUBLE R, LPINT mask, int, LPLPINT cells,
	   LPINT mcells, int, LPINT edges, LPINT hnodes, FILE *)
{
  std::cout<<"Sarting readgr"<<std::endl;
  int i;
 //add error messages where format can be inconsistent

  //Find the boundary nodes
  std::vector<bool> on_boundary;
  MeshTools::find_boundary_nodes(_mesh, on_boundary);

  //Grab node coordinates and set mask
  {
    MeshBase::const_node_iterator       it  = _mesh.nodes_begin();
    const MeshBase::const_node_iterator end = _mesh.nodes_end();

    //Only compute the node to elem map once
    std::vector<std::vector<const Elem*> > nodes_to_elem_map;
    MeshTools::build_nodes_to_elem_map(_mesh, nodes_to_elem_map);

    i=0;
    for (; it != end; ++it)
    {
      //For each node grab its X Y [Z] coordinates
      for(unsigned int j=0;j<_dim;j++)
	R[i][j]=(*(*it))(j);

      //Set the Proper Mask
      //Internal nodes are 0
      //Immovable boundary nodes are 1
      //Movable boundary nodes are 2
      if(on_boundary[i])
      {
	//Only look for sliding edge nodes in 2D
	if(_dim == 2)
	{
	  //Find all the nodal neighbors... that is the nodes directly connected
	  //to this node through one edge
	  std::vector<const Node*> neighbors;
	  MeshTools::find_nodal_neighbors(_mesh, *(*it), nodes_to_elem_map, neighbors);

	  std::vector<const Node*>::const_iterator ne = neighbors.begin();
	  std::vector<const Node*>::const_iterator ne_end = neighbors.end();

	  //Grab the x,y coordinates
	  Real x = (*(*it))(0);
	  Real y = (*(*it))(1);

	  //Theta will represent the atan2 angle (meaning with the proper quadrant in mind)
	  //of the neighbor node in a system where the current node is at the origin
	  Real theta = 0;
	  std::vector<Real> thetas;

	  //Calculate the thetas
	  for(; ne != ne_end; ne++)
	  {
	    //Note that the x and y values of this node are subtracted off
	    //this centers the system around this node
	    theta = atan2((*(*ne))(1)-y,(*(*ne))(0)-x);

	    //Save it for later
	    thetas.push_back(theta);
	  }

	  //Assume the node is immovable... then prove otherwise
	  mask[i]=1;

	  //Search through neighbor nodes looking for two that form a straight line with this node
	  for(unsigned int a=0;a<thetas.size()-1;a++)
	  {
	    //Only try each pairing once
	    for(unsigned int b=a+1;b<thetas.size();b++)
	    {
	      //Find if the two neighbor nodes angles are 180 degrees (pi) off of eachother (withing a tolerance)
	      //In order to make this a true movable boundary node... the two that forma  straight line with
	      //it must also be on the boundary
	      if(on_boundary[neighbors[a]->id()] &&
		on_boundary[neighbors[b]->id()] &&
		(
		  (fabs(thetas[a]-(thetas[b] + (libMesh::pi))) < .001) ||
		  (fabs(thetas[a]-(thetas[b] - (libMesh::pi))) < .001)
		)
		)
	      {
		//if((*(*it))(1) > 0.25 || (*(*it))(0) > .7 || (*(*it))(0) < .01)
		  mask[i]=2;
	      }

	    }
	  }
	}
	else //In 3D set all boundary nodes to be fixed
	  mask[i]=1;
      }
      else
	mask[i]=0;  //Internal Node

      //std::cout<<"Node: "<<i<<"  Mask: "<<mask[i]<<std::endl;
      i++;
    }
  }

  // Grab the connectivity
  // FIXME: Generalize this!
  {
    MeshBase::const_element_iterator it  = _mesh.active_elements_begin();
    const MeshBase::const_element_iterator end = _mesh.active_elements_end();

    i=0;
    for (; it != end; ++it)
    {
      //Keep track of the number of nodes
      //there must be 6 for 2D
      //10 for 3D
      //If there are actually less than that -1 must be used
      //to fill out the rest
      int num=0;
      int num_necessary = 0;

      if (_dim == 2)
	num_necessary = 6;
      else
	num_necessary = 10;

      if(_dim == 2)
      {
	switch((*it)->n_vertices())
	{
	  //Grab nodes that do exist
	  case 3:  //Tri
	    for(unsigned int k=0;k<(*it)->n_vertices();k++)
	      cells[i][k]=(*it)->node(k);
	    num=(*it)->n_vertices();
	    break;
	  case 4:  //Quad 4
	    cells[i][0]=(*it)->node(0);
	    cells[i][1]=(*it)->node(1);
	    cells[i][2]=(*it)->node(3); //Note that 2 and 3 are switched!
	    cells[i][3]=(*it)->node(2);
	    num=4;
	    break;
	  default:
	    libmesh_assert(false);
	}
      }
      else
      {
	//Grab nodes that do exist
	switch((*it)->n_vertices())
	{
	  case 4:  //Tet 4
	    for(unsigned int k=0;k<(*it)->n_vertices();k++)
	      cells[i][k]=(*it)->node(k);
	    num=(*it)->n_vertices();
	    break;
	  case 8:  //Hex 8f
	    cells[i][0]=(*it)->node(0);
	    cells[i][1]=(*it)->node(1);
	    cells[i][2]=(*it)->node(3); //Note that 2 and 3 are switched!
	    cells[i][3]=(*it)->node(2);

	    cells[i][4]=(*it)->node(4);
	    cells[i][5]=(*it)->node(5);
	    cells[i][6]=(*it)->node(7); //Note that 6 and 7 are switched!
	    cells[i][7]=(*it)->node(6);
	    num=8;
	  default:
	    libmesh_assert(false);
	}
      }

      //Fill in the rest with -1
      for(int j=num;j<3*n+n%2;j++)
	cells[i][j]=-1;

      //Mask it with 0 to state that this is an active element
      //FIXME: Could be something other than zero
      mcells[i]=0;

      i++;
    }
  }

  //Grab hanging node connectivity
  {
    std::map<unsigned int, std::vector<unsigned int> >::iterator it = _hanging_nodes.begin();
    std::map<unsigned int, std::vector<unsigned int> >::iterator end = _hanging_nodes.end();

    for(i=0;it!=end;it++)