system_projection.c

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

// $Id: system_projection.C 2789 2008-04-13 02:24:40Z roystgnr $

// 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 <vector>

// Local includes
#include "system.h"
#include "mesh.h"
#include "elem.h"
#include "dof_map.h"
#include "fe_interface.h"
#include "numeric_vector.h"
#include "libmesh_logging.h"

#include "dense_matrix.h"
#include "fe_base.h"
#include "quadrature_gauss.h"
#include "dense_vector.h"
#include "threads.h"


// ------------------------------------------------------------
// System implementation
void System::project_vector (NumericVector<Number>& vector) const
{
  // Create a copy of the vector, which currently
  // contains the old data.
  AutoPtr<NumericVector<Number> >
    old_vector (vector.clone());

  // Project the old vector to the new vector
  this->project_vector (*old_vector, vector);
}


/**
 * This method projects the vector
 * via L2 projections or nodal
 * interpolations on each element.
 */
void System::project_vector (const NumericVector<Number>& old_v,
			     NumericVector<Number>& new_v) const
{
  START_LOG ("project_vector()", "System");

  /**
   * This method projects a solution from an old mesh to a current, refined
   * mesh.  The input vector \p old_v gives the solution on the
   * old mesh, while the \p new_v gives the solution (to be computed)
   * on the new mesh.
   */
  new_v.clear();

  // Resize the new vector and get a serial version.
  NumericVector<Number> *new_vector_ptr;
  AutoPtr<NumericVector<Number> > new_vector_built;
  NumericVector<Number> *local_old_vector;
  AutoPtr<NumericVector<Number> > local_old_vector_built;
  const NumericVector<Number> *old_vector_ptr;

  // If the old vector was uniprocessor, make the new
  // vector uniprocessor
  if (old_v.size() == old_v.local_size())
    {
      new_v.init (this->n_dofs(),
		  this->n_dofs());
      new_vector_ptr = &new_v;
      old_vector_ptr = &old_v;
    }

  // Otherwise it is a parallel, distributed vector, which
  // we need to localize.
  else
    {
      new_v.init (this->n_dofs(),
		  this->n_local_dofs());
      new_vector_built = NumericVector<Number>::build();
      local_old_vector_built = NumericVector<Number>::build();
      new_vector_ptr = new_vector_built.get();
      local_old_vector = local_old_vector_built.get();
      new_vector_ptr->init(this->n_dofs(), this->n_dofs());
      local_old_vector->init(old_v.size(), old_v.size());
      old_v.localize(*local_old_vector);
      local_old_vector->close();
      old_vector_ptr = local_old_vector;
    }
  // Note that the above will have zeroed the new_vector

  NumericVector<Number> &new_vector = *new_vector_ptr;
  const NumericVector<Number> &old_vector = *old_vector_ptr;
    
#ifdef ENABLE_AMR 

  Threads::parallel_for (ConstElemRange (this->get_mesh().active_local_elements_begin(),
					 this->get_mesh().active_local_elements_end(),
					 1000),
			 ProjectVector(*this,
				       old_vector,
				       new_vector)
			 );

  new_vector.close();

  // If the old vector was serial, we probably need to send our values
  // to other processors
  if (old_v.size() == old_v.local_size())
    {
      AutoPtr<NumericVector<Number> > dist_v = NumericVector<Number>::build();
      dist_v->init(this->n_dofs(), this->n_local_dofs());
      dist_v->close();
    
      for (unsigned int i=0; i!=dist_v->size(); i++)
        if (new_vector(i) != 0.0)
          dist_v->set(i, new_vector(i));

      dist_v->close();

      dist_v->localize (new_v);
      new_v.close();
    }
  // If the old vector was parallel, we need to update it
  // and free the localized copies
  else
    {
      // We may have to set dof values that this processor doesn't
      // own in certain special cases, like LAGRANGE FIRST or
      // HERMITE THIRD elements on second-order meshes
      for (unsigned int i=0; i!=new_v.size(); i++)
        if (new_vector(i) != 0.0)
          new_v.set(i, new_vector(i));
      new_v.close();
    }

  this->get_dof_map().enforce_constraints_exactly(*this, &new_v);

#else

  // AMR is disabled: simply copy the vector
  new_vector = old_vector;
  
#endif // #ifdef ENABLE_AMR

  STOP_LOG("project_vector()", "System");
}



/**
 * This method projects an analytic function onto the solution via L2
 * projections and nodal interpolations on each element.
 */
void System::project_solution (Number fptr(const Point& p,
                                           const Parameters& parameters,
                                           const std::string& sys_name,
                                           const std::string& unknown_name),
                               Gradient gptr(const Point& p,
                                             const Parameters& parameters,
                                             const std::string& sys_name,
                                             const std::string& unknown_name),
                               Parameters& parameters) const
{
  this->project_vector(fptr, gptr, parameters, *solution);

  solution->localize(*current_local_solution);
}



/**
 * This method projects an analytic function via L2 projections and
 * nodal interpolations on each element.
 */
void System::project_vector (Number fptr(const Point& p,
                                         const Parameters& parameters,
                                         const std::string& sys_name,
                                         const std::string& unknown_name),
                             Gradient gptr(const Point& p,
                                           const Parameters& parameters,
                                           const std::string& sys_name,
                                           const std::string& unknown_name),
                             Parameters& parameters,
                             NumericVector<Number>& new_vector) const
{
  START_LOG ("project_vector()", "System");

  Threads::parallel_for (ConstElemRange (this->get_mesh().active_local_elements_begin(),
					 this->get_mesh().active_local_elements_end(),
					 1000),
			 ProjectSolution(*this,
					 fptr,
					 gptr,
					 parameters,
					 new_vector)
			 );


  new_vector.close();

#ifdef ENABLE_AMR
  this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
#endif

  STOP_LOG("project_vector()", "System");
}


#ifndef ENABLE_AMR
void System::ProjectVector::operator()(const ConstElemRange &) const
{
  libmesh_error();
#else
void System::ProjectVector::operator()(const ConstElemRange &range) const
{
  // A vector for Lagrange element interpolation, indicating if we
  // have visited a DOF yet.  Note that this is thread-local storage,
  // hence shared DOFS that live on thread boundaries may be doubly
  // computed.  It is expected that this will be more efficient
  // than locking a thread-global version of already_done, though.
  //
  // FIXME: we should use this for non-Lagrange coarsening, too
  std::vector<bool> already_done (system.n_dofs(), false);

  // The number of variables in this system
  const unsigned int n_variables = system.n_vars();

  // The dimensionality of the current mesh
  const unsigned int dim = system.get_mesh().mesh_dimension();

  // The DofMap for this system
  const DofMap& dof_map = system.get_dof_map();

  // The element matrix and RHS for projections.
  // Note that Ke is always real-valued, whereas
  // Fe may be complex valued if complex number
  // support is enabled
  DenseMatrix<Real> Ke;
  DenseVector<Number> Fe;
  // The new element coefficients
  DenseVector<Number> Ue;


  // Loop over all the variables in the system
  for (unsigned int var=0; var<n_variables; var++)
    {
      // Get FE objects of the appropriate type
      const FEType& base_fe_type = dof_map.variable_type(var);     
      AutoPtr<FEBase> fe (FEBase::build(dim, base_fe_type));      
      AutoPtr<FEBase> fe_coarse (FEBase::build(dim, base_fe_type));      

      // Create FE objects with potentially different p_level
      FEType fe_type, temp_fe_type;

      // Prepare variables for non-Lagrange projection
      AutoPtr<QBase> qrule     (base_fe_type.default_quadrature_rule(dim));
      AutoPtr<QBase> qedgerule (base_fe_type.default_quadrature_rule(1));
      AutoPtr<QBase> qsiderule (base_fe_type.default_quadrature_rule(dim-1));
      std::vector<Point> coarse_qpoints;

      // The values of the shape functions at the quadrature
      // points
      const std::vector<std::vector<Real> >& phi_values =
	fe->get_phi();
      const std::vector<std::vector<Real> >& phi_coarse =
	fe_coarse->get_phi();

      // The gradients of the shape functions at the quadrature
      // points on the child element.
      const std::vector<std::vector<RealGradient> > *dphi_values =
        NULL;
      const std::vector<std::vector<RealGradient> > *dphi_coarse =
        NULL;

      const FEContinuity cont = fe->get_continuity();

      if (cont == C_ONE)
        {
          const std::vector<std::vector<RealGradient> >&
            ref_dphi_values = fe->get_dphi();
          dphi_values = &ref_dphi_values;
          const std::vector<std::vector<RealGradient> >&
            ref_dphi_coarse = fe_coarse->get_dphi();
          dphi_coarse = &ref_dphi_coarse;
        }

      // The Jacobian * quadrature weight at the quadrature points
      const std::vector<Real>& JxW =
	fe->get_JxW();
     
      // The XYZ locations of the quadrature points on the
      // child element
      const std::vector<Point>& xyz_values =
	fe->get_xyz();


      // The global DOF indices
      std::vector<unsigned int> new_dof_indices, old_dof_indices;
      // Side/edge DOF indices
      std::vector<unsigned int> new_side_dofs, old_side_dofs;
   
      // Iterate over the elements in the range
      for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
	{
	  const Elem* elem = *elem_it;
	  const Elem* parent = elem->parent();

          // Adjust the FE type for p-refined elements
          fe_type = base_fe_type;
          fe_type.order = static_cast<Order>(fe_type.order +
                                             elem->p_level());

          // We may need to remember the parent's p_level
          unsigned int old_parent_level = 0;

	  // Update the DOF indices for this element based on
          // the new mesh
	  dof_map.dof_indices (elem, new_dof_indices, var);

	  // The number of DOFs on the new element
	  const unsigned int new_n_dofs = new_dof_indices.size();

          // Fixed vs. free DoFs on edge/face projections
          std::vector<char> dof_is_fixed(new_n_dofs, false); // bools
          std::vector<int> free_dof(new_n_dofs, 0);

	  // The element type
	  const ElemType elem_type = elem->type();

	  // The number of nodes on the new element
	  const unsigned int n_nodes = elem->n_nodes();

          // Zero the interpolated values
          Ue.resize (new_n_dofs); Ue.zero();

	  // Update the DOF indices based on the old mesh.
	  // This is done in one of three ways:
	  // 1.) If the child was just refined then it was not
	  //     active in the previous mesh & hence has no solution
	  //     values on it.  In this case simply project or
	  //     interpolate the solution from the parent, who was
          //     active in the previous mesh
	  // 2.) If the child was just coarsened, obtaining a
	  //     well-defined solution may require doing independent
	  //     projections on nodes, edges, faces, and interiors
	  // 3.) If the child was active in the previous
	  //     mesh, we can just copy coefficients directly
	  if (elem->refinement_flag() == Elem::JUST_REFINED)
	    {
	      libmesh_assert (parent != NULL);
              old_parent_level = parent->p_level();

              // We may have done p refinement or coarsening as well;
              // if so then we need to reset the parent's p level
              // so we can get the right DoFs from it
              if (elem->p_refinement_flag() == Elem::JUST_REFINED)
                {
                  libmesh_assert(elem->p_level() > 0);
                  (const_cast<Elem *>(parent))->hack_p_level(elem->p_level() - 1);
                }
              else if (elem->p_refinement_flag() == Elem::JUST_COARSENED)
                {
                  (const_cast<Elem *>(parent))->hack_p_level(elem->p_level() + 1);
                }
	 
	      dof_map.old_dof_indices (parent, old_dof_indices, var);
            }
	  else
	    {
	      dof_map.old_dof_indices (elem, old_dof_indices, var);

              if (elem->p_refinement_flag() == Elem::DO_NOTHING)
	        libmesh_assert (old_dof_indices.size() == new_n_dofs);
              if (elem->p_refinement_flag() == Elem::JUST_COARSENED)
	        libmesh_assert (elem->has_children());
	    }

	  unsigned int old_n_dofs = old_dof_indices.size();

          if (fe_type.family != LAGRANGE) {