system_projection.c

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

	    // For refined non-Lagrange elements, we do an L2
            // projection
            // FIXME: this will be a suboptimal and ill-defined
            // result if we're using non-nested finite element
            // spaces or if we're on a p-coarsened element!
	    if (elem->refinement_flag() == Elem::JUST_REFINED)
	      {
	        // Update the fe object based on the current child
                fe->attach_quadrature_rule (qrule.get());
	        fe->reinit (elem);
	  
	        // The number of quadrature points on the child
	        const unsigned int n_qp = qrule->n_points();

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

                fe_coarse->reinit(parent, &coarse_qpoints);

	        // Reinitialize the element matrix and vector for
	        // the current element.  Note that this will zero them
	        // before they are summed.
	        Ke.resize (new_n_dofs, new_n_dofs); Ke.zero();
	        Fe.resize (new_n_dofs); Fe.zero();
	  
	  
	        // Loop over the quadrature points
	        for (unsigned int qp=0; qp<n_qp; qp++)
	          {
	            // The solution value at the quadrature point	      
	            Number val = libMesh::zero;

		    // Sum the function values * the DOF values
		    // at the point of interest to get the function value
		    // (Note that the # of DOFs on the parent need not be the
		    //  same as on the child!)
		    for (unsigned int i=0; i<old_n_dofs; i++)
		      {
		        val += (old_vector(old_dof_indices[i])*
			        phi_coarse[i][qp]);
		      }

	            // Now \p val contains the solution value of variable
	            // \p var at the qp'th quadrature point on the child
	            // element.  It has been interpolated from the parent
	            // in case the child was just refined.  Next we will
	            // construct the L2-projection matrix for the element.

	            // Construct the Mass Matrix
	            for (unsigned int i=0; i<new_n_dofs; i++)
		      for (unsigned int j=0; j<new_n_dofs; j++)
		        Ke(i,j) += JxW[qp]*phi_values[i][qp]*phi_values[j][qp];

	            // Construct the RHS
	            for (unsigned int i=0; i<new_n_dofs; i++)
		      Fe(i) += JxW[qp]*phi_values[i][qp]*val;
	      
	          } // end qp loop

                Ke.cholesky_solve(Fe, Ue);

                // Fix up the parent's p level in case we changed it
                (const_cast<Elem *>(parent))->hack_p_level(old_parent_level);
	      }
            else if (elem->refinement_flag() == Elem::JUST_COARSENED)
	      {
                FEBase::coarsened_dof_values(old_vector, dof_map,
					     elem, Ue, var, true);
              }
	    // For unrefined uncoarsened elements, we just copy DoFs
	    else
	      {
                // FIXME - I'm sure this function would be about half
                // the size if anyone ever figures out how to improve
                // the DofMap interface... - RHS
                if (elem->p_refinement_flag() == Elem::JUST_REFINED)
                  {
                    libmesh_assert (elem->p_level() > 0);
                    // P refinement of non-hierarchic bases will
                    // require a whole separate code path
                    libmesh_assert (fe->is_hierarchic());
                    temp_fe_type = fe_type;
                    temp_fe_type.order =
                      static_cast<Order>(temp_fe_type.order - 1);
                    unsigned int old_index = 0, new_index = 0;
                    for (unsigned int n=0; n != n_nodes; ++n)
                      {
		        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(new_index + i) =
                              old_vector(old_dof_indices[old_index++]);
                          }
                        new_index +=
		          FEInterface::n_dofs_at_node (dim, fe_type,
                                                       elem_type, n);
                      }
                  }
                else if (elem->p_refinement_flag() ==
                         Elem::JUST_COARSENED)
                  {
                    // P coarsening of non-hierarchic bases will
                    // require a whole separate code path
                    libmesh_assert (fe->is_hierarchic());
                    temp_fe_type = fe_type;
                    temp_fe_type.order =
                      static_cast<Order>(temp_fe_type.order + 1);
                    unsigned int old_index = 0, new_index = 0;
                    for (unsigned int n=0; n != n_nodes; ++n)
                      {
		        const unsigned int nc =
		          FEInterface::n_dofs_at_node (dim, fe_type,
                                                       elem_type, n);
                        for (unsigned int i=0; i != nc; ++i)
                          {
                            Ue(new_index++) =
                              old_vector(old_dof_indices[old_index+i]);
                          }
                        old_index +=
		          FEInterface::n_dofs_at_node (dim, temp_fe_type,
                                                       elem_type, n);
                      }
                  }
                else
                  // If there's no p refinement, we can copy every DoF
		  for (unsigned int i=0; i<new_n_dofs; i++)
		    Ue(i) = old_vector(old_dof_indices[i]);
	      }
          } 
	  else { // fe type is Lagrange
	    // Loop over the DOFs on the element
	    for (unsigned int new_local_dof=0;
	         new_local_dof<new_n_dofs; new_local_dof++)
	      {
	        // The global DOF number for this local DOF
	        const unsigned int new_global_dof =
		  new_dof_indices[new_local_dof];

	        // The global DOF might lie outside of the bounds of a
	        // distributed vector.  Check for that and possibly continue
	        // the loop
	        if ((new_global_dof <  new_vector.first_local_index()) ||
		    (new_global_dof >= new_vector.last_local_index()))
		  continue;

	        // We might have already computed the solution for this DOF.
	        // This is likely in the case of a shared node, particularly
	        // at the corners of an element.  Check to see if that is the
	        // case
	        if (already_done[new_global_dof] == true)
		  continue;

		already_done[new_global_dof] = true;
	      
	        if (elem->refinement_flag() == Elem::JUST_REFINED)
		  {
		    // The location of the child's node on the parent element
		    const Point point =
		      FEInterface::inverse_map (dim, fe_type, parent,
					        elem->point(new_local_dof));
		  
		    // Sum the function values * the DOF values
		    // at the point of interest to get the function value
		    // (Note that the # of DOFs on the parent need not be the
		    //  same as on the child!)
		    for (unsigned int old_local_dof=0;
		         old_local_dof<old_n_dofs; old_local_dof++)
		      {
		        const unsigned int old_global_dof =
			  old_dof_indices[old_local_dof];
		      
		        Ue(new_local_dof) += 
                          (old_vector(old_global_dof)*
			  FEInterface::shape(dim, fe_type, parent,
					     old_local_dof, point));
		      }
		  }
	        else
		  {
		    // Get the old global DOF index
		    const unsigned int old_global_dof =
		      old_dof_indices[new_local_dof];

		    Ue(new_local_dof) = old_vector(old_global_dof);
		  }
              } // end local DOF loop

            // We may have to clean up a parent's p_level
	    if (elem->refinement_flag() == Elem::JUST_REFINED)
              (const_cast<Elem *>(parent))->hack_p_level(old_parent_level);
          }  // end fe_type if()

	  // Lock the new_vector since it is shared among threads.
	  {
	    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);

	    for (unsigned int i = 0; i < new_n_dofs; i++) 
	      if (Ue(i) != 0.)
		new_vector.set(new_dof_indices[i], Ue(i));
	  }
        }  // end elem loop
    } // end variables loop
#endif // ENABLE_AMR
}




void System::ProjectSolution::operator()(const ConstElemRange &range) const
{
  // We need data to project
  libmesh_assert(fptr);

  /**
   * This method projects an analytic solution to the current
   * mesh.  The input function \p fptr gives the analytic solution,
   * while the \p new_vector (which should already be correctly sized)
   * gives the solution (to be computed) on the current mesh.
   */

  // 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& fe_type = dof_map.variable_type(var);     
      AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));      

      // Prepare variables for projection
      AutoPtr<QBase> qrule     (fe_type.default_quadrature_rule(dim));
      AutoPtr<QBase> qedgerule (fe_type.default_quadrature_rule(1));
      AutoPtr<QBase> qsiderule (fe_type.default_quadrature_rule(dim-1));

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

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

      const FEContinuity cont = fe->get_continuity();

      if (cont == C_ONE)
        {
          // We'll need gradient data for a C1 projection
          libmesh_assert(gptr);

          const std::vector<std::vector<RealGradient> >&
            ref_dphi = fe->get_dphi();
          dphi = &ref_dphi;
        }

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

      // The global DOF indices
      std::vector<unsigned int> dof_indices;
      // Side/edge DOF indices
      std::vector<unsigned int> side_dofs;
   
      // Iterate over all the elements in the range
      for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
	{
	  const Elem* elem = *elem_it;

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

	  // The number of DOFs on the element
	  const unsigned int n_dofs = dof_indices.size();

          // Fixed vs. free DoFs on edge/face projections
          std::vector<char> dof_is_fixed(n_dofs, false); // bools
          std::vector<int> free_dof(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 (n_dofs); Ue.zero();

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

          // Interpolate node values first
          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 nc =
		FEInterface::n_dofs_at_node (dim, fe_type, elem_type,
                                             n);
              if (!elem->is_vertex(n))
                {
                  current_dof += nc;
                  continue;
                }
              if (cont == DISCONTINUOUS)
                {
                  libmesh_assert(nc == 0);
                }
              // Assume that C_ZERO elements have a single nodal
              // value shape function
              else if (cont == C_ZERO)
                {
                  libmesh_assert(nc == 1);
		  Ue(current_dof) = fptr(elem->point(n),
                                         parameters,
                                         system.name(),
                                         system.variable_name(var));
                  dof_is_fixed[current_dof] = true;
                  current_dof++;
                }
              // The hermite element vertex shape functions are weird
              else if (fe_type.family == HERMITE)
                {
                  Ue(current_dof) = fptr(elem->point(n),
                                         parameters,
                                         system.name(),
                                         system.variable_name(var));
                  dof_is_fixed[current_dof] = true;
                  current_dof++;
                  Gradient g = gptr(elem->point(n),
                                    parameters,
                                    system.name(),
                                    system.variable_name(var));
                  // x derivative
                  Ue(current_dof) = g(0);
                  dof_is_fixed[current_dof] = true;
                  current_dof++;
                  if (dim > 1)
                    {
                      // We'll finite difference mixed derivatives
                      Point nxminus = elem->point(n),
                            nxplus = elem->point(n);
                      nxminus(0) -= TOLERANCE;
                      nxplus(0) += TOLERANCE;
                      Gradient gxminus = gptr(nxminus,
                                              parameters,
                                              system.name(),
                                              system.variable_name(var));
                      Gradient gxplus = gptr(nxminus,
                                             parameters,
                                             system.name(),
                                             system.variable_name(var));
                      // y derivative
                      Ue(current_dof) = g(1);
                      dof_is_fixed[current_dof] = true;
                      current_dof++;
                      // xy derivative
                      Ue(current_dof) = (gxplus(1) - gxminus(1))
                                        / 2. / TOLERANCE;
                      dof_is_fixed[current_dof] = true;
                      current_dof++;

                      if (dim > 2)
                        {
                          // z derivative
                          Ue(current_dof) = g(2);
                          dof_is_fixed[current_dof] = true;
                          current_dof++;