fe_base.c

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

	    (eta  >= 0.-eps) &&
	    (zeta >= 0.-eps) &&
	    ((xi + eta + zeta) <= 1.+eps))
	  return true;
		
	return false;
      }

      
    case HEX8:
    case HEX20:
    case HEX27:
      {
	/*
	  if ((xi   >= -1.) &&
	  (xi   <=  1.) &&
	  (eta  >= -1.) &&
	  (eta  <=  1.) &&
	  (zeta >= -1.) &&
	  (zeta <=  1.))
	  return true;
	*/
	
	// The reference hexahedral element is [-1,1]^3.
	if ((xi   >= -1.-eps) &&
	    (xi   <=  1.+eps) &&
	    (eta  >= -1.-eps) &&
	    (eta  <=  1.+eps) &&
	    (zeta >= -1.-eps) &&
	    (zeta <=  1.+eps))
	  {
	    //	    std::cout << "Strange Point:\n";
	    //	    p.print();
	    return true;
	  }

	return false;
      }

    case PRISM6:
    case PRISM15:
    case PRISM18:
      {
	// Figure this one out...
	// inside the reference triange with zeta in [-1,1]
	if ((xi   >=  0.-eps) &&
	    (eta  >=  0.-eps) &&
	    (zeta >= -1.-eps) &&
	    (zeta <=  1.+eps) &&
	    ((xi + eta) <= 1.+eps))
	  return true;

	return false;
      }


    case PYRAMID5:
      {
	std::cerr << "BEN: Implement this you lazy bastard!"
		  << std::endl;
	libmesh_error();

	return false;
      }

#ifdef ENABLE_INFINITE_ELEMENTS
    case INFHEX8:
      {      
	// The reference infhex8 is a [-1,1]^3.
	if ((xi   >= -1.-eps) &&
	    (xi   <=  1.+eps) &&
	    (eta  >= -1.-eps) &&
	    (eta  <=  1.+eps) &&
	    (zeta >= -1.-eps) &&
	    (zeta <=  1.+eps))
	  {
	    return true;
	  }
	return false;
      }

    case INFPRISM6:
      {      
	// inside the reference triange with zeta in [-1,1]
	if ((xi   >=  0.-eps) &&
	    (eta  >=  0.-eps) &&
	    (zeta >= -1.-eps) &&
	    (zeta <=  1.+eps) &&
	    ((xi + eta) <= 1.+eps))
	  {
	    return true;
	  }

	return false;
      }
#endif

    default:
      std::cerr << "ERROR: Unknown element type " << t << std::endl;
      libmesh_error();
    }

  // If we get here then the point is _not_ in the
  // reference element.   Better return false.
  
  return false;
}




void FEBase::print_JxW(std::ostream& os) const
{
  for (unsigned int i=0; i<JxW.size(); ++i)
    os << JxW[i] << std::endl;
}




void FEBase::print_phi(std::ostream& os) const
{
  for (unsigned int i=0; i<phi.size(); ++i)
    for (unsigned int j=0; j<phi[i].size(); ++j)
      os << " phi[" << i << "][" << j << "]=" << phi[i][j] << std::endl;
}




void FEBase::print_dphi(std::ostream& os) const
{
  for (unsigned int i=0; i<dphi.size(); ++i)
    for (unsigned int j=0; j<dphi[i].size(); ++j)
      os << " dphi[" << i << "][" << j << "]=" << dphi[i][j];
}



#ifdef ENABLE_SECOND_DERIVATIVES


void FEBase::print_d2phi(std::ostream& os) const
{
  for (unsigned int i=0; i<dphi.size(); ++i)
    for (unsigned int j=0; j<dphi[i].size(); ++j)
      os << " d2phi[" << i << "][" << j << "]=" << d2phi[i][j];
}

#endif




void FEBase::print_xyz(std::ostream& os) const
{
  for (unsigned int i=0; i<xyz.size(); ++i)
    os << xyz[i];
}




void FEBase::print_info(std::ostream& os) const
{
  os << "Shape functions at the Gauss pts." << std::endl;
  this->print_phi(os);
  
  os << "Shape function gradients at the Gauss pts." << std::endl;
  this->print_dphi(os);
  
  os << "XYZ locations of the Gauss pts." << std::endl;
  this->print_xyz(os);
  
  os << "Values of JxW at the Gauss pts." << std::endl;
  this->print_JxW(os);
}




std::ostream& operator << (std::ostream& os, const FEBase& fe)
{
  fe.print_info(os);
  return os;
}



#ifdef ENABLE_AMR

void FEBase::coarsened_dof_values(const NumericVector<Number> &old_vector,
                                  const DofMap &dof_map,
                                  const Elem *elem,
                                  DenseVector<Number> &Ue,
                                  const unsigned int var,
                                  const bool use_old_dof_indices)
{
  // Side/edge DOF indices
  std::vector<unsigned int> new_side_dofs, old_side_dofs;

  // FIXME: what about 2D shells in 3D space?
  unsigned int dim = elem->dim();

  // We use local FE objects for now
  // FIXME: we should use more, external objects instead for efficiency
  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));

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



  FEType fe_type = base_fe_type, temp_fe_type;
  const ElemType elem_type = elem->type();
  fe_type.order = static_cast<Order>(fe_type.order +
                                     elem->max_descendant_p_level());

  // Number of nodes on parent element
  const unsigned int n_nodes = elem->n_nodes();

  // Number of dofs on parent element
  const unsigned int new_n_dofs =
    FEInterface::n_dofs(dim, fe_type, elem_type);

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

  DenseMatrix<Real> Ke;
  DenseVector<Number> Fe;
  Ue.resize(new_n_dofs); Ue.zero();


  // When coarsening, in general, we need a series of
  // projections to ensure a unique and continuous
  // solution.  We start by interpolating nodes, then
  // hold those fixed and project edges, then
  // hold those fixed and project faces, then
  // hold those fixed and project interiors

  // Copy node values first
  {
  std::vector<unsigned int> node_dof_indices;
  if (use_old_dof_indices)
    dof_map.old_dof_indices (elem, node_dof_indices, var);
  else
    dof_map.dof_indices (elem, node_dof_indices, var);

  unsigned int current_dof = 0;
  for (unsigned int n=0; n!= n_nodes; ++n)
    {
      // FIXME: this should go through the DofMap,
      // not duplicate dof_indices code badly!
      const unsigned int my_nc =
        FEInterface::n_dofs_at_node (dim, fe_type,
                                     elem_type, n);
      if (!elem->is_vertex(n))
        {
          current_dof += my_nc;
          continue;
        }

      temp_fe_type = base_fe_type;
      // We're assuming here that child n shares vertex n,
      // which is wrong on non-simplices right now
      // ... but this code isn't necessary except on elements
      // where p refinement creates more vertex dofs; we have
      // no such elements yet.
/*
      if (elem->child(n)->p_level() < elem->p_level())
        {
          temp_fe_type.order = 
            static_cast<Order>(temp_fe_type.order +
                               elem->child(n)->p_level());
        }
*/
      const unsigned int nc =
        FEInterface::n_dofs_at_node (dim, temp_fe_type,
                                     elem_type, n);
      for (unsigned int i=0; i!= nc; ++i)
        {
          Ue(current_dof) =
            old_vector(node_dof_indices[current_dof]);
          dof_is_fixed[current_dof] = true;
          current_dof++;
        }
    }
  }

  // In 3D, project any edge values next
  if (dim > 2 && cont != DISCONTINUOUS)
    for (unsigned int e=0; e != elem->n_edges(); ++e)
      {
        FEInterface::dofs_on_edge(elem, dim, fe_type,
                                  e, new_side_dofs);

        // Some edge dofs are on nodes and already
        // fixed, others are free to calculate
        unsigned int free_dofs = 0;
        for (unsigned int i=0; i !=
             new_side_dofs.size(); ++i)
          if (!dof_is_fixed[new_side_dofs[i]])
            free_dof[free_dofs++] = i;
        Ke.resize (free_dofs, free_dofs); Ke.zero();
        Fe.resize (free_dofs); Fe.zero();
        // The new edge coefficients
        DenseVector<Number> Uedge(free_dofs);

        // Add projection terms from each child sharing
        // this edge
        for (unsigned int c=0; c != elem->n_children();
             ++c)
          {
            if (!elem->is_child_on_edge(c,e))
              continue;
            Elem *child = elem->child(c);

            std::vector<unsigned int> child_dof_indices;
            if (use_old_dof_indices)
              dof_map.old_dof_indices (child,
                child_dof_indices, var);
            else
              dof_map.dof_indices (child,
                child_dof_indices, var);
            const unsigned int child_n_dofs = child_dof_indices.size();

            temp_fe_type = base_fe_type;
            temp_fe_type.order = 
              static_cast<Order>(temp_fe_type.order +
                                 child->p_level());

            FEInterface::dofs_on_edge(child, dim,
              temp_fe_type, e, old_side_dofs);

            // Initialize both child and parent FE data
            // on the child's edge
            fe->attach_quadrature_rule (qedgerule.get());
            fe->edge_reinit (child, e);
            const unsigned int n_qp = qedgerule->n_points();

            FEInterface::inverse_map (dim, fe_type, elem,
                            xyz_values, coarse_qpoints);

            fe_coarse->reinit(elem, &coarse_qpoints);

            // Loop over the quadrature points
            for (unsigned int qp=0; qp<n_qp; qp++)
              {
                // solution value at the quadrature point
                Number fineval = libMesh::zero;
                // solution grad at the quadrature point
                Gradient finegrad;

                // Sum the solution values * the DOF
                // values at the quadrature point to
                // get the solution value and gradient.
                for (unsigned int i=0; i<child_n_dofs;
                     i++)
                  {
                    fineval +=
                      (old_vector(child_dof_indices[i])*
                      phi_values[i][qp]);
                    if (cont == C_ONE)
                      finegrad.add_scaled((*dphi_values)[i][qp],
                                          old_vector(child_dof_indices[i]));
                  }

                // Form edge projection matrix
                for (unsigned int sidei=0, freei=0; 
                     sidei != new_side_dofs.size();
                     ++sidei)
                  {
                    unsigned int i = new_side_dofs[sidei];
                    // fixed DoFs aren't test functions
                    if (dof_is_fixed[i])
                      continue;
                    for (unsigned int sidej=0, freej=0;
                         sidej != new_side_dofs.size();
                         ++sidej)
                      {
                        unsigned int j =