dof_map.c

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

	        {
	          const unsigned int old_node_dofs =
	            node->n_comp(this->sys_number(), var);

		  const unsigned int vertex_dofs =
		    std::max(FEInterface::n_dofs_at_node(dim, fe_type,
                                                         type, n),
                             old_node_dofs);
		  
		  // Some discontinuous FEs have no vertex dofs
		  if (vertex_dofs > old_node_dofs)
		    {
		      node->set_n_comp(this->sys_number(), var,
				       vertex_dofs);
		      // Abusing dof_number to set a "this is a
		      // vertex" flag
		      node->set_dof_number(this->sys_number(),
					   var, 0, vertex_dofs);
		    }
	        }
	    }
	}

      elem_it = mesh.active_elements_begin();

      for ( ; elem_it != elem_end; ++elem_it)
	{
	  Elem*    elem       = *elem_it;
	  const ElemType type = elem->type();
	     
          FEType fe_type = base_fe_type;
          fe_type.order = static_cast<Order>(fe_type.order +
                                             elem->p_level());

	  // Allocate the edge and face DOFs
	  for (unsigned int n=0; n<elem->n_nodes(); n++)
	    {
	      Node* node = elem->get_node(n);

	      const unsigned int old_node_dofs =
	        node->n_comp(this->sys_number(), var);

              const unsigned int vertex_dofs = old_node_dofs?
                node->dof_number (this->sys_number(),var,0):0;

	      const unsigned int new_node_dofs =
		FEInterface::n_dofs_at_node(dim, fe_type, type, n);

	      // We've already allocated vertex DOFs
	      if (elem->is_vertex(n))
	        {
		  libmesh_assert(old_node_dofs >= vertex_dofs);
		  libmesh_assert(vertex_dofs >= new_node_dofs);
		}
	      // We need to allocate the rest
	      else
	        {
		  // If this has no dofs yet, it needs no vertex
		  // dofs, so we just give it edge or face dofs
		  if (!old_node_dofs)
                    {
		      node->set_n_comp(this->sys_number(), var,
				       new_node_dofs);
		      // Abusing dof_number to set a "this has no
		      // vertex dofs" flag
                      if (new_node_dofs)
		        node->set_dof_number(this->sys_number(),
					     var, 0, 0);
                    }

		  // If this has dofs, but has no vertex dofs,
		  // it may still need more edge or face dofs if
                  // we're p-refined.
		  else if (vertex_dofs == 0)
		    {
                      if (new_node_dofs > old_node_dofs)
                        {
		          node->set_n_comp(this->sys_number(), var,
				           new_node_dofs);
		          node->set_dof_number(this->sys_number(),
					       var, 0, vertex_dofs);
                        }
		    }
		  // If this is another element's vertex,
		  // add more (non-overlapping) edge/face dofs if
		  // necessary
		  else if (extra_hanging_dofs)
		    {
                      if (new_node_dofs > old_node_dofs - vertex_dofs)
                        {
		          node->set_n_comp(this->sys_number(), var,
				           vertex_dofs + new_node_dofs);
		          node->set_dof_number(this->sys_number(),
					       var, 0, vertex_dofs);
                        }
		    }
		  // If this is another element's vertex, add any
		  // (overlapping) edge/face dofs if necessary
		  else
		    {
		      libmesh_assert(old_node_dofs >= vertex_dofs);
                      if (new_node_dofs > old_node_dofs)
                        {
		          node->set_n_comp(this->sys_number(), var,
				           new_node_dofs);
		          node->set_dof_number(this->sys_number(),
					       var, 0, vertex_dofs);
                        }
		    }
		}
	    }
          // Allocate the element DOFs
	  const unsigned int dofs_per_elem =
			  FEInterface::n_dofs_per_elem(dim, fe_type,
						       type);

	  elem->set_n_comp(this->sys_number(), var, dofs_per_elem);

	}
    }

  // Calling DofMap::reinit() by itself makes little sense,
  // so we won't bother with nonlocal DofObjects.
  // Those will be fixed by distribute_dofs
/*
  //------------------------------------------------------------
  // At this point, all n_comp and dof_number values on local
  // DofObjects should be correct, but a ParallelMesh might have
  // incorrect values on non-local DofObjects.  Let's request the
  // correct values from each other processor.

  this->set_nonlocal_n_comps(mesh.nodes_begin(),
                             mesh.nodes_end(),
                             mesh,
                             &DofMap::node_ptr);

  this->set_nonlocal_n_comps(mesh.elements_begin(),
                             mesh.elements_end(),
                             mesh,
                             &DofMap::elem_ptr);
*/

  //------------------------------------------------------------
  // Finally, clear all the current DOF indices
  // (distribute_dofs expects them cleared!)
  this->invalidate_dofs(mesh);
  
  STOP_LOG("reinit()", "DofMap");
}



void DofMap::invalidate_dofs(MeshBase& mesh) const
{
  const unsigned int sys_num = this->sys_number();

  // 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)->invalidate_dofs(sys_num);
  
  // All the 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)->invalidate_dofs(sys_num);
}



void DofMap::clear()
{
  // we don't want to clear
  // the coupling matrix!
  // It should not change...
  //_dof_coupling->clear();

  // Since we allocate memory when adding entries to the _variable_types
  // vector, we must free that memory here!
  for (unsigned int i=0; i<_variable_types.size(); ++i)
    delete _variable_types[i];

  
  _variable_types.clear();  
  _first_df.clear();
  _end_df.clear();
  _send_list.clear();
  _n_nz.clear();
  _n_oz.clear();


#ifdef ENABLE_AMR

  _dof_constraints.clear();
  _n_old_dfs = 0;

#endif

  _matrices.clear();

  _n_dfs = 0;
}



void DofMap::distribute_dofs (MeshBase& mesh)
{
  // This function must be run on all processors at once
  parallel_only();

  // Log how long it takes to distribute the degrees of freedom
  START_LOG("distribute_dofs()", "DofMap");

  libmesh_assert (mesh.is_prepared());

  const unsigned int proc_id = libMesh::processor_id();
  const unsigned int n_proc  = libMesh::n_processors();
  
  libmesh_assert (this->n_variables() > 0);
  libmesh_assert (proc_id < n_proc);
  
  // re-init in case the mesh has changed
  this->reinit(mesh);
  
  // By default distribute variables in a
  // var-major fashion, but allow run-time
  // specification
  bool node_major_dofs = libMesh::on_command_line ("--node_major_dofs");

  // The DOF counter, will be incremented as we encounter
  // new degrees of freedom
  unsigned int next_free_dof = 0;

  // Set temporary DOF indices on this processor
  if (node_major_dofs)
    this->distribute_local_dofs_node_major (next_free_dof, mesh, false);
  else
    this->distribute_local_dofs_var_major (next_free_dof, mesh, false);

  // Get DOF counts on all processors
  std::vector<unsigned int> dofs_on_proc(n_proc, 0);
  Parallel::allgather(next_free_dof, dofs_on_proc);

  // Resize the _first_df and _end_df arrays
  _first_df.resize(n_proc);
  _end_df.resize (n_proc);

  // Get DOF offsets
  _first_df[0] = 0;
  for (unsigned int i=1; i < n_proc; ++i)
    _first_df[i] = _end_df[i-1] = _first_df[i-1] + dofs_on_proc[i-1];
  _end_df[n_proc-1] = _first_df[n_proc-1] + dofs_on_proc[n_proc-1];

  // Clear all the current DOF indices
  // (distribute_dofs expects them cleared!)
  this->invalidate_dofs(mesh);

  next_free_dof = _first_df[proc_id];

  // Set permanent DOF indices on this processor
  if (node_major_dofs)
    this->distribute_local_dofs_node_major (next_free_dof, mesh, true);
  else
    this->distribute_local_dofs_var_major (next_free_dof, mesh, true);

  libmesh_assert(next_free_dof == _end_df[proc_id]);

  //------------------------------------------------------------
  // At this point, all n_comp and dof_number values on local
  // DofObjects should be correct, but a ParallelMesh might have
  // incorrect values on non-local DofObjects.  Let's request the
  // correct values from each other processor.

  if (libMesh::n_processors() > 1)
    {
      this->set_nonlocal_dof_objects(mesh.nodes_begin(),
                                     mesh.nodes_end(),
                                     mesh, &DofMap::node_ptr);

      this->set_nonlocal_dof_objects(mesh.elements_begin(),
                                     mesh.elements_end(),
                                     mesh, &DofMap::elem_ptr);
    }
  
  // Set the total number of degrees of freedom
#ifdef ENABLE_AMR
  _n_old_dfs = _n_dfs;
#endif
  _n_dfs = _end_df[n_proc-1];

  STOP_LOG("distribute_dofs()", "DofMap");

  // Note that in the add_neighbors_to_send_list nodes on processor
  // boundaries that are shared by multiple elements are added for
  // each element.
  this->add_neighbors_to_send_list(mesh);
  
  // Here we need to clean up that data structure
  this->sort_send_list ();
}


void DofMap::distribute_local_dofs_node_major(unsigned int &next_free_dof,
                                              MeshBase& mesh,
                                              const bool build_send_list)
{
  const unsigned int sys_num = this->sys_number();
  const unsigned int n_vars  = this->n_variables();

  // Number of elements to add to send_list
  unsigned int send_list_size = 0;
  
  //-------------------------------------------------------------------------
  // First count and assign temporary numbers to local dofs
  MeshBase::element_iterator       elem_it  = mesh.active_local_elements_begin();
  const MeshBase::element_iterator elem_end = mesh.active_local_elements_end();

  for ( ; elem_it != elem_end; ++elem_it)
    {
      // Only number dofs connected to active
      // elements on this processor.
      Elem* elem                 = *elem_it;
      const unsigned int n_nodes = elem->n_nodes();
      
      // First number the nodal DOFS
      for (unsigned int n=0; n<n_nodes; n++)
        {
          Node* node = elem->get_node(n);
              
          for (unsigned var=0; var<n_vars; var++)
            {
              // assign dof numbers (all at once) if this is
              // our node and if they aren't already there
              if ((node->n_comp(sys_num,var) > 0) &&
                  (node->processor_id() == libMesh::processor_id()) &&
                  (node->dof_number(sys_num,var,0) ==
                   DofObject::invalid_id))
                {
                  node->set_dof_number(sys_num,
                                       var,
                                       0,
                                       next_free_dof);
                  next_free_dof += node->n_comp(sys_num,var);

		  // If we've been requested to, add local dofs to the
		  // _send_list
                  if (build_send_list)
                    for (unsigned int index=0; index<node->n_comp(sys_num,var);
                         index++)
                      _send_list.push_back(node->dof_number(sys_num,
                                                            var,
                                                            index));
                  else
                    send_list_size += node->n_comp(sys_num,var);
                }
            }
        }
                  
      // Now number the element DOFS
      for (unsigned var=0; var<n_vars; var++)
        if (elem->n_comp(sys_num,var) > 0)
          {
            libmesh_assert (elem->dof_number(sys_num,var,0) ==
                    DofObject::invalid_id);

            elem->set_dof_number(sys_num,
                                 var,
                                 0,
                                 next_free_dof);

            next_free_dof += elem->n_comp(sys_num,var);
              
	    // If we've been requested to, add local dofs to the
	    // _send_list
            if (build_send_list)
              for (unsigned int index=0; index<elem->n_comp(sys_num,var);
                   index++)
                _send_list.push_back(elem->dof_number(sys_num,
                                                      var,
                                                      index));
            else
              send_list_size += elem->n_comp(sys_num,var);
          }
    }

#ifdef DEBUG
// Make sure we didn't miss any nodes
  MeshTools::libmesh_assert_valid_node_procids(mesh);

  MeshBase::node_iterator       node_it  = mesh.local_nodes_begin();
  const MeshBase::node_iterator node_end = mesh.local_nodes_end();
  for (; node_it != node_end; ++node_it)
    {
      Node *obj = *node_it;
      libmesh_assert(obj);
      unsigned int n_variables = obj->n_vars(this->sys_number());
      for (unsigned int v=0; v != n_variables; ++v)
        {