hp_coarsentest.c

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

          const unsigned int n_dofs = dof_indices.size();

          // The number of nodes on the fine element
          const unsigned int n_nodes = elem->n_nodes();

	  // The average element value (used as an ugly hack
          // when we have nothing p-coarsened to compare to)
          // Real average_val = 0.;
          Number average_val = 0.;

	  // Calculate this variable's contribution to the p
	  // refinement error

	  if (elem->p_level() == 0)
	    {
              unsigned int n_vertices = 0;
              for (unsigned int n = 0; n != n_nodes; ++n)
                if (elem->is_vertex(n))
                  {
                    n_vertices++;
                    const Node * const node = elem->get_node(n);
                    average_val += system.current_solution
                      (node->dof_number(sys_num,var,0));
                  }
              average_val /= n_vertices;
	    }
          else
	    {
              unsigned int old_elem_level = elem->p_level();
              (const_cast<Elem *>(elem))->hack_p_level(old_elem_level - 1);
              
              fe_coarse->reinit(elem, &(qrule->get_points()));

              (const_cast<Elem *>(elem))->hack_p_level(old_elem_level);

              Ke.resize(phi_coarse->size(), phi_coarse->size());
              Ke.zero();
              Fe.resize(phi_coarse->size());
              Fe.zero();

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

                      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]));
                        }
                    }
                }

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

          // loop over the integration points on the fine element
          for (unsigned int qp=0; qp<n_qp; qp++)
            {
              Number value_error = 0.;
              Gradient grad_error;
              Tensor hessian_error;
              for (unsigned int i=0; i<n_dofs; i++)
                {
                  const unsigned int dof_num = dof_indices[i];
                  value_error += (*phi)[i][qp] *
                    system.current_solution(dof_num);
                  if (cont == C_ZERO || cont == C_ONE)
		    grad_error.add_scaled((*dphi)[i][qp], system.current_solution(dof_num));
		    // grad_error += (*dphi)[i][qp] *
                    //  system.current_solution(dof_num);
                  if (cont == C_ONE)
                    hessian_error.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
		  // hessian_error += (*d2phi)[i][qp] *
                  //    system.current_solution(dof_num);
                }
	      if (elem->p_level() == 0)
	        {
                  value_error -= average_val;
	        }
	      else
	        {
                  for (unsigned int i=0; i<Up.size(); i++)
                    {
                      value_error -= (*phi_coarse)[i][qp] * Up(i);
                      if (cont == C_ZERO || cont == C_ONE)
			grad_error.subtract_scaled((*dphi_coarse)[i][qp], Up(i));
                        // grad_error -= (*dphi_coarse)[i][qp] * Up(i);
                      if (cont == C_ONE)
			hessian_error.subtract_scaled((*d2phi_coarse)[i][qp], Up(i));
                        // hessian_error -= (*d2phi_coarse)[i][qp] * Up(i);
                    }
                }

	      p_error_per_cell[e_id] += component_scale[var] * 
		(*JxW)[qp] * libmesh_norm(value_error);
              if (cont == C_ZERO || cont == C_ONE)
	        p_error_per_cell[e_id] += component_scale[var] *
		  (*JxW)[qp] * grad_error.size_sq();
              if (cont == C_ONE)
	        p_error_per_cell[e_id] += component_scale[var] *
		  (*JxW)[qp] * hessian_error.size_sq();
            }

	  // Calculate this variable's contribution to the h
	  // refinement error

          if (!elem->parent())
            {
              // For now, we'll always start with an h refinement
              h_error_per_cell[e_id] =
                std::numeric_limits<ErrorVectorReal>::max() / 2;
            }
          else
            {
              FEInterface::inverse_map (dim, fe_type, coarse,
                *xyz_values, coarse_qpoints);

              unsigned int old_parent_level = coarse->p_level();
              (const_cast<Elem *>(coarse))->hack_p_level(elem->p_level());

              fe_coarse->reinit(coarse, &coarse_qpoints);

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

              // The number of DOFS on the coarse element
              unsigned int n_coarse_dofs = phi_coarse->size();

              // Loop over the quadrature points
              for (unsigned int qp=0; qp<n_qp; qp++)
                {
                  // The solution difference at the quadrature point
                  Number value_error = libMesh::zero;
                  Gradient grad_error;
                  Tensor hessian_error;

                  for (unsigned int i=0; i != n_dofs; ++i)
                    {
                      const unsigned int dof_num = dof_indices[i];
                      value_error += (*phi)[i][qp] *
                        system.current_solution(dof_num);
                      if (cont == C_ZERO || cont == C_ONE)
			grad_error.add_scaled((*dphi)[i][qp], system.current_solution(dof_num));
                        // grad_error += (*dphi)[i][qp] *
                        //  system.current_solution(dof_num);
                      if (cont == C_ONE)
			hessian_error.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
			// hessian_error += (*d2phi)[i][qp] *
                        //  system.current_solution(dof_num);
                    }

                  for (unsigned int i=0; i != n_coarse_dofs; ++i)
                    {
                      value_error -= (*phi_coarse)[i][qp] * Uc(i);
                      if (cont == C_ZERO || cont == C_ONE)
			// grad_error -= (*dphi_coarse)[i][qp] * Uc(i);
			grad_error.subtract_scaled((*dphi_coarse)[i][qp], Uc(i));
                      if (cont == C_ONE)
			hessian_error.subtract_scaled((*d2phi_coarse)[i][qp], Uc(i));
			// hessian_error -= (*d2phi_coarse)[i][qp] * Uc(i);
                    }

	          h_error_per_cell[e_id] += component_scale[var] * 
		    (*JxW)[qp] * libmesh_norm(value_error);
                  if (cont == C_ZERO || cont == C_ONE)
	            h_error_per_cell[e_id] += component_scale[var] * 
		      (*JxW)[qp] * grad_error.size_sq();
                  if (cont == C_ONE)
	            h_error_per_cell[e_id] += component_scale[var] * 
		      (*JxW)[qp] * hessian_error.size_sq();
                }

            }
	}
    }

  // Now that we've got our approximations for p_error and h_error, let's see
  // if we want to switch any h refinement flags to p refinement

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

      unsigned int dofs_per_elem = 0, dofs_per_p_elem = 0;

      // Loop over all the variables in the system
      for (unsigned int var=0; var<n_vars; var++)
        {
          // The type of finite element to use for this variable
          const FEType& fe_type = dof_map.variable_type (var);

          // FIXME: we're overestimating the number of DOFs added by h
          // refinement
          FEType elem_fe_type = fe_type;
          elem_fe_type.order =
            static_cast<Order>(fe_type.order + elem->p_level());
          dofs_per_elem += 
            FEInterface::n_dofs(dim, elem_fe_type, elem->type());

          elem_fe_type.order =
            static_cast<Order>(fe_type.order + elem->p_level() + 1);
          dofs_per_p_elem += 
            FEInterface::n_dofs(dim, elem_fe_type, elem->type());
        }

      const unsigned int new_h_dofs = dofs_per_elem *
        (elem->n_children() - 1);

      const unsigned int new_p_dofs = dofs_per_p_elem -
        dofs_per_elem;
      
std::cerr << "Cell " << e_id << ": h = " << elem->hmax()
          << ", p = " << elem->p_level() + 1 << "," << std::endl 
          << "     h_error = " << h_error_per_cell[e_id] 
          << ", p_error = " << p_error_per_cell[e_id] << std::endl
          << "     new_h_dofs = " << new_h_dofs
          << ", new_p_dofs = " << new_p_dofs << std::endl;
        
      if ((std::sqrt(p_error_per_cell[e_id]) * p_weight / new_p_dofs) 
          > (std::sqrt(h_error_per_cell[e_id]) / new_h_dofs))
/*
      if (std::sqrt(p_error_per_cell[e_id]) * p_weight
          > std::sqrt(h_error_per_cell[e_id]))
*/
        {
          elem->set_p_refinement_flag(Elem::REFINE);
//          elem->set_refinement_flag(Elem::COARSEN);
          elem->set_refinement_flag(Elem::DO_NOTHING);
        }
    }

  STOP_LOG("select_refinement()", "HPCoarsenTest");
}

#endif // #ifdef ENABLE_AMR