patch_recovery_error_estimator.c

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

// $Id: patch_recovery_error_estimator.C 2946 2008-07-30 22:09:25Z jwpeterson $

// 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 <algorithm> // for std::fill
#include <cmath>     // for std::sqrt std::pow std::abs


// Local Includes
#include "libmesh_common.h"
#include "patch_recovery_error_estimator.h"
#include "dof_map.h"
#include "fe.h"
#include "dense_matrix.h"
#include "dense_vector.h"
#include "error_vector.h"
#include "quadrature_gauss.h"
#include "libmesh_logging.h"
#include "elem.h"
#include "patch.h"
#include "quadrature_grid.h"
#include "system.h"
#include "mesh.h"
#include "tensor_value.h"


//-----------------------------------------------------------------
// PatchRecoveryErrorEstimator implementations
std::vector<Real> PatchRecoveryErrorEstimator::specpoly(const unsigned int dim,
							const Order order,
							const Point p,
							const unsigned int matsize)
{
  std::vector<Real> psi;
  psi.reserve(matsize);
  Real x = p(0), y=0., z=0.;
  if (dim > 1)
    y = p(1);
  if (dim > 2)
    z = p(2);
    
  // builds psi vector of form 1 x y z x^2 xy xz y^2 yz z^2 etc..
  // I haven't added 1D support here
  for (unsigned int poly_deg=0; poly_deg <= static_cast<unsigned int>(order) ; poly_deg++)
    { // loop over all polynomials of total degreee = poly_deg

      switch (dim)
	{
	  // 3D spectral polynomial basis functions
	case 3:
	  {	
	    for (int xexp=poly_deg; xexp >= 0; xexp--) // use an int for xexp since we -- it
	      for (int yexp=poly_deg-xexp; yexp >= 0; yexp--)
                {
                  int zexp = poly_deg - xexp - yexp;
		  psi.push_back(std::pow(x,static_cast<Real>(xexp))*
				std::pow(y,static_cast<Real>(yexp))*
				std::pow(z,static_cast<Real>(zexp)));
                }
	    break;
	  }

	  // 2D spectral polynomial basis functions
	case 2:
	  {
	    for (int xexp=poly_deg; xexp >= 0; xexp--) // use an int for xexp since we -- it
              {
                int yexp = poly_deg - xexp;
		psi.push_back(std::pow(x,static_cast<Real>(xexp))*
			      std::pow(y,static_cast<Real>(yexp)));
              }
	    break;
	  }

	  // 1D spectral polynomial basis functions
	case 1:
	  {
            int xexp = poly_deg;
	    psi.push_back(std::pow(x,static_cast<Real>(xexp)));
	    break;
	  }
	  
	default:
	  libmesh_error();
	}
    }

  return psi;
}
    
  

void PatchRecoveryErrorEstimator::estimate_error (const System& system,
						  ErrorVector& error_per_cell,
					          bool)
{
  START_LOG("estimate_error()", "PatchRecoveryErrorEstimator");

#ifdef ENABLE_SECOND_DERIVATIVES
  libmesh_assert (_sobolev_order == 1 || _sobolev_order == 2);
#else
  libmesh_assert (_sobolev_order == 1);
#endif

  // The current mesh
  const MeshBase& mesh = system.get_mesh();
  
  // The number of variables in the system
  const unsigned int n_vars = system.n_vars();

  // Resize the error_per_cell vector to be
  // the number of elements, initialize it to 0.
  error_per_cell.resize (mesh.max_elem_id());
  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);

  // Check for the use of component_mask
  this->convert_component_mask_to_scale();

  // Check for a valid component_scale
  if (!component_scale.empty())
    if (component_scale.size() != n_vars)
      {
	std::cerr << "ERROR: component_scale is the wrong size:"
		  << std::endl
		  << " component_scale.size()=" << component_scale.size()
		  << std::endl
		  << ", n_vars=" << n_vars
		  << std::endl;
	libmesh_error();
      }


  //------------------------------------------------------------
  // Iterate over all the active elements in the mesh
  // that live on this processor.
  Threads::parallel_for (ConstElemRange(mesh.active_local_elements_begin(),
					mesh.active_local_elements_end(),
					200),
			 EstimateError(system,
				       *this,
				       error_per_cell)
			 );

  // Each processor has now computed the error contribuions
  // for its local elements, and error_per_cell contains 0 for all the
  // non-local elements.  Summing the vector will provide the L_oo value
  // for each element, local or remote
  this->reduce_error(error_per_cell);
  
  STOP_LOG("estimate_error()", "PatchRecoveryErrorEstimator");
}



void PatchRecoveryErrorEstimator::EstimateError::operator()(const ConstElemRange &range) const
{
#ifdef ENABLE_SECOND_DERIVATIVES
  libmesh_assert (error_estimator._sobolev_order == 1 || error_estimator._sobolev_order == 2);
#else
  libmesh_assert (error_estimator._sobolev_order == 1);
#endif

  // The current mesh
  const MeshBase& mesh = system.get_mesh();

  // The dimensionality of the mesh
  const unsigned int dim = mesh.mesh_dimension();
  
  // The number of variables in the system
  const unsigned int n_vars = system.n_vars();
  
  // The DofMap for this system
  const DofMap& dof_map = system.get_dof_map();


  //------------------------------------------------------------
  // Iterate over all the elements in the range.
  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it!=range.end(); ++elem_it)
    {
      // elem is necessarily an active element on the local processor
      const Elem* elem = *elem_it;

      // Build a patch containing the current element
      // and its neighbors on the local processor
      Patch patch;

      // Use user specified patch size and growth strategy
      patch.build_around_element (elem, error_estimator.target_patch_size,
				  error_estimator.patch_growth_strategy);

      //------------------------------------------------------------
      // Process each variable in the system using the current patch
      for (unsigned int var=0; var<n_vars; var++)
	{
	  // Possibly skip this variable
	  if (!error_estimator.component_scale.empty())
	    if (error_estimator.component_scale[var] == 0.0) continue;
	  
	  // The type of finite element to use for this variable
	  const FEType& fe_type = dof_map.variable_type (var);

          const Order element_order  = static_cast<Order>
            (fe_type.order + elem->p_level());
      
	  // Finite element object for use in this patch
	  AutoPtr<FEBase> fe (FEBase::build (dim, fe_type));
	  
	  // Build an appropriate Gaussian quadrature rule
	  AutoPtr<QBase> qrule (fe_type.default_quadrature_rule(dim));

	  // Tell the finite element about the quadrature rule.
	  fe->attach_quadrature_rule (qrule.get());
      
	  // Get Jacobian values, etc..
	  const std::vector<Real>&                       JxW     = fe->get_JxW();
	  const std::vector<Point>&                      q_point = fe->get_xyz();
#ifndef NDEBUG
	  // We only use phi below to assert that the dof_indices
	  // vector has the correct size.  So we avoid declaring it