mesh_refinement.c

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

			//         was originally flagged for refinement. If it
			//         wasn't flagged already we need to repeat
			//         this process.			 
			else if ((neighbor->level()+1) == my_level)
			  {
			    if (neighbor->refinement_flag() != Elem::REFINE)
			      {
				neighbor->set_refinement_flag(Elem::REFINE);
                                if (neighbor->parent())
				  neighbor->parent()->set_refinement_flag(Elem::INACTIVE);
				compatible_with_coarsening = false;
				level_one_satisfied = false; 
			      }
			  }
#ifdef DEBUG		
			
			// Sanity check. We should never get into a
			// case when our neighbot is more than one
			// level away.			 
			else if ((neighbor->level()+1) < my_level)
			  {
			    libmesh_error();
			  }
			
			
			// Note that the only other possibility is that the
			// neighbor is already refined, in which case it isn't
			// active and we should never get here.			 
			else
			  {
			    libmesh_error();
			  }
#endif		
		      } 
                  }
	      }
	    if (elem->p_refinement_flag() == Elem::REFINE)  // If the element is active and the
                                                            // p refinement flag is set    
	      {
		const unsigned int my_p_level = elem->p_level();

		for (unsigned int side=0; side != elem->n_sides(); side++)
                  {
                    Elem *neighbor = elem->neighbor(side);
		    if (neighbor != NULL &&      // I have a
		        neighbor != remote_elem) // neighbor here
		      if (neighbor->active()) // and it is active
		        {
                          if (neighbor->p_level() < my_p_level &&
                               neighbor->p_refinement_flag() != Elem::REFINE)
			    {
			      neighbor->set_p_refinement_flag(Elem::REFINE);
			      level_one_satisfied = false;
			      compatible_with_coarsening = false;
			    }
                           if (neighbor->p_level() == my_p_level &&
                              neighbor->p_refinement_flag() == Elem::COARSEN)
			    {
			      neighbor->set_p_refinement_flag(Elem::DO_NOTHING);
			      level_one_satisfied = false;
			      compatible_with_coarsening = false;
			    }
		        } 
		      else // I have an inactive neighbor
		        {
                           libmesh_assert(neighbor->has_children());
	                   for (unsigned int c=0; c!=neighbor->n_children(); c++)
                             {
                               Elem *subneighbor = neighbor->child(c);
                               if (subneighbor == remote_elem)
                                 continue;
                               if (subneighbor->active() &&
                                   subneighbor->has_neighbor(elem))
                                 if (subneighbor->p_level() < my_p_level &&
                                     subneighbor->p_refinement_flag() != Elem::REFINE)
			           {
                                     // We should already be level one
                                     // compatible
                                     libmesh_assert(subneighbor->p_level() + 2u > 
                                            my_p_level);
			             subneighbor->set_p_refinement_flag(Elem::REFINE);
			             level_one_satisfied = false;
			             compatible_with_coarsening = false;
			           }
                                 if (subneighbor->p_level() == my_p_level &&
                                     subneighbor->p_refinement_flag() == Elem::COARSEN)
			           {
			             subneighbor->set_p_refinement_flag(Elem::DO_NOTHING);
			             level_one_satisfied = false;
			             compatible_with_coarsening = false;
			           }
                             }
		        } 
                  }
	      }
            }
	}
      
      while (!level_one_satisfied);
    } // end if (_maintain_level_one)

  // If we're not compatible on one processor, we're globally not
  // compatible
  Parallel::min(compatible_with_coarsening);
  
  STOP_LOG ("make_refinement_compatible()", "MeshRefinement");

  return compatible_with_coarsening;
}




bool MeshRefinement::_coarsen_elements ()
{
  // This function must be run on all processors at once
  parallel_only();

  START_LOG ("_coarsen_elements()", "MeshRefinement");

  // Flag indicating if this call actually changes the mesh
  bool mesh_changed = false;
  
  // Clear the unused_elements data structure.
  // The elements have been packed since it was built,
  // so there are _no_ unused elements.  We cannot trust
  // any iterators currently in this data structure.
  // _unused_elements.clear();

  MeshBase::element_iterator       it  = _mesh.elements_begin();
  const MeshBase::element_iterator end = _mesh.elements_end();

  // Loop over the elements.   
  for ( ; it != end; ++it)
    {
      Elem* elem = *it;

      // Not necessary when using elem_iterator
      // libmesh_assert (elem != NULL);
      
      // active elements flagged for coarsening will
      // no longer be deleted until MeshRefinement::contract()
      if (elem->refinement_flag() == Elem::COARSEN)
	{
	  // Huh?  no level-0 element should be active
	  // and flagged for coarsening.
	  libmesh_assert (elem->level() != 0);
	  
	  // Remove this element from any neighbor
	  // lists that point to it.
	  elem->nullify_neighbors();

	  // Remove any boundary information associated
	  // with this element
	  _mesh.boundary_info->remove (elem);
	      
	  // Add this iterator to the _unused_elements
	  // data structure so we might fill it.
	  // The _unused_elements optimization is currently off.
	  // _unused_elements.push_back (it);

	  // Don't delete the element until
	  // MeshRefinement::contract()
	  // _mesh.delete_elem(elem);

	  // the mesh has certainly changed
	  mesh_changed = true;
        }
      
      // inactive elements flagged for coarsening
      // will become active
      else if (elem->refinement_flag() == Elem::COARSEN_INACTIVE)
	{
	  elem->coarsen();
	  libmesh_assert (elem->active());

	  // the mesh has certainly changed
	  mesh_changed = true;
	}
      if (elem->p_refinement_flag() == Elem::COARSEN)
        {
          if (elem->p_level() > 0)
            {
              elem->set_p_refinement_flag(Elem::JUST_COARSENED);
              elem->set_p_level(elem->p_level() - 1);
	      mesh_changed = true;
            }
          else
            {
              elem->set_p_refinement_flag(Elem::DO_NOTHING);
            }
        }
    }  

  // If the mesh changed on any processor, it changed globally
  Parallel::max(mesh_changed);
  // And we may need to update ParallelMesh values reflecting the changes
  if (mesh_changed)
    _mesh.update_parallel_id_counts();

  // Node processor ids may need to change if an element of that id
  // was coarsened away
  if (mesh_changed && !_mesh.is_serial())
    {
      // Update the _new_nodes_map so that processors can
      // find requested nodes
      this->update_nodes_map ();

      MeshCommunication().make_nodes_parallel_consistent
        (_mesh, _new_nodes_map);

      // Clear the _new_nodes_map
      this->clear();
  
#ifdef DEBUG
      MeshTools::libmesh_assert_valid_node_procids(_mesh);
#endif
    }

  STOP_LOG ("_coarsen_elements()", "MeshRefinement");

  return mesh_changed;
}



bool MeshRefinement::_refine_elements ()
{
  // This function must be run on all processors at once
  parallel_only();

  // Update the _new_nodes_map so that elements can
  // find nodes to connect to.
  this->update_nodes_map ();

  START_LOG ("_refine_elements()", "MeshRefinement");

  // Iterate over the elements, counting the elements
  // flagged for h refinement.
  unsigned int n_elems_flagged = 0;

  MeshBase::element_iterator       it  = _mesh.elements_begin();
  const MeshBase::element_iterator end = _mesh.elements_end();

  for (; it != end; ++it)
    {
      Elem* elem = *it;
      if (elem->refinement_flag() == Elem::REFINE)
	n_elems_flagged++;
    }

  // Construct a local vector of Elem* which have been
  // previously marked for refinement.  We reserve enough
  // space to allow for every element to be refined.
  std::vector<Elem*> local_copy_of_elements;
  local_copy_of_elements.reserve(n_elems_flagged);

  // Iterate over the elements, looking for elements
  // flagged for refinement.
  for (it = _mesh.elements_begin(); it != end; ++it)
    {
      Elem* elem = *it;
      if (elem->refinement_flag() == Elem::REFINE)
	local_copy_of_elements.push_back(elem);
      if (elem->p_refinement_flag() == Elem::REFINE &&
          elem->active())
        {
	  elem->set_p_level(elem->p_level()+1);
	  elem->set_p_refinement_flag(Elem::JUST_REFINED);
        }
    }

  // Now iterate over the local copies and refine each one.
  // This may resize the mesh's internal container and invalidate
  // any existing iterators.
  
  for (unsigned int e = 0; e != local_copy_of_elements.size(); ++e)
    local_copy_of_elements[e]->refine(*this);

  // The mesh changed if there were elements h refined
  bool mesh_changed = !local_copy_of_elements.empty();

  // If the mesh changed on any processor, it changed globally
  Parallel::max(mesh_changed);

  // And we may need to update ParallelMesh values reflecting the changes
  if (mesh_changed)
    _mesh.update_parallel_id_counts();

  if (mesh_changed && !_mesh.is_serial())
    {
      MeshCommunication().make_elems_parallel_consistent (_mesh);
      MeshCommunication().make_nodes_parallel_consistent
        (_mesh, _new_nodes_map);
#ifdef DEBUG
      ParallelMesh *pmesh = dynamic_cast<ParallelMesh *>(&_mesh);
      if (pmesh)
        pmesh->libmesh_assert_valid_parallel_ids();
#endif
    }
  
  // Clear the _new_nodes_map and _unused_elements data structures.
  this->clear();
  
  STOP_LOG ("_refine_elements()", "MeshRefinement");

  return mesh_changed;
}



void MeshRefinement::uniformly_p_refine (unsigned int n)
{
  // Refine n times
  for (unsigned int rstep=0; rstep<n; rstep++)
    {
      // P refine all the active elements
      MeshBase::element_iterator       elem_it  = _mesh.active_elements_begin();
      const MeshBase::element_iterator elem_end = _mesh.active_elements_end(); 

      for ( ; elem_it != elem_end; ++elem_it)
        {
	  (*elem_it)->set_p_level((*elem_it)->p_level()+1);
	  (*elem_it)->set_p_refinement_flag(Elem::JUST_REFINED);
        }
    }
}



void MeshRefinement::uniformly_p_coarsen (unsigned int n)
{
  // Coarsen p times
  for (unsigned int rstep=0; rstep<n; rstep++)
    {
      // P coarsen all the active elements
      MeshBase::element_iterator       elem_it  = _mesh.active_elements_begin();
      const MeshBase::element_iterator elem_end = _mesh.active_elements_end(); 

      for ( ; elem_it != elem_end; ++elem_it)
        {
          if ((*elem_it)->p_level() > 0)
            {
	      (*elem_it)->set_p_level((*elem_it)->p_level()-1);
	      (*elem_it)->set_p_refinement_flag(Elem::JUST_COARSENED);
            }
        }
    }
}



void MeshRefinement::uniformly_refine (unsigned int n)
{
  // Refine n times
  // FIXME - this won't work if n>1 and the mesh 
  // has already been attached to an equation system
  for (unsigned int rstep=0; rstep<n; rstep++)
    {
      // Clean up the refinement flags
      this->clean_refinement_flags();
      
      // Flag all the active elements for refinement.       
      MeshBase::element_iterator       elem_it  = _mesh.active_elements_begin();
      const MeshBase::element_iterator elem_end = _mesh.active_elements_end(); 

      for ( ; elem_it != elem_end; ++elem_it)
	(*elem_it)->set_refinement_flag(Elem::REFINE);

      // Refine all the elements we just flagged.
      this->_refine_elements();
    }
  
  // Finally, the new mesh needs to be prepared for use
  _mesh.prepare_for_use ();
}



void MeshRefinement::uniformly_coarsen (unsigned int n)
{
  // Coarsen n times
  for (unsigned int rstep=0; rstep<n; rstep++)
    {
      // Clean up the refinement flags
      this->clean_refinement_flags();
      
      // Flag all the active elements for coarsening
      MeshBase::element_iterator       elem_it  = _mesh.active_elements_begin();
      const MeshBase::element_iterator elem_end = _mesh.active_elements_end(); 

      for ( ; elem_it != elem_end; ++elem_it)
        {
	  (*elem_it)->set_refinement_flag(Elem::COARSEN);
          if ((*elem_it)->parent()) 
	    (*elem_it)->parent()->set_refinement_flag(Elem::COARSEN_INACTIVE);
        }

      // Coarsen all the elements we just flagged.
      this->_coarsen_elements();
    }
    
  
  // Finally, the new mesh needs to be prepared for use
  _mesh.prepare_for_use ();
}


#endif