mesh_modification.c

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

  /*
   * form a vector that will hold the node id's of
   * the vertices that are adjacent to the son-th
   * second-order node.  Pull this outside of the
   * loop so that silly compilers don't repeatedly
   * create and destroy the vector.
   */
  std::vector<unsigned int> adjacent_vertices_ids;

  /**
   * Loop over the low-ordered elements in the _elements vector.
   * First make sure they _are_ indeed low-order, and then replace
   * them with an equivalent second-order element.  Don't
   * forget to delete the low-order element, or else it will leak!
   */
  const_element_iterator endit = elements_end();
  for (const_element_iterator it = elements_begin();
       it != endit; ++it)
    {
      // the linear-order element
      Elem* lo_elem = *it;

      libmesh_assert (lo_elem != NULL);

      // make sure it is linear order
      if (lo_elem->default_order() != FIRST)
        {	  
	  std::cerr << "ERROR: This is not a linear element: type=" 
		    << lo_elem->type() << std::endl;
	  libmesh_error();
	}

      // this does _not_ work for refined elements
      libmesh_assert (lo_elem->level () == 0);

      /*
       * build the second-order equivalent, add to
       * the new_elements list.  Note that this here
       * is the only point where \p full_ordered
       * is necessary.  The remaining code works well
       * for either type of seconrd-order equivalent, e.g.
       * Hex20 or Hex27, as equivalents for Hex8
       */
      Elem* so_elem = 
	Elem::build (Elem::second_order_equivalent_type(lo_elem->type(), 
							full_ordered) ).release();

      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 < lo_elem->n_vertices(); v++)
	so_elem->set_node(v) = lo_elem->get_node(v);

      /*
       * Now handle the additional mid-side nodes.  This
       * is simply handled through a map that remembers
       * the already-added nodes.  This map maps the global
       * ids of the vertices (that uniquely define this 
       * higher-order node) to the new node. 
       * Notation: son = second-order node
       */
      const unsigned int son_begin = so_elem->n_vertices();
      const unsigned int son_end   = so_elem->n_nodes();
      

      for (unsigned int son=son_begin; son<son_end; son++)
        {
	  const unsigned int n_adjacent_vertices =
	    so_elem->n_second_order_adjacent_vertices(son);

	  adjacent_vertices_ids.resize(n_adjacent_vertices);
	  
	  for (unsigned int v=0; v<n_adjacent_vertices; v++)
	    adjacent_vertices_ids[v] =
	      so_elem->node( so_elem->second_order_adjacent_vertex(son,v) );

	  /*
	   * \p adjacent_vertices_ids is now in order of the current
	   * side.  sort it, so that comparisons  with the 
	   * \p adjacent_vertices_ids created through other elements' 
	   * sides can match
	   */
	  std::sort(adjacent_vertices_ids.begin(),
		    adjacent_vertices_ids.end());


	  // does this set of vertices already has a mid-node added?
	  std::pair<std::map<std::vector<unsigned int>, Node*>::iterator,
                    std::map<std::vector<unsigned int>, Node*>::iterator>	    
	    pos = adj_vertices_to_so_nodes.equal_range (adjacent_vertices_ids);

	  // no, not added yet
	  if (pos.first == pos.second)
	    {
	      /*
	       * for this set of vertices, there is no 
	       * second_order node yet.  Add it.
	       *
	       * compute the location of the new node as
	       * the average over the adjacent vertices.
	       */
	      Point new_location = this->point(adjacent_vertices_ids[0]);
	      for (unsigned int v=1; v<n_adjacent_vertices; v++)
		new_location += this->point(adjacent_vertices_ids[v]);

	      new_location /= static_cast<Real>(n_adjacent_vertices);
	      
	      /* Add the new point to the mesh, giving it a globally
               * well-defined processor id.
	       */
	      Node* so_node = this->add_point
                (new_location, DofObject::invalid_id,
                this->node(adjacent_vertices_ids[0]).processor_id());

	      /* 
	       * insert the new node with its defining vertex
	       * set into the map, and relocate pos to this
	       * new entry, so that the so_elem can use
	       * \p pos for inserting the node
	       */
	      adj_vertices_to_so_nodes.insert(pos.first,
					      std::make_pair(adjacent_vertices_ids,
							     so_node));

	      so_elem->set_node(son) = so_node;
	    }
	  // yes, already added.
	  else
	    {
	      libmesh_assert (pos.first->second != NULL);
	      
	      so_elem->set_node(son) = pos.first->second;
	    }
	}


      /**
       * If the linear element had any boundary conditions they
       * should be transfered to the second-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<lo_elem->n_sides(); s++)
	{
	  const short int boundary_id =
	    this->boundary_info->boundary_id (lo_elem, s);
	    
	  if (boundary_id != this->boundary_info->invalid_id)
	    this->boundary_info->add_side (so_elem, s, boundary_id);
	}

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

  // we can clear the map
  adj_vertices_to_so_nodes.clear();


  STOP_LOG("all_second_order()", "Mesh");

  // In a ParallelMesh our ghost node processor ids may be bad and
  // the ids of nodes touching remote elements may be inconsistent.
  // Fix them.
  if (!this->is_serial())
    {
      LocationMap<Node> loc_map;
      MeshCommunication().make_nodes_parallel_consistent
        (*this, loc_map);
    }

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






void MeshTools::Modification::all_tri (MeshBase& mesh)
{
  // We store pointers to the newly created elements in a vector
  // until they are ready to be added to the mesh.  This is because
  // adding new elements on the fly can cause reallocation and invalidation
  // of existing iterators.
  std::vector<Elem*> new_elements;
  new_elements.reserve (2*mesh.n_active_elem());

  // If the original mesh has boundary data, we carry that over
  // to the new mesh with triangular elements.
  const bool mesh_has_boundary_data = (mesh.boundary_info->n_boundary_ids() > 0);

  // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
  std::vector<Elem*> new_bndry_elements;
  std::vector<unsigned short int> new_bndry_sides;
  std::vector<short int> new_bndry_ids;
  
  // Iterate over the elements, splitting QUADS into
  // pairs of conforming triangles.
  {
    MeshBase::element_iterator       el  = mesh.elements_begin();
    const MeshBase::element_iterator end = mesh.elements_end(); 

    for (; el!=end; ++el)
      {
	const ElemType etype = (*el)->type();

	// all_tri currently only works on coarse meshes
	libmesh_assert ((*el)->parent() == NULL);

	// We split the quads using the shorter of the two diagonals
	// to maintain the best angle properties.
	bool edge_swap = false;

	// True if we actually split the current element.
	bool split_elem = false;

	// The two new triangular elements we will split the quad into.
	Elem* tri0 = NULL;
	Elem* tri1 = NULL;

	
	switch (etype)
	  {
	  case QUAD4:
	    {
	      split_elem = true;
	      
	      tri0 = new Tri3;
	      tri1 = new Tri3;
	      
	      // Check for possible edge swap
	      if (((*el)->point(0) - (*el)->point(2)).size() <
		  ((*el)->point(1) - (*el)->point(3)).size())
		{	      
		  tri0->set_node(0) = (*el)->get_node(0);
		  tri0->set_node(1) = (*el)->get_node(1);
		  tri0->set_node(2) = (*el)->get_node(2);
		
		  tri1->set_node(0) = (*el)->get_node(0);
		  tri1->set_node(1) = (*el)->get_node(2);
		  tri1->set_node(2) = (*el)->get_node(3);
		}
	    
	      else
		{
		  edge_swap=true;
		  
		  tri0->set_node(0) = (*el)->get_node(0);
		  tri0->set_node(1) = (*el)->get_node(1);
		  tri0->set_node(2) = (*el)->get_node(3);
		
		  tri1->set_node(0) = (*el)->get_node(1);
		  tri1->set_node(1) = (*el)->get_node(2);
		  tri1->set_node(2) = (*el)->get_node(3);
		}

	      
	      break;
	    }
      
	  case QUAD8:
	    {
	      split_elem =  true;
	      
	      tri0 = new Tri6;
	      tri1 = new Tri6;
	  
	      Node* new_node = mesh.add_point((mesh.node((*el)->node(0)) +
					       mesh.node((*el)->node(1)) +
					       mesh.node((*el)->node(2)) +
					       mesh.node((*el)->node(3)))*.25
					       );
	  
	      // Check for possible edge swap
	      if (((*el)->point(0) - (*el)->point(2)).size() <
		  ((*el)->point(1) - (*el)->point(3)).size())
		{	      
		  tri0->set_node(0) = (*el)->get_node(0);
		  tri0->set_node(1) = (*el)->get_node(1);
		  tri0->set_node(2) = (*el)->get_node(2);
		  tri0->set_node(3) = (*el)->get_node(4);
		  tri0->set_node(4) = (*el)->get_node(5);
		  tri0->set_node(5) = new_node;
	      
		  tri1->set_node(0) = (*el)->get_node(0);
		  tri1->set_node(1) = (*el)->get_node(2);
		  tri1->set_node(2) = (*el)->get_node(3);
		  tri1->set_node(3) = new_node;
		  tri1->set_node(4) = (*el)->get_node(6);
		  tri1->set_node(5) = (*el)->get_node(7);

		}
	  
	      else
		{
		  edge_swap=true;
		  
		  tri0->set_node(0) = (*el)->get_node(3);
		  tri0->set_node(1) = (*el)->get_node(0);
		  tri0->set_node(2) = (*el)->get_node(1);
		  tri0->set_node(3) = (*el)->get_node(7);
		  tri0->set_node(4) = (*el)->get_node(4);
		  tri0->set_node(5) = new_node;
	      
		  tri1->set_node(0) = (*el)->get_node(1);
		  tri1->set_node(1) = (*el)->get_node(2);
		  tri1->set_node(2) = (*el)->get_node(3);
		  tri1->set_node(3) = (*el)->get_node(5);
		  tri1->set_node(4) = (*el)->get_node(6);
		  tri1->set_node(5) = new_node;
		}
	  
	      break;
	    }
      
	  case QUAD9:
	    {
	      split_elem =  true;
	      
	      tri0 = new Tri6;
	      tri1 = new Tri6;

	      // Check for possible edge swap
	      if (((*el)->point(0) - (*el)->point(2)).size() <
		  ((*el)->point(1) - (*el)->point(3)).size())
		{	      
		  tri0->set_node(0) = (*el)->get_node(0);
		  tri0->set_node(1) = (*el)->get_node(1);
		  tri0->set_node(2) = (*el)->get_node(2);
		  tri0->set_node(3) = (*el)->get_node(4);
		  tri0->set_node(4) = (*el)->get_node(5);
		  tri0->set_node(5) = (*el)->get_node(8);
	      
		  tri1->set_node(0) = (*el)->get_node(0);
		  tri1->set_node(1) = (*el)->get_node(2);
		  tri1->set_node(2) = (*el)->get_node(3);
		  tri1->set_node(3) = (*el)->get_node(8);
		  tri1->set_node(4) = (*el)->get_node(6);
		  tri1->set_node(5) = (*el)->get_node(7);
		}

	      else
		{
		  edge_swap=true;
		  
		  tri0->set_node(0) = (*el)->get_node(0);
		  tri0->set_node(1) = (*el)->get_node(1);
		  tri0->set_node(2) = (*el)->get_node(3);
		  tri0->set_node(3) = (*el)->get_node(4);
		  tri0->set_node(4) = (*el)->get_node(8);
		  tri0->set_node(5) = (*el)->get_node(7);
	      
		  tri1->set_node(0) = (*el)->get_node(1);
		  tri1->set_node(1) = (*el)->get_node(2);
		  tri1->set_node(2) = (*el)->get_node(3);
		  tri1->set_node(3) = (*el)->get_node(5);
		  tri1->set_node(4) = (*el)->get_node(6);
		  tri1->set_node(5) = (*el)->get_node(8);
		}
	    
	      break;
	    }
          // No need to split elements that are already triangles
          case TRI3:
          case TRI6:
            continue;
          // Try to ignore non-2D elements for now
	  default:
	    {
	      std::cerr << "Warning, encountered non-2D element "
                        << Utility::enum_to_string<ElemType>(etype)
			<< " in MeshTools::Modification::all_tri(), hope that's OK..."
			<< std::endl;
	    }
	  } // end switch (etype)

	

	if (split_elem)
	  {
	    // If this mesh has boundary data, be sure the correct
	    // ID's are also set for tri0 and tri1.
	    if (mesh_has_boundary_data)
	      {
		for (unsigned int sn=0; sn<(*el)->n_sides(); ++sn)
		  {
		    short int b_id = mesh.boundary_info->boundary_id(*el, sn);

		    if (b_id != BoundaryInfo::invalid_id)
		      {
			// Add the boundary ID to the list of new boundary ids
			new_bndry_ids.push_back(b_id);
			  
			// Convert the boundary side information of the old element to
			// boundary side information for the new element.
			if (!edge_swap)
			  {
			    switch (sn)
			      {
			      case 0:
				{