fe_base.c

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

                          new_side_dofs[sidej];
                        if (dof_is_fixed[j])
                          Fe(freei) -=
                            phi_coarse[i][qp] *
                            phi_coarse[j][qp] * JxW[qp] *
                            Ue(j);
                        else
                          Ke(freei,freej) +=
                            phi_coarse[i][qp] *
                            phi_coarse[j][qp] * JxW[qp];
                        if (cont == C_ONE)
                          {
                            if (dof_is_fixed[j])
                              Fe(freei) -=
                                ((*dphi_coarse)[i][qp] *
                                 (*dphi_coarse)[j][qp]) *
                                JxW[qp] *
                                Ue(j);
                            else
                              Ke(freei,freej) +=
                                ((*dphi_coarse)[i][qp] *
                                 (*dphi_coarse)[j][qp])
                                * JxW[qp];
                          }
                        if (!dof_is_fixed[j])
                          freej++;
                      }
                    Fe(freei) += phi_coarse[i][qp] *
                                 fineval * JxW[qp];
                    if (cont == C_ONE)
                      Fe(freei) +=
                        (finegrad * (*dphi_coarse)[i][qp]) * JxW[qp];
                    freei++;
                  }
              }
          }
        Ke.cholesky_solve(Fe, Uedge);

        // Transfer new edge solutions to element
        for (unsigned int i=0; i != free_dofs; ++i)
          {
            Number &ui = Ue(new_side_dofs[free_dof[i]]);
            libmesh_assert(std::abs(ui) < TOLERANCE ||
                   std::abs(ui - Uedge(i)) < TOLERANCE);
            ui = Uedge(i);
            dof_is_fixed[new_side_dofs[free_dof[i]]] =
              true;
          }
      }
   
  // Project any side values (edges in 2D, faces in 3D)
  if (dim > 1 && cont != DISCONTINUOUS)
    for (unsigned int s=0; s != elem->n_sides(); ++s)
      {
        FEInterface::dofs_on_side(elem, dim, fe_type,
                                  s, new_side_dofs);

        // Some side dofs are on nodes/edges 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 side coefficients
        DenseVector<Number> Uside(free_dofs);

        // Add projection terms from each child sharing
        // this side
        for (unsigned int c=0; c != elem->n_children();
             ++c)
          {
            if (!elem->is_child_on_side(c,s))
              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_side(child, dim,
              temp_fe_type, s, old_side_dofs);

            // Initialize both child and parent FE data
            // on the child's side
            fe->attach_quadrature_rule (qsiderule.get());
            fe->reinit (child, s);
            const unsigned int n_qp = qsiderule->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 side 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 =
                          new_side_dofs[sidej];
                        if (dof_is_fixed[j])
                          Fe(freei) -=
                            phi_coarse[i][qp] *
                            phi_coarse[j][qp] * JxW[qp] *
                            Ue(j);
                        else
                          Ke(freei,freej) +=
                            phi_coarse[i][qp] *
                            phi_coarse[j][qp] * JxW[qp];
                        if (cont == C_ONE)
                          {
                            if (dof_is_fixed[j])
                              Fe(freei) -=
                                ((*dphi_coarse)[i][qp] *
                                 (*dphi_coarse)[j][qp]) *
                                JxW[qp] *
                                Ue(j);
                            else
                              Ke(freei,freej) +=
                                ((*dphi_coarse)[i][qp] *
                                 (*dphi_coarse)[j][qp])
                                * JxW[qp];
                          }
                        if (!dof_is_fixed[j])
                          freej++;
                      }
                    Fe(freei) += (fineval * phi_coarse[i][qp]) * JxW[qp];
                    if (cont == C_ONE)
                      Fe(freei) +=
                        (finegrad * (*dphi_coarse)[i][qp]) * JxW[qp];
                    freei++;
                  }
              }
          }
        Ke.cholesky_solve(Fe, Uside);

        // Transfer new side solutions to element
        for (unsigned int i=0; i != free_dofs; ++i)
          {
            Number &ui = Ue(new_side_dofs[free_dof[i]]);
            libmesh_assert(std::abs(ui) < TOLERANCE ||
                   std::abs(ui - Uside(i)) < TOLERANCE);
            ui = Uside(i);
            dof_is_fixed[new_side_dofs[free_dof[i]]] =
              true;
          }
      }

  // Project the interior values, finally

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

  // Add projection terms from each child
  for (unsigned int c=0; c != elem->n_children(); ++c)
    {
      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();

      // Initialize both child and parent FE data
      // on the child's quadrature points
      fe->attach_quadrature_rule (qrule.get());
      fe->reinit (child);
      const unsigned int n_qp = qrule->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 interior projection matrix
          for (unsigned int i=0, freei=0;
               i != new_n_dofs; ++i)
            {
              // fixed DoFs aren't test functions
              if (dof_is_fixed[i])
                continue;
              for (unsigned int j=0, freej=0; j !=
                   new_n_dofs; ++j)
                {
                  if (dof_is_fixed[j])
                    Fe(freei) -=
                      phi_coarse[i][qp] *
                      phi_coarse[j][qp] * JxW[qp] *
                      Ue(j);
                  else
                    Ke(freei,freej) +=
                      phi_coarse[i][qp] *
                      phi_coarse[j][qp] * JxW[qp];
                  if (cont == C_ONE)
                    {
                      if (dof_is_fixed[j])
                        Fe(freei) -=
                          ((*dphi_coarse)[i][qp] *
                           (*dphi_coarse)[j][qp]) *
                          JxW[qp] * Ue(j);
                      else
                        Ke(freei,freej) +=
                          ((*dphi_coarse)[i][qp] *
                           (*dphi_coarse)[j][qp]) * JxW[qp];
                    }
                  if (!dof_is_fixed[j])
                    freej++;
                }
              Fe(freei) += phi_coarse[i][qp] * fineval *
                           JxW[qp];
              if (cont == C_ONE)
                Fe(freei) += (finegrad * (*dphi_coarse)[i][qp]) * JxW[qp];
              freei++;
            }
        }
    }
  Ke.cholesky_solve(Fe, Uint);

  // Transfer new interior solutions to element
  for (unsigned int i=0; i != free_dofs; ++i)
    {
      Number &ui = Ue(free_dof[i]);
      libmesh_assert(std::abs(ui) < TOLERANCE ||
             std::abs(ui - Uint(i)) < TOLERANCE);
      ui = Uint(i);
      dof_is_fixed[free_dof[i]] = true;
    }

  // Make sure every DoF got reached!
  for (unsigned int i=0; i != new_n_dofs; ++i)
    libmesh_assert(dof_is_fixed[i]);
}



void FEBase::compute_proj_constraints (DofConstraints &constraints,
				       DofMap &dof_map,
				       const unsigned int variable_number,
				       const Elem* elem)
{
  libmesh_assert (elem != NULL);

  const unsigned int Dim = elem->dim();

  // Only constrain elements in 2,3D.
  if (Dim == 1)
    return;

  // Only constrain active elements with this method
  if (!elem->active())
    return;

  const FEType& base_fe_type = dof_map.variable_type(variable_number);

  // Construct FE objects for this element and its neighbors.
  AutoPtr<FEBase> my_fe (FEBase::build(Dim, base_fe_type));
  const FEContinuity cont = my_fe->get_continuity();

  // We don't need to constrain discontinuous elements
  if (cont == DISCONTINUOUS)
    return;
  libmesh_assert (cont == C_ZERO || cont == C_ONE);

  AutoPtr<FEBase> neigh_fe (FEBase::build(Dim, base_fe_type));

  QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
  my_fe->attach_quadrature_rule (&my_qface);
  std::vector<Point> neigh_qface;

  const std::vector<Real>& JxW = my_fe->get_JxW();
  const std::vector<Point>& q_point = my_fe->get_xyz();
  const std::vector<std::vector<Real> >& phi = my_fe->get_phi();
  const std::vector<std::vector<Real> >& neigh_phi =
		  neigh_fe->get_phi();
  const std::vector<Point> *face_normals = NULL;
  const std::vector<std::vector<RealGradient> > *dphi = NULL;
  const std::vector<std::vector<RealGradient> > *neigh_dphi = NULL;

  std::vector<unsigned int> my_dof_indices, neigh_dof_indices;
  std::vector<unsigned int> my_side_dofs, neigh_side_dofs;

  if (cont != C_ZERO)
    {
      const std::vector<Point>& ref_face_normals =
        my_fe->get_normals();
      face_normals = &ref_face_normals;
      const std::vector<std::vector<RealGradient> >& ref_dphi =
	my_fe->get_dphi();
      dphi = &ref_dphi;
      const std::vector<std::vector<RealGradient> >& ref_neigh_dphi =
	neigh_fe->get_dphi();
      neigh_dphi = &ref_neigh_dphi;
    }

  DenseMatrix<Real> Ke;
  DenseVector<Real> Fe;
  std::vector<DenseVector<Real> > Ue;

  // Look at the element faces.  Check to see if we need to
  // build constraints.
  for (unsigned int s=0; s<elem->n_sides(); s++)
    if (elem->neighbor(s) != NULL)
      {
        // Get pointers to the element's neighbor.
        const Elem* neigh = elem->neighbor(s);

        // h refinement constraints:
        // constrain dofs shared between
        // this element and ones coarser
        // than this element.
        if (neigh->level() < elem->level()) 
          {
	    unsigned int s_neigh = neigh->which_neighbor_am_i(elem);
            libmesh_assert (s_neigh < neigh->n_neighbors());

            // Find the minimum p level; we build the h constraint
            // matrix with this and then constrain away all higher p
            // DoFs.
            libmesh_assert(neigh->active());
            const unsigned int min_p_level =
              std::min(elem->p_level(), neigh->p_level());

            // we may need to make the FE objects reinit with the
            // minimum shared p_level
            // FIXME - I hate using const_cast<> and avoiding
            // accessor functions; there's got to be a
            // better way to do this!
            const unsigned int old_elem_level = elem->p_level();
            if (old_elem_level != min_p_level)
              (const_cast<Elem *>(elem))->hack_p_level(min_p_level);
            const unsigned int old_neigh_level = neigh->p_level();
            if (old_neigh_level != min_p_level)
              (const_cast<Elem *>(neigh))->hack_p_level(min_p_level);

	    my_fe->reinit(elem, s);
	    
	    // This function gets called element-by-element, so there
	    // will be a lot of memory allocation going on.  We can 
	    // at least minimize this for the case of the dof indices
	    // by efficiently preallocating the requisite storage.
	    // n_nodes is not necessarily n_dofs, but it is better
	    // than nothing!
	    my_dof_indices.reserve    (elem->n_nodes());
	    neigh_dof_indices.reserve (neigh->n_nodes());