mesh_refinement.c

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

    (_mesh.elements_begin(), _mesh.elements_end(), psync);

  // If we weren't consistent in both h and p on every processor then
  // we weren't globally consistent
  bool parallel_consistent = hsync.parallel_consistent &&
                             psync.parallel_consistent;
  Parallel::min(parallel_consistent);

  STOP_LOG ("make_flags_parallel_consistent()", "MeshRefinement");

  return parallel_consistent;
}



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

  START_LOG ("make_coarsening_compatible()", "MeshRefinement");
  
  bool _maintain_level_one = maintain_level_one;

  // If the user used non-default parameters, let's warn that they're
  // deprecated
  if (!maintain_level_one)
    {
      deprecated();
    }
  else
    _maintain_level_one = _face_level_mismatch_limit;

  
  // Unless we encounter a specific situation level-one
  // will be satisfied after executing this loop just once   
  bool level_one_satisfied = true;

  
  // Unless we encounter a specific situation we will be compatible
  // with any selected refinement flags
  bool compatible_with_refinement = true;

  
  // find the maximum h and p levels in the mesh   
  unsigned int max_level = 0;
  unsigned int max_p_level = 0;
    
  
  // First we look at all the active level-0 elements.  Since it doesn't make
  // sense to coarsen them we must un-set their coarsen flags if
  // they are set.   
  MeshBase::element_iterator       el     = _mesh.active_elements_begin();
  const MeshBase::element_iterator end_el = _mesh.active_elements_end(); 

  for (; el != end_el; ++el)
    {      
      Elem *elem = *el;
      max_level = std::max(max_level, elem->level());
      max_p_level =
        std::max(max_p_level,
		 static_cast<unsigned int>(elem->p_level()));
      
      if ((elem->level() == 0) &&
	  (elem->refinement_flag() == Elem::COARSEN))
	elem->set_refinement_flag(Elem::DO_NOTHING);

      if ((elem->p_level() == 0) &&
	  (elem->p_refinement_flag() == Elem::COARSEN))
	elem->set_p_refinement_flag(Elem::DO_NOTHING);
    }
  
  // if there are no refined elements on this processor then
  // there is no work for us to do   
  if (max_level == 0 && max_p_level == 0)
    {
      STOP_LOG ("make_coarsening_compatible()", "MeshRefinement");

      // But we still have to check with other processors
      Parallel::min(compatible_with_refinement);

      return compatible_with_refinement;
    }

  
  
  // Loop over all the active elements.  If an element is marked
  // for coarsening we better check its neighbors.  If ANY of these neighbors
  // are marked for refinement AND are at the same level then there is a
  // conflict.  By convention refinement wins, so we un-mark the element for
  // coarsening.  Level-one would be violated in this case so we need to re-run
  // the loop.   
  if (_maintain_level_one)
    {
      
    repeat:
      level_one_satisfied = true;
  
      do
	{
	  level_one_satisfied = true;

	  MeshBase::element_iterator       el     = _mesh.active_elements_begin();
	  const MeshBase::element_iterator end_el = _mesh.active_elements_end(); 

	  for (; el != end_el; ++el)
	    {
	      Elem* elem = *el;
              bool my_flag_changed = false;
	      
	      if (elem->refinement_flag() == Elem::COARSEN) // If the element is active and 
		// the coarsen flag is set
		{
		  const unsigned int my_level = elem->level();
		  
		  for (unsigned int n=0; n<elem->n_neighbors(); n++)
		    if (elem->neighbor(n) != NULL &&      // I have a
		        elem->neighbor(n) != remote_elem) // neighbor here
		      if (elem->neighbor(n)->active()) // and it is active
			{
			  const Elem* neighbor = elem->neighbor(n);
			  
			  if ((neighbor->level() == my_level) &&
			      (neighbor->refinement_flag() == Elem::REFINE)) // the neighbor is at my level
	        	                                                     // and wants to be refined
			    {
			      elem->set_refinement_flag(Elem::DO_NOTHING);
                              my_flag_changed = true;
                              break;
			    }
			}
		      else // I have a neighbor and it is not active. That means it has children.
			{  // While it _may_ be possible to coarsen us if all the children of
			   // that element want to be coarsened, it is impossible to know at this
			   // stage.  Forget about it for the moment...  This can be handled in
			   // two steps.
			  elem->set_refinement_flag(Elem::DO_NOTHING);
			  my_flag_changed = true;
                          break;
			}
		}
	      if (elem->p_refinement_flag() == Elem::COARSEN) // If
                // the element is active and the order reduction flag is set
		{
		  const unsigned int my_p_level = elem->p_level();

		  for (unsigned int n=0; n<elem->n_neighbors(); n++)
		    if (elem->neighbor(n) != NULL &&      // I have a
		        elem->neighbor(n) != remote_elem) // neighbor here
		      if (elem->neighbor(n)->active()) // and it is active
			{
			  const Elem* neighbor = elem->neighbor(n);

                          if ((neighbor->p_level() > my_p_level &&
                               neighbor->p_refinement_flag() != Elem::COARSEN)
                              || (neighbor->p_level() == my_p_level &&
                               neighbor->p_refinement_flag() == Elem::REFINE))
			    {
			      elem->set_p_refinement_flag(Elem::DO_NOTHING);
			      my_flag_changed = true;
                              break;
			    }
			}
		      else // I have a neighbor and it is not active.
			{  // We need to find which of its children
                           // have me as a neighbor, and maintain
                           // level one p compatibility with them.
                           // Because we currently have level one h
                           // compatibility, we don't need to check
                           // grandchildren
			  const Elem* neighbor = elem->neighbor(n);

                           libmesh_assert(neighbor->has_children());
	                   for (unsigned int c=0; c!=neighbor->n_children(); c++)
                             {
                               Elem *subneighbor = neighbor->child(c);
                               if (subneighbor != remote_elem &&
                                   subneighbor->active() &&
                                   subneighbor->has_neighbor(elem))
                                 if ((subneighbor->p_level() > my_p_level &&
                                     subneighbor->p_refinement_flag() != Elem::COARSEN)
                                     || (subneighbor->p_level() == my_p_level &&
                                     subneighbor->p_refinement_flag() == Elem::REFINE))
			           {
			             elem->set_p_refinement_flag(Elem::DO_NOTHING);
			             my_flag_changed = true;
                                     break;
			           }
                             }
                           if (my_flag_changed)
                             break;
			}
		}
		  
	      // If the current element's flag changed, we hadn't
	      // satisfied the level one rule.
	      if (my_flag_changed)
                level_one_satisfied = false;
	      
	      // Additionally, if it has non-local neighbors, and
	      // we're not in serial, then we'll eventually have to
	      // return compatible_with_refinement = false, because
	      // our change has to propagate to neighboring
	      // processors.
	      if (my_flag_changed && !_mesh.is_serial())
              for (unsigned int n=0; n != elem->n_neighbors(); ++n)
                {
                  Elem *neigh = elem->neighbor(n);
                  if (!neigh)
                    continue;
                  if (neigh == remote_elem ||
                      neigh->processor_id() !=
                      libMesh::processor_id())
                    {
                      compatible_with_refinement = false;
                      break;
                    }
                  // FIXME - for non-level one meshes we should
                  // test all descendants
                  if (neigh->has_children())
                    for (unsigned int c=0; c != neigh->n_children(); ++c)
                      if (neigh->child(c) == remote_elem ||
                          neigh->child(c)->processor_id() !=
                          libMesh::processor_id())
                        {
                          compatible_with_refinement = false;
                          break;
                        }
                }
	    }      
	}
      while (!level_one_satisfied);
      
    } // end if (_maintain_level_one)
  
  
  // Next we look at all of the ancestor cells.
  // If there is a parent cell with all of its children
  // wanting to be unrefined then the element is a candidate
  // for unrefinement.  If all the children don't
  // all want to be unrefined then ALL of them need to have their
  // unrefinement flags cleared.    
  for (int level=(max_level); level >= 0; level--)
    {
      MeshBase::element_iterator       el     = _mesh.level_elements_begin(level);
      const MeshBase::element_iterator end_el = _mesh.level_elements_end(level); 
      for (; el != end_el; ++el)
        {
          Elem *elem = *el;
          if (elem->ancestor())
	    {
	  
	      // right now the element hasn't been disqualified
	      // as a candidate for unrefinement	   
	      bool is_a_candidate = true;
	      bool found_remote_child = false;
	      
	      for (unsigned int c=0; c<elem->n_children(); c++)
                {
                  Elem *child = elem->child(c);
                  if (child == remote_elem)
                    found_remote_child = true;
	          else if ((child->refinement_flag() != Elem::COARSEN) ||
		           !child->active() )
	            is_a_candidate = false;
                }
	      
	      if (!is_a_candidate && !found_remote_child)
	        {
	          elem->set_refinement_flag(Elem::INACTIVE);
		  
	          for (unsigned int c=0; c<elem->n_children(); c++)
                    {
                      Elem *child = elem->child(c);
                      if (child == remote_elem)
                        continue;
		      if (child->refinement_flag() == Elem::COARSEN)
		        {
		          level_one_satisfied = false;
		          child->set_refinement_flag(Elem::DO_NOTHING);
		        }
                    }
	        }
	    }
         }
     }
     
  if (!level_one_satisfied && _maintain_level_one) goto repeat;
  
  
  // If all the children of a parent are set to be coarsened
  // then flag the parent so that they can kill thier kids...   
  MeshBase::element_iterator       all_el     = _mesh.elements_begin();
  const MeshBase::element_iterator all_el_end = _mesh.elements_end(); 
  for (; all_el != all_el_end; ++all_el)
    {
      Elem *elem = *all_el;
        if (elem->ancestor())
          {
	
	    // Presume all the children are local and flagged for
	    // coarsening and then look for a contradiction	 
	    bool all_children_flagged_for_coarsening = true;
	    bool found_remote_child = false;
	
	    for (unsigned int c=0; c<elem->n_children(); c++)
              {
                Elem *child = elem->child(c);
	        if (child == remote_elem)
                  found_remote_child = true;
	        else if (child->refinement_flag() != Elem::COARSEN)
	          all_children_flagged_for_coarsening = false;
              }
	
	    if (!found_remote_child &&
                all_children_flagged_for_coarsening)
	      elem->set_refinement_flag(Elem::COARSEN_INACTIVE);
            else if (!found_remote_child)
	      elem->set_refinement_flag(Elem::INACTIVE);
          }
    }
	
  STOP_LOG ("make_coarsening_compatible()", "MeshRefinement");

  // If one processor finds an incompatibility, we're globally
  // incompatible
  Parallel::min(compatible_with_refinement);
  
  return compatible_with_refinement;
}








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

  START_LOG ("make_refinement_compatible()", "MeshRefinement");
  
  bool _maintain_level_one = maintain_level_one;

  // If the user used non-default parameters, let's warn that they're
  // deprecated
  if (!maintain_level_one)
    {
      deprecated();
    }
  else
    _maintain_level_one = _face_level_mismatch_limit;

  // Unless we encounter a specific situation level-one
  // will be satisfied after executing this loop just once
  bool level_one_satisfied = true;

  // Unless we encounter a specific situation we will be compatible
  // with any selected coarsening flags
  bool compatible_with_coarsening = true;

  // This loop enforces the level-1 rule.  We should only
  // execute it if the user indeed wants level-1 satisfied!   
  if (_maintain_level_one)
    {  
      do
	{
	  level_one_satisfied = true;
	  
	  MeshBase::element_iterator       el     = _mesh.active_elements_begin();
	  const MeshBase::element_iterator end_el = _mesh.active_elements_end(); 

	  for (; el != end_el; ++el)
            {
            Elem *elem = *el;
	    if (elem->refinement_flag() == Elem::REFINE)  // If the element is active and the
                                                          // h refinement flag is set    
	      {
		const unsigned int my_level = elem->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
		        neighbor->active()) // and it is active
		      {
			
			
			// Case 1:  The neighbor is at the same level I am.
			//        1a: The neighbor will be refined       -> NO PROBLEM
			//        1b: The neighbor won't be refined      -> NO PROBLEM
			//        1c: The neighbor wants to be coarsened -> PROBLEM			 
			if (neighbor->level() == my_level)
			  {
			    if (neighbor->refinement_flag() == Elem::COARSEN)
			      {
				neighbor->set_refinement_flag(Elem::DO_NOTHING);
                                if (neighbor->parent())
                                  neighbor->parent()->set_refinement_flag(Elem::INACTIVE);
				compatible_with_coarsening = false;
				level_one_satisfied = false;
			      }
			  }
			
			
			// Case 2: The neighbor is one level lower than I am.
			//         The neighbor thus MUST be refined to satisfy
			//         the level-one rule, regardless of whether it