dof_map.c

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

// $Id: dof_map.C 2972 2008-08-11 19:54:46Z 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



// C++ Includes -------------------------------------
#include <set>
#include <algorithm> // for std::fill, std::equal_range, std::max, std::lower_bound, etc.

// Local Includes -----------------------------------
#include "coupling_matrix.h"
#include "dense_matrix.h"
#include "dense_vector_base.h"
#include "dof_map.h"
#include "elem.h"
#include "fe_interface.h"
#include "fe_type.h"
#include "fe_base.h" // FEBase::build() for continuity test
#include "libmesh_logging.h"
#include "mesh_base.h"
#include "mesh_tools.h"
#include "numeric_vector.h"
#include "parallel.h"
#include "sparse_matrix.h"
#include "string_to_enum.h"
#include "threads_allocators.h"


// ------------------------------------------------------------
// DofMap member functions

// Destructor
DofMap::~DofMap()
{
  this->clear();
}



// Add a copy of type to the vector.  Be sure to delete
// this memory in the destructor!
void DofMap::add_variable (const FEType& type)
{
  _variable_types.push_back (new FEType(type) );
}




Order DofMap::variable_order (const unsigned int c) const
{
  return _variable_types[c]->order;
}



void DofMap::attach_matrix (SparseMatrix<Number>& matrix)
{
  _matrices.push_back(&matrix);
  
  matrix.attach_dof_map (*this);
}



DofObject* DofMap::node_ptr(MeshBase& mesh, unsigned int i) const
{
  return mesh.node_ptr(i);
}



DofObject* DofMap::elem_ptr(MeshBase& mesh, unsigned int i) const
{
  return mesh.elem(i);
}



template <typename iterator_type>
void DofMap::set_nonlocal_dof_objects(iterator_type objects_begin,
                                      iterator_type objects_end,
                                      MeshBase &mesh,
                                      dofobject_accessor objects)
{
  // This function must be run on all processors at once
  parallel_only();

  // First, iterate over local objects to find out how many
  // are on each processor
  std::vector<unsigned int>
    ghost_objects_from_proc(libMesh::n_processors(), 0);

  iterator_type it  = objects_begin;

  for (; it != objects_end; ++it)
    {
      DofObject *obj = *it;

      if (obj)
        {
          unsigned int obj_procid = obj->processor_id();
          // We'd better be completely partitioned by now
          libmesh_assert(obj_procid != DofObject::invalid_processor_id);
          ghost_objects_from_proc[obj_procid]++;
        }
    }

  std::vector<unsigned int> objects_on_proc(libMesh::n_processors(), 0);
  Parallel::allgather(ghost_objects_from_proc[libMesh::processor_id()],
                      objects_on_proc);

#ifdef DEBUG
  for (unsigned int p=0; p != libMesh::n_processors(); ++p)
    libmesh_assert(ghost_objects_from_proc[p] <= objects_on_proc[p]);
#endif

  // Request sets to send to each processor
  std::vector<std::vector<unsigned int> >
    requested_ids(libMesh::n_processors());

  // We know how many of our objects live on each processor, so
  // reserve() space for requests from each.
  for (unsigned int p=0; p != libMesh::n_processors(); ++p)
    if (p != libMesh::processor_id())
      requested_ids[p].reserve(ghost_objects_from_proc[p]);

  for (it = objects_begin; it != objects_end; ++it)
    {
      DofObject *obj = *it;
      if (obj->processor_id() != DofObject::invalid_processor_id)
        requested_ids[obj->processor_id()].push_back(obj->id());
    }
#ifdef DEBUG
  for (unsigned int p=0; p != libMesh::n_processors(); ++p)
    libmesh_assert(requested_ids[p].size() == ghost_objects_from_proc[p]);
#endif

  // Next set ghost object n_comps from other processors
  for (unsigned int p=1; p != libMesh::n_processors(); ++p)
    {
      // Trade my requests with processor procup and procdown
      unsigned int procup = (libMesh::processor_id() + p) %
                             libMesh::n_processors();
      unsigned int procdown = (libMesh::n_processors() +
                               libMesh::processor_id() - p) %
                               libMesh::n_processors();
      std::vector<unsigned int> request_to_fill;
      Parallel::send_receive(procup, requested_ids[procup],
                             procdown, request_to_fill);

      // Fill those requests
      const unsigned int n_variables = this->n_variables();

      std::vector<unsigned int> ghost_data
        (request_to_fill.size() * 2 * n_variables);

      for (unsigned int i=0; i != request_to_fill.size(); ++i)
        {
          DofObject *requested = (this->*objects)(mesh, request_to_fill[i]);
          libmesh_assert(requested);
          libmesh_assert(requested->processor_id() == libMesh::processor_id());
          libmesh_assert(requested->n_vars(this->sys_number()) == n_variables);
          for (unsigned int v=0; v != n_variables; ++v)
            {
              unsigned int n_comp =
                requested->n_comp(this->sys_number(), v);
              ghost_data[i*2*n_variables+v] = n_comp;
              unsigned int first_dof = n_comp ?
                requested->dof_number(this->sys_number(), v, 0) : 0;
              libmesh_assert(first_dof != DofObject::invalid_id);
              ghost_data[i*2*n_variables+n_variables+v] = first_dof;
            }
        }

      // Trade back the results
      std::vector<unsigned int> filled_request;
      Parallel::send_receive(procdown, ghost_data,
                             procup, filled_request);

      // And copy the id changes we've now been informed of
      libmesh_assert(filled_request.size() ==
             requested_ids[procup].size() * 2 * n_variables);
      for (unsigned int i=0; i != requested_ids[procup].size(); ++i)
        {
          DofObject *requested = (this->*objects)(mesh, requested_ids[procup][i]);
          libmesh_assert(requested);
          libmesh_assert (requested->processor_id() == procup);
          for (unsigned int v=0; v != n_variables; ++v)
            {
              unsigned int n_comp = filled_request[i*2*n_variables+v];
              requested->set_n_comp(this->sys_number(), v, n_comp);
              if (n_comp)
                {
                  unsigned int first_dof = 
                    filled_request[i*2*n_variables+n_variables+v];
                  libmesh_assert(first_dof != DofObject::invalid_id);
                  requested->set_dof_number
                    (this->sys_number(), v, 0, first_dof);

		  // don't forget to add these remote dofs to the _send_list
		  for (unsigned int comp=0; comp!=n_comp; ++comp)
		    _send_list.push_back(first_dof+comp);
                }
            }
        }
    }

#ifdef DEBUG
  // Double check for invalid dofs
  for (it = objects_begin; it != objects_end; ++it)
    {
      DofObject *obj = *it;
      libmesh_assert (obj);
      unsigned int n_variables = obj->n_vars(this->sys_number());
      for (unsigned int v=0; v != n_variables; ++v)
        {
          unsigned int n_comp =
            obj->n_comp(this->sys_number(), v);
          unsigned int first_dof = n_comp ?
            obj->dof_number(this->sys_number(), v, 0) : 0;
          libmesh_assert(first_dof != DofObject::invalid_id);
        }
    }
#endif
}



void DofMap::reinit(MeshBase& mesh)
{
  libmesh_assert (mesh.is_prepared());
  
  START_LOG("reinit()", "DofMap");
  
  //this->clear();

  const unsigned int n_var = this->n_variables();
  const unsigned int dim   = mesh.mesh_dimension();

#ifdef ENABLE_AMR
  
  //------------------------------------------------------------
  // Clear the old_dof_objects for all the nodes
  // and elements so that we can overwrite them
  {
    MeshBase::node_iterator       node_it  = mesh.nodes_begin();
    const MeshBase::node_iterator node_end = mesh.nodes_end();
    
    for ( ; node_it != node_end; ++node_it)
      {
	(*node_it)->clear_old_dof_object();
	libmesh_assert ((*node_it)->old_dof_object == NULL);
      }
    
    MeshBase::element_iterator       elem_it  = mesh.elements_begin();
    const MeshBase::element_iterator elem_end = mesh.elements_end();
    
    for ( ; elem_it != elem_end; ++elem_it)
      {
	(*elem_it)->clear_old_dof_object();
	libmesh_assert ((*elem_it)->old_dof_object == NULL);
      }
  }

  
  //------------------------------------------------------------
  // Set the old_dof_objects for the elements that
  // weren't just created, if these old dof objects
  // had variables
  {
    MeshBase::element_iterator       elem_it  = mesh.elements_begin();
    const MeshBase::element_iterator elem_end = mesh.elements_end();
    
    for ( ; elem_it != elem_end; ++elem_it)
      {
	Elem* elem = *elem_it;

        // Skip the elements that were just refined
	if (elem->refinement_flag() == Elem::JUST_REFINED) continue;

	for (unsigned int n=0; n<elem->n_nodes(); n++)
	  {
	    Node* node = elem->get_node(n);

	    if (node->old_dof_object == NULL)
	      if (node->has_dofs())
		node->set_old_dof_object();
	  }

	libmesh_assert (elem->old_dof_object == NULL);

	if (elem->has_dofs())
	  elem->set_old_dof_object();
      }
  }

#endif // #ifdef ENABLE_AMR

  
  //------------------------------------------------------------
  // Then set the number of variables for each \p DofObject
  // equal to n_variables() for this system.  This will
  // handle new \p DofObjects that may have just been created
  {
    // All the nodes
    MeshBase::node_iterator       node_it  = mesh.nodes_begin();
    const MeshBase::node_iterator node_end = mesh.nodes_end();

    for ( ; node_it != node_end; ++node_it)
      (*node_it)->set_n_vars(this->sys_number(),n_var);

    // All the elements
    MeshBase::element_iterator       elem_it  = mesh.elements_begin();
    const MeshBase::element_iterator elem_end = mesh.elements_end();

    for ( ; elem_it != elem_end; ++elem_it)
      (*elem_it)->set_n_vars(this->sys_number(),n_var);
  }
  

  //------------------------------------------------------------
  // Next allocate space for the DOF indices
  for (unsigned int var=0; var<this->n_variables(); var++)
    {
      const FEType& base_fe_type = this->variable_type(var);
	     
      // This should be constant even on p-refined elements
      const bool extra_hanging_dofs =
	FEInterface::extra_hanging_dofs(base_fe_type);

      // For all the active elements
      MeshBase::element_iterator       elem_it  = mesh.active_elements_begin();
      const MeshBase::element_iterator elem_end = mesh.active_elements_end(); 
 
      // Count vertex degrees of freedom first
      for ( ; elem_it != elem_end; ++elem_it)
	{
	  Elem*    elem       = *elem_it;
	  const ElemType type = elem->type();

          FEType fe_type = base_fe_type;

#ifdef ENABLE_AMR
          // Make sure we haven't done more p refinement than we can
          // handle
          if (elem->p_level() + base_fe_type.order >
              FEInterface::max_order(base_fe_type, type))
            {
#  ifdef DEBUG
              if (FEInterface::max_order(base_fe_type,type) <
                  static_cast<unsigned int>(base_fe_type.order))
                {
                  std::cerr << "ERROR: Finite element "
                    << Utility::enum_to_string(base_fe_type.family)
                    << " on geometric element "
                    << Utility::enum_to_string(type) << std::endl
                    << "only supports FEInterface::max_order = " 
                    << FEInterface::max_order(base_fe_type,type)
                    << ", not fe_type.order = " << base_fe_type.order
                    << std::endl;
                }

              std::cerr << "WARNING: Finite element "
                    << Utility::enum_to_string(base_fe_type.family)
                    << " on geometric element "
                    << Utility::enum_to_string(type) << std::endl
                    << "could not be p refined past FEInterface::max_order = " 
                    << FEInterface::max_order(base_fe_type,type)
                    << std::endl;
#  endif
              elem->set_p_level(FEInterface::max_order(base_fe_type,type)
                                - base_fe_type.order);
            }
#endif

	  fe_type.order = static_cast<Order>(fe_type.order +
                                             elem->p_level());
	     
	  // Allocate the vertex DOFs
	  for (unsigned int n=0; n<elem->n_nodes(); n++)
	    {
	      Node* node = elem->get_node(n);

	      if (elem->is_vertex(n))