hp_coarsentest.c

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

// $Id: hp_coarsentest.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 <limits> // for std::numeric_limits::max
#include <math.h>    // for sqrt


// Local Includes
#include "hp_coarsentest.h"
#include "dense_matrix.h"
#include "dense_vector.h"
#include "dof_map.h"
#include "fe_base.h"
#include "fe_interface.h"
#include "libmesh_logging.h"
#include "elem.h"
#include "error_vector.h"
#include "mesh.h"
#include "quadrature.h"
#include "system.h"
#include "tensor_value.h"

#ifdef ENABLE_AMR

//-----------------------------------------------------------------
// HPCoarsenTest implementations

void HPCoarsenTest::add_projection(const System &system,
                                const Elem *elem,
                                unsigned int var)
{
  // If we have children, we need to add their projections instead
  if (!elem->active())
    {
      libmesh_assert(!elem->subactive());
      for (unsigned int c = 0; c != elem->n_children(); ++c)
        this->add_projection(system, elem->child(c), var);
      return;
    }

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

  // The type of finite element to use for this variable
  const FEType& fe_type = dof_map.variable_type (var);

  const FEContinuity cont = fe->get_continuity();

  fe->reinit(elem);

  dof_map.dof_indices(elem, dof_indices, var);

  FEInterface::inverse_map (system.get_mesh().mesh_dimension(),
    fe_type, coarse, *xyz_values, coarse_qpoints);

  fe_coarse->reinit(coarse, &coarse_qpoints);

  if (Uc.size() == 0)
    {
      Ke.resize(phi_coarse->size(), phi_coarse->size());
      Ke.zero();
      Fe.resize(phi_coarse->size());
      Fe.zero();
      Uc.resize(phi_coarse->size());
      Uc.zero();
    }
  libmesh_assert(Uc.size() == phi_coarse->size());

  // Loop over the quadrature points
  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
    {
      // The solution value at the quadrature point
      Number val = libMesh::zero;
      Gradient grad;
      Tensor hess;

      for (unsigned int i=0; i != dof_indices.size(); i++)
        {
          unsigned int dof_num = dof_indices[i];
          val += (*phi)[i][qp] *
            system.current_solution(dof_num);
          if (cont == C_ZERO || cont == C_ONE)
	    grad.add_scaled((*dphi)[i][qp],system.current_solution(dof_num));
           // grad += (*dphi)[i][qp] *
            //  system.current_solution(dof_num);
          if (cont == C_ONE)
	    hess.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
	    // hess += (*d2phi)[i][qp] *
            //  system.current_solution(dof_num);
        }

      // The projection matrix and vector
      for (unsigned int i=0; i != Fe.size(); ++i)
        {
          Fe(i) += (*JxW)[qp] * 
            (*phi_coarse)[i][qp]*val;
          if (cont == C_ZERO || cont == C_ONE)
            Fe(i) += (*JxW)[qp] *
              (grad*(*dphi_coarse)[i][qp]);
          if (cont == C_ONE)
	    Fe(i) += (*JxW)[qp] *
              hess.contract((*d2phi_coarse)[i][qp]);
	    // Fe(i) += (*JxW)[qp] *
            //  (*d2phi_coarse)[i][qp].contract(hess);

          for (unsigned int j=0; j != Fe.size(); ++j)
            {
              Ke(i,j) += (*JxW)[qp] *
                (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
              if (cont == C_ZERO || cont == C_ONE)
                Ke(i,j) += (*JxW)[qp] *
                  (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
              if (cont == C_ONE)
                Ke(i,j) += (*JxW)[qp] *
                  ((*d2phi_coarse)[i][qp].contract((*d2phi_coarse)[j][qp]));
            }
        }
    }
}

void HPCoarsenTest::select_refinement (System &system)
{
  START_LOG("select_refinement()", "HPCoarsenTest");

  // 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();

  // The system number (for doing bad hackery)
  const unsigned int sys_num = system.number();

  // 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();
	}
    }
  else
    {
      // No specified scaling.  Scale all variables by one.
      component_scale.resize (n_vars, 1.0);
    }

  // Resize the error_per_cell vectors to handle
  // the number of elements, initialize them to 0.
  std::vector<ErrorVectorReal> h_error_per_cell(mesh.n_elem(), 0.);
  std::vector<ErrorVectorReal> p_error_per_cell(mesh.n_elem(), 0.);
  
  // Loop over all the variables in the system
  for (unsigned int var=0; var<n_vars; var++)
    {
      // Possibly skip this variable
      if (!component_scale.empty())
	if (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);

      FEType low_p_fe_type = fe_type;
      FEType high_p_fe_type = fe_type;

      // Finite element objects for a fine (and probably a coarse)
      // element will be needed
      fe = FEBase::build (dim, fe_type);
      fe_coarse = FEBase::build (dim, fe_type);

      // Any cached coarse element results have expired
      coarse = NULL;
      unsigned int cached_coarse_p_level = 0;

      const FEContinuity cont = fe->get_continuity();
      libmesh_assert (cont == DISCONTINUOUS || cont == C_ZERO || 
	      cont == C_ONE);

      // Build an appropriate quadrature rule
      qrule = fe_type.default_quadrature_rule(dim);

      // Tell the refined finite element about the quadrature
      // rule.  The coarse finite element need not know about it
      fe->attach_quadrature_rule (qrule.get());

      // We will always do the integration
      // on the fine elements.  Get their Jacobian values, etc..
      JxW = &(fe->get_JxW());
      xyz_values = &(fe->get_xyz());

      // The shape functions
      phi = &(fe->get_phi());
      phi_coarse = &(fe_coarse->get_phi());

      // The shape function derivatives
      if (cont == C_ZERO || cont == C_ONE)
	{
	  dphi = &(fe->get_dphi());
	  dphi_coarse = &(fe_coarse->get_dphi());
	}

#ifdef ENABLE_SECOND_DERIVATIVES
      // The shape function second derivatives
      if (cont == C_ONE)
	{
	  d2phi = &(fe->get_d2phi());
	  d2phi_coarse = &(fe_coarse->get_d2phi());
	}
#endif // defined (ENABLE_SECOND_DERIVATIVES)

      // Iterate over all the active elements in the mesh
      // that live on this processor.

      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;

	  // We're only checking elements that are already flagged for h
          // refinement
          if (elem->refinement_flag() != Elem::REFINE)
            continue;

          const unsigned int e_id = elem->id();

          // Find the projection onto the parent element,
          // if necessary
          if (elem->parent() &&
              (coarse != elem->parent() ||
               cached_coarse_p_level != elem->p_level()))
            {
	      Uc.resize(0);

              coarse = elem->parent();
              cached_coarse_p_level = elem->p_level();

              unsigned int old_parent_level = coarse->p_level();
              (const_cast<Elem *>(coarse))->hack_p_level(elem->p_level());
              
              this->add_projection(system, coarse, var);

              (const_cast<Elem *>(coarse))->hack_p_level(old_parent_level);

              // Solve the h-coarsening projection problem
              Ke.cholesky_solve(Fe, Uc);
            }

	  fe->reinit(elem);
	  
          // Get the DOF indices for the fine element
          dof_map.dof_indices (elem, dof_indices, var);

          // The number of quadrature points
          const unsigned int n_qp = qrule->n_points();

          // The number of DOFS on the fine element