dof_map_constraints.c

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

  this->build_constraint_matrix (C, orig_col_dofs);

  START_LOG("constrain_elem_matrix()", "DofMap");
  
  row_dofs = orig_row_dofs;
  col_dofs = orig_col_dofs;
  
    
  // It is possible that the matrix is not constrained at all.
  if ((R.m() == matrix.m()) &&
      (R.n() == row_dofs.size()) &&
      (C.m() == matrix.n()) &&
      (C.n() == col_dofs.size())) // If the matrix is constrained
    {
      // K_constrained = R^T K C
      matrix.left_multiply_transpose  (R);
      matrix.right_multiply (C);
      
      
      libmesh_assert (matrix.m() == row_dofs.size());
      libmesh_assert (matrix.n() == col_dofs.size());
      
      
      for (unsigned int i=0; i<row_dofs.size(); i++)
	if (this->is_constrained_dof(row_dofs[i]))
	  {
	    for (unsigned int j=0; j<matrix.n(); j++)
	      matrix(i,j) = 0.;
	  
	    // If the DOF is constrained
	    matrix(i,i) = 1.;
	    
            if (asymmetric_constraint_rows)
              {
	        DofConstraints::const_iterator
	          pos = _dof_constraints.find(row_dofs[i]);
	    
	        libmesh_assert (pos != _dof_constraints.end());
	    
	        const DofConstraintRow& constraint_row = pos->second;
	    
	        libmesh_assert (!constraint_row.empty());
	    
	        for (DofConstraintRow::const_iterator
		       it=constraint_row.begin(); it != constraint_row.end();
		     ++it)
	          for (unsigned int j=0; j<col_dofs.size(); j++)
		    if (col_dofs[j] == it->first)
		      matrix(i,j) = -it->second;	
              }
	  }
    } // end if is constrained...
  
  STOP_LOG("constrain_elem_matrix()", "DofMap");  
}



void DofMap::constrain_element_vector (DenseVector<Number>&       rhs,
				       std::vector<unsigned int>& row_dofs,
				       bool) const
{
  libmesh_assert (rhs.size() == row_dofs.size());

  // check for easy return
  if (this->n_constrained_dofs() == 0)
    return;
  
  // The constrained RHS is built up as R^T F.  
  DenseMatrix<Number> R;

  this->build_constraint_matrix (R, row_dofs);

  START_LOG("constrain_elem_vector()", "DofMap");
    
  // It is possible that the vector is not constrained at all.
  if ((R.m() == rhs.size()) &&
      (R.n() == row_dofs.size())) // if the RHS is constrained
    {
      // Compute the matrix-vector product
      DenseVector<Number> old_rhs(rhs);

      // resize RHS & zero before summation
      rhs.resize(row_dofs.size());
      rhs.zero();

      // compute matrix/vector product
      for (unsigned int i=0; i<row_dofs.size(); i++)
	for (unsigned int j=0; j<old_rhs.size(); j++)
	  rhs(i) += R.transpose(i,j)*old_rhs(j);
      
      libmesh_assert (row_dofs.size() == rhs.size());

      for (unsigned int i=0; i<row_dofs.size(); i++)
	if (this->is_constrained_dof(row_dofs[i]))
	  {	
	    // If the DOF is constrained
	    rhs(i) = 0.;
	  }
    } // end if the RHS is constrained.
  
  STOP_LOG("constrain_elem_vector()", "DofMap");  
}



void DofMap::enforce_constraints_exactly (const System &system,
                                          NumericVector<Number> *v) const
{
  parallel_only();

  unsigned int local_constraints = this->n_constrained_dofs();
  Parallel::max(local_constraints);
  if (!local_constraints)
    return;

  if (!v)
    v = system.solution.get();

  NumericVector<Number> *v_local;
  NumericVector<Number> *v_global;
  AutoPtr<NumericVector<Number> > v_built;
  if (v->size() == v->local_size())
    {
      v_built = NumericVector<Number>::build();
      v_built->init(this->n_dofs(), this->n_local_dofs());
      v_built->close();

      for (unsigned int i=v_built->first_local_index();
           i<v_built->last_local_index(); i++)
        v_built->set(i, (*v)(i));
      v_built->close();
      v_global = v_built.get();

      v_local = v;
      libmesh_assert (v_local->closed());
    }
  else
    {
      v_built = NumericVector<Number>::build();
      v_built->init (v->size(), v->size());
      v->localize(*v_built);
      v_built->close();
      v_local = v_built.get();

      v_global = v;
    }

  const MeshBase &mesh = system.get_mesh();

  libmesh_assert (this == &(system.get_dof_map()));

  // indices on each element
  std::vector<unsigned int> local_dof_indices;

  MeshBase::const_element_iterator       elem_it  =
    mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator elem_end =
    mesh.active_local_elements_end(); 
      
  for ( ; elem_it != elem_end; ++elem_it)
    {
      const Elem* elem = *elem_it;

      this->dof_indices(elem, local_dof_indices);
      std::vector<unsigned int> raw_dof_indices = local_dof_indices;

      // Constraint matrix for each element
      DenseMatrix<Number> C;

      this->build_constraint_matrix (C, local_dof_indices);

      // Continue if the element is unconstrained
      if (!C.m())
        continue;

      libmesh_assert(C.m() == raw_dof_indices.size());
      libmesh_assert(C.n() == local_dof_indices.size());

      for (unsigned int i=0; i!=C.m(); ++i)
        {
          // Recalculate any constrained dof owned by this processor
          unsigned int global_dof = raw_dof_indices[i];
          if (this->is_constrained_dof(global_dof))
          {
            Number exact_value = 0;
            for (unsigned int j=0; j!=C.n(); ++j)
              {
                if (local_dof_indices[j] != global_dof)
                  exact_value += C(i,j) * 
                    (*v_local)(local_dof_indices[j]);
              }
            v_global->set(global_dof, exact_value);
          }
        }
    }

  // If the old vector was serial, we probably need to send our values
  // to other processors
  if (v->size() == v->local_size())
    {
      v_global->localize (*v);
    }
  v->close();
}



std::pair<Real, Real>
DofMap::max_constraint_error (const System &system,
                              NumericVector<Number> *v) const
{
  if (!v)
    v = system.solution.get();
  NumericVector<Number> &vec = *v;

  // We'll assume the vector is closed
  libmesh_assert (vec.closed());

  Real max_absolute_error = 0., max_relative_error = 0.;

  const MeshBase &mesh = system.get_mesh();

  libmesh_assert (this == &(system.get_dof_map()));

  // indices on each element
  std::vector<unsigned int> local_dof_indices;

  MeshBase::const_element_iterator       elem_it  =
    mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator elem_end =
    mesh.active_local_elements_end(); 
      
  for ( ; elem_it != elem_end; ++elem_it)
    {
      const Elem* elem = *elem_it;

      this->dof_indices(elem, local_dof_indices);
      std::vector<unsigned int> raw_dof_indices = local_dof_indices;

      // Constraint matrix for each element
      DenseMatrix<Number> C;

      this->build_constraint_matrix (C, local_dof_indices);

      // Continue if the element is unconstrained
      if (!C.m())
        continue;

      libmesh_assert(C.m() == raw_dof_indices.size());
      libmesh_assert(C.n() == local_dof_indices.size());

      for (unsigned int i=0; i!=C.m(); ++i)
        {
          // Recalculate any constrained dof owned by this processor
          unsigned int global_dof = raw_dof_indices[i];
          if (this->is_constrained_dof(global_dof) &&
              global_dof >= vec.first_local_index() &&
              global_dof < vec.last_local_index())
          {
            Number exact_value = 0;
            for (unsigned int j=0; j!=C.n(); ++j)
              {
                if (local_dof_indices[j] != global_dof)
                  exact_value += C(i,j) * 
                    vec(local_dof_indices[j]);
              }

            max_absolute_error = std::max(max_absolute_error,
              std::abs(vec(global_dof) - exact_value));
            max_relative_error = std::max(max_relative_error,
              std::abs(vec(global_dof) - exact_value)
              / std::abs(exact_value));
          }
        }
    }

  return std::pair<Real, Real>(max_absolute_error, max_relative_error);
}



void DofMap::build_constraint_matrix (DenseMatrix<Number>& C,
				      std::vector<unsigned int>& elem_dofs,
				      const bool called_recursively) const
{
  if (!called_recursively) START_LOG("build_constraint_matrix()", "DofMap");

  // Create a set containing the DOFs we already depend on
  typedef std::set<unsigned int> RCSet;
  RCSet dof_set;

  bool we_have_constraints = false;

  // Next insert any other dofs the current dofs might be constrained
  // in terms of.  Note that in this case we may not be done: Those
  // may in turn depend on others.  So, we need to repeat this process
  // in that case until the system depends only on unconstrained
  // degrees of freedom.
  for (unsigned int i=0; i<elem_dofs.size(); i++)
    if (this->is_constrained_dof(elem_dofs[i]))
      {
        we_have_constraints = true;

	// If the DOF is constrained
	DofConstraints::const_iterator
	  pos = _dof_constraints.find(elem_dofs[i]);
	
	libmesh_assert (pos != _dof_constraints.end());
	
	const DofConstraintRow& constraint_row = pos->second;
	
// Constraint rows in p refinement may be empty
//	libmesh_assert (!constraint_row.empty());
	
	for (DofConstraintRow::const_iterator
	       it=constraint_row.begin(); it != constraint_row.end();
	     ++it)
	  dof_set.insert (it->first);
      }

  // May be safe to return at this point
  // (but remember to stop the perflog)
  if (!we_have_constraints)
    {
      STOP_LOG("build_constraint_matrix()", "DofMap");
      return;
    }

  // delay inserting elem_dofs for efficiency in the case of
  // no constraints.  In that case we don't get here!
  dof_set.insert (elem_dofs.begin(),
		  elem_dofs.end());

  // If we added any DOFS then we need to do this recursively.
  // It is possible that we just added a DOF that is also
  // constrained!
  //
  // Also, we need to handle the special case of an element having DOFs
  // constrained in terms of other, local DOFs
  if ((dof_set.size() != elem_dofs.size()) || // case 1: constrained in terms of other DOFs
      !called_recursively)                    // case 2: constrained in terms of our own DOFs
    {
      // Create a new list of element DOFs containing the
      // contents of the current dof_set.
      std::vector<unsigned int> new_elem_dofs (dof_set.begin(),
					       dof_set.end());
      
      // Now we can build the constraint matrix.
      // Note that resize also zeros for a DenseMatrix<Number>.
      C.resize (elem_dofs.size(), new_elem_dofs.size());
      
      // Create the C constraint matrix.
      for (unsigned int i=0; i<elem_dofs.size(); i++)
	if (this->is_constrained_dof(elem_dofs[i]))
	  {
	    // If the DOF is constrained
	    DofConstraints::const_iterator
	      pos = _dof_constraints.find(elem_dofs[i]);
	    
	    libmesh_assert (pos != _dof_constraints.end());
	    
	    const DofConstraintRow& constraint_row = pos->second;
	    
// p refinement creates empty constraint rows
//	    libmesh_assert (!constraint_row.empty());
	    
	    for (DofConstraintRow::const_iterator
		   it=constraint_row.begin(); it != constraint_row.end();
		 ++it)
	      for (unsigned int j=0; j<new_elem_dofs.size(); j++)
		if (new_elem_dofs[j] == it->first)
		  C(i,j) = it->second;
	  }
	else
	  {
	    for (unsigned int j=0; j<new_elem_dofs.size(); j++)
	      if (new_elem_dofs[j] == elem_dofs[i])
		C(i,j) =  1.;
	  }	

      // May need to do this recursively.  It is possible
      // that we just replaced a constrained DOF with another
      // constrained DOF.
      elem_dofs = new_elem_dofs;
      
      DenseMatrix<Number> Cnew;
      
      this->build_constraint_matrix (Cnew, elem_dofs, true);

      if ((C.n() == Cnew.m()) &&
	  (Cnew.n() == elem_dofs.size())) // If the constraint matrix	                                 
	C.right_multiply(Cnew);           // is constrained...
      
      libmesh_assert (C.n() == elem_dofs.size());
    }
  
  if (!called_recursively) STOP_LOG("build_constraint_matrix()", "DofMap");  
}



void DofMap::allgather_recursive_constraints()
{
  // This function must be run on all processors at once
  parallel_only();

  // Return immediately if there's nothing to gather
  if (libMesh::n_processors() == 1)
    return;