mesh_modification.c

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

// $Id: mesh_modification.C 2822 2008-05-01 20:43:09Z roystgnr $

// 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



// C++ includes
#include <cmath> // for std::acos()
#include <algorithm>
#include <map>

// Local includes
#include "boundary_info.h"
#include "face_tri3.h"
#include "face_tri6.h"
#include "libmesh_logging.h"
#include "location_maps.h"
#include "mesh_communication.h"
#include "mesh_modification.h"
#include "mesh_tools.h"
#include "parallel.h"
#include "string_to_enum.h"
#include "unstructured_mesh.h"


// ------------------------------------------------------------
// MeshTools::Modification functions for mesh modification
void MeshTools::Modification::distort (MeshBase& mesh,
				       const Real factor,
				       const bool perturb_boundary)
{
  libmesh_assert (mesh.mesh_dimension() != 1);
  libmesh_assert (mesh.n_nodes());
  libmesh_assert (mesh.n_elem());
  libmesh_assert ((factor >= 0.) && (factor <= 1.));

  START_LOG("distort()", "MeshTools::Modification");



  // First find nodes on the boundary and flag them
  // so that we don't move them
  // on_boundary holds false (not on boundary) and true (on boundary)
  std::vector<bool> on_boundary (mesh.n_nodes(), false);
  
  if (!perturb_boundary) MeshTools::find_boundary_nodes (mesh, on_boundary);

  // Now calculate the minimum distance to
  // neighboring nodes for each node.
  // hmin holds these distances.
  std::vector<float> hmin (mesh.n_nodes(), 1.e20);

  MeshBase::element_iterator       el  = mesh.active_elements_begin();
  const MeshBase::element_iterator end = mesh.active_elements_end(); 

  for (; el!=end; ++el)
    for (unsigned int n=0; n<(*el)->n_nodes(); n++)
      hmin[(*el)->node(n)] = std::min(hmin[(*el)->node(n)],
                                      static_cast<float>((*el)->hmin()));
  
  
  // Now actually move the nodes
  {
    const unsigned int seed = 123456;
    
    // seed the random number generator
    std::srand(seed);
    
    // If the node is on the boundary or
    // the node is not used by any element (hmin[n]<1.e20)
    // then we should not move it.
    // [Note: Testing for (in)equality might be wrong
    // (different types, namely float and double)]
    for (unsigned int n=0; n<mesh.n_nodes(); n++)
      if (!on_boundary[n] && (hmin[n] < 1.e20) )
	{
	  // the direction, random but unit normalized
	  
	  Point dir( static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
		     static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
		     ((mesh.mesh_dimension() == 3) ?
		      static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) :
		      0.)
		     );
	  
	  dir(0) = (dir(0)-.5)*2.;
	  dir(1) = (dir(1)-.5)*2.;
	  if (mesh.mesh_dimension() == 3)
	    dir(2) = (dir(2)-.5)*2.;
	  
	  dir = dir.unit();

          mesh.node(n)(0) += dir(0)*factor*hmin[n];
          mesh.node(n)(1) += dir(1)*factor*hmin[n];
          if (mesh.mesh_dimension() == 3)
            mesh.node(n)(2) += dir(2)*factor*hmin[n];
	}
  }


  // All done  
  STOP_LOG("distort()", "MeshTools::Modification");
}



void MeshTools::Modification::translate (MeshBase& mesh,
					 const Real xt,
					 const Real yt,
					 const Real zt)
{
  const Point p(xt, yt, zt);

  for (unsigned int n=0; n<mesh.n_nodes(); n++)
    mesh.node(n) += p;
}


// void MeshTools::Modification::rotate2D (MeshBase& mesh,
//                                         const Real alpha)
// {
//   libmesh_assert (mesh.mesh_dimension() != 1);

//   const Real pi = std::acos(-1);
//   const Real  a = alpha/180.*pi;
//   for (unsigned int n=0; n<mesh.n_nodes(); n++)
//     {
//       const Point p = mesh.node(n);
//       const Real  x = p(0);
//       const Real  y = p(1);
//       const Real  z = p(2);
//       mesh.node(n) = Point(std::cos(a)*x - std::sin(a)*y,
//                            std::sin(a)*x + std::cos(a)*y,
//                            z);
//     }

// }



void MeshTools::Modification::rotate (MeshBase& mesh,
				      const Real phi,
				      const Real theta,
				      const Real psi)
{
  libmesh_assert (mesh.mesh_dimension() != 1);

  const Real pi = std::acos(-1.);
  const Real  p = -phi/180.*pi;
  const Real  t = -theta/180.*pi;
  const Real  s = -psi/180.*pi;
  const Real sp = std::sin(p), cp = std::cos(p);
  const Real st = std::sin(t), ct = std::cos(t);
  const Real ss = std::sin(s), cs = std::cos(s);

  for (unsigned int n=0; n<mesh.n_nodes(); n++)
    {
      const Point p = mesh.node(n);
      const Real  x = p(0);
      const Real  y = p(1);
      const Real  z = p(2);
      mesh.node(n) = Point(( cp*cs-sp*ct*ss)*x + ( sp*cs+cp*ct*ss)*y + (st*ss)*z,
                           (-cp*ss-sp*ct*cs)*x + (-sp*st+cp*ct*cs)*y + (st*cs)*z,
                           ( sp*st)*x          + (-cp*st)*y          + (ct)*z   );
    }
}


void MeshTools::Modification::scale (MeshBase& mesh,
				     const Real xs,
				     const Real ys,
				     const Real zs)
{
  const Real x_scale = xs;
  Real y_scale       = ys;
  Real z_scale       = zs;
  
  if (ys == 0.)
    {
      libmesh_assert (zs == 0.);

      y_scale = z_scale = x_scale;
    }

  // Scale the x coordinate in all dimensions
  for (unsigned int n=0; n<mesh.n_nodes(); n++)
    mesh.node(n)(0) *= x_scale;


  // Only scale the y coordinate in 2 and 3D
  if (mesh.spatial_dimension() > 1)
    {

      for (unsigned int n=0; n<mesh.n_nodes(); n++)
	mesh.node(n)(1) *= y_scale;

      // Only scale the z coordinate in 3D
      if (mesh.spatial_dimension() == 3)
	{
	  for (unsigned int n=0; n<mesh.n_nodes(); n++)
	    mesh.node(n)(2) *= z_scale;
	}
    }
}




// ------------------------------------------------------------
// UnstructuredMesh class member functions for mesh modification
void UnstructuredMesh::all_first_order ()
{
  /*
   * when the mesh is not prepared,
   * at least renumber the nodes and 
   * elements, so that the node ids
   * are correct
   */
  if (!this->_is_prepared)
    this->renumber_nodes_and_elements ();

  START_LOG("all_first_order()", "Mesh");

  /**
   * Loop over the high-ordered elements.
   * First make sure they _are_ indeed high-order, and then replace
   * them with an equivalent first-order element.
   */
  const_element_iterator endit = elements_end();
  for (const_element_iterator it = elements_begin();
       it != endit; ++it)
    {
      Elem* so_elem = *it;

      libmesh_assert (so_elem != NULL);

      /*
       * build the first-order equivalent, add to
       * the new_elements list.
       */
      Elem *newparent = so_elem->parent();
      Elem* lo_elem = Elem::build
        (Elem::first_order_equivalent_type
          (so_elem->type()), newparent).release();

#ifdef ENABLE_AMR
      /*
       * Add this element to it's parent if it has one
       */
      if (newparent)
        newparent->add_child(lo_elem);

      /*
       * Reset the parent links of any child elements
       */
      if (so_elem->has_children())
        {
          for (unsigned int c=0; c != so_elem->n_children(); ++c)
            so_elem->child(c)->set_parent(lo_elem);
        }

      /*
       * Copy as much data to the new element as makes sense
       */
      lo_elem->set_p_level(so_elem->p_level());
      lo_elem->set_refinement_flag(so_elem->refinement_flag());
      lo_elem->set_p_refinement_flag(so_elem->p_refinement_flag());
#endif

      libmesh_assert (lo_elem->n_vertices() == so_elem->n_vertices());

      /*
       * By definition the vertices of the linear and
       * second order element are identically numbered.
       * transfer these.
       */
      for (unsigned int v=0; v < so_elem->n_vertices(); v++)
	lo_elem->set_node(v) = so_elem->get_node(v);

      /**
       * If the second order element had any boundary conditions they
       * should be transfered to the first-order element.  The old
       * boundary conditions will be removed from the BoundaryInfo
       * data structure by insert_elem.
       */
      libmesh_assert (lo_elem->n_sides() == so_elem->n_sides());
	
      for (unsigned int s=0; s<so_elem->n_sides(); s++)
	{
	  const short int boundary_id =
	    this->boundary_info->boundary_id (so_elem, s);
	    
	  if (boundary_id != this->boundary_info->invalid_id)
	    this->boundary_info->add_side (lo_elem, s, boundary_id);
	}

      /*
       * The new first-order element is ready.
       * Inserting it into the mesh will replace and delete
       * the second-order element.
       */
      lo_elem->set_id(so_elem->id());
      this->insert_elem(lo_elem);
    }

  STOP_LOG("all_first_order()", "Mesh");

  // delete or renumber nodes, etc
  this->prepare_for_use();
}



void UnstructuredMesh::all_second_order (const bool full_ordered)
{
  // This function must be run on all processors at once
  parallel_only();

  /*
   * when the mesh is not prepared,
   * at least renumber the nodes and 
   * elements, so that the node ids
   * are correct
   */
  if (!this->_is_prepared)
    this->renumber_nodes_and_elements ();

  /*
   * If the mesh is empty
   * then we have nothing to do
   */
  if (!this->n_elem())
    return;
 
  /*
   * If the mesh is already second order
   * then we have nothing to do.
   * We have to test for this in a round-about way to avoid
   * a bug on distributed parallel meshes with more processors
   * than elements.
   */
  bool already_second_order = false;
  if (this->elements_begin() != this->elements_end() &&
      (*(this->elements_begin()))->default_order() != FIRST)
    already_second_order = true;
  Parallel::max(already_second_order);
  if (already_second_order)
    return;

  START_LOG("all_second_order()", "Mesh");

  /*
   * this map helps in identifying second order
   * nodes.  Namely, a second-order node:
   * - edge node
   * - face node
   * - bubble node
   * is uniquely defined through a set of adjacent
   * vertices.  This set of adjacent vertices is
   * used to identify already added higher-order
   * nodes.  We are safe to use node id's since we
   * make sure that these are correctly numbered.
   */
  std::map<std::vector<unsigned int>, Node*> adj_vertices_to_so_nodes;

  /*
   * for speed-up of the \p add_point() method, we
   * can reserve memory.  Guess the number of additional
   * nodes for different dimensions
   */
  switch (this->mesh_dimension())
  {
    case 1:
      /*
       * in 1D, there can only be order-increase from Edge2
       * to Edge3.  Something like 1/2 of n_nodes() have
       * to be added
       */
      this->reserve_nodes(static_cast<unsigned int>(1.5*this->n_nodes()));
      break;

    case 2:
      /*
       * in 2D, either refine from Tri3 to Tri6 (double the nodes)
       * or from Quad4 to Quad8 (again, double) or Quad9 (2.25 that much)
       */
      this->reserve_nodes(static_cast<unsigned int>(2*this->n_nodes()));
      break;


    case 3:
      /*
       * in 3D, either refine from Tet4 to Tet10 (factor = 2.5) up to
       * Hex8 to Hex27 (something  > 3).  Since in 3D there _are_ already
       * quite some nodes, and since we do not want to overburden the memory by
       * a too conservative guess, use the lower bound
       */
      this->reserve_nodes(static_cast<unsigned int>(2.5*this->n_nodes()));
      break;
	
    default:
      // Hm?
      libmesh_error();
  }