C++ DiffContext::get_elem_solution_derivative方法代码示例

bool L2System::element_time_derivative (bool request_jacobian,
                                        DiffContext & context)
{
  FEMContext & c = cast_ref<FEMContext &>(context);

  // First we get some references to cell-specific data that
  // will be used to assemble the linear system.

  // Element Jacobian * quadrature weights for interior integration
  const std::vector<Real> & JxW = c.get_element_fe(0)->get_JxW();

  const std::vector<std::vector<Real>> & phi = c.get_element_fe(0)->get_phi();

  const std::vector<Point> & xyz = c.get_element_fe(0)->get_xyz();

  // The number of local degrees of freedom in each variable
  const unsigned int n_u_dofs = c.n_dof_indices(0);

  // The subvectors and submatrices we need to fill:
  DenseSubMatrix<Number> & K = c.get_elem_jacobian(0, 0);
  DenseSubVector<Number> & F = c.get_elem_residual(0);

  unsigned int n_qpoints = c.get_element_qrule().n_points();

  libmesh_assert (input_contexts.find(&c) != input_contexts.end());

  FEMContext & input_c = *input_contexts[&c];
  input_c.pre_fe_reinit(*input_system, &c.get_elem());
  input_c.elem_fe_reinit();

  for (unsigned int qp=0; qp != n_qpoints; qp++)
    {
      Number u = c.interior_value(0, qp);

      Number ufunc = (*goal_func)(input_c, xyz[qp]);

      for (unsigned int i=0; i != n_u_dofs; i++)
        F(i) += JxW[qp] * ((u - ufunc) * phi[i][qp]);
      if (request_jacobian)
        {
          const Number JxWxD = JxW[qp] *
            context.get_elem_solution_derivative();

          for (unsigned int i=0; i != n_u_dofs; i++)
            for (unsigned int j=0; j != n_u_dofs; ++j)
              K(i,j) += JxWxD * (phi[i][qp] * phi[j][qp]);
        }
    } // end of the quadrature point qp-loop

  return request_jacobian;
}

//for non-Dirichlet boundary conditions and the bit from diffusion term
bool ConvDiff_AuxSadjSys::side_time_derivative(bool request_jacobian, DiffContext & context)
{
    const unsigned int dim = this->get_mesh().mesh_dimension();

    FEMContext &ctxt = cast_ref<FEMContext&>(context);

    // First we get some references to cell-specific data that
    // will be used to assemble the linear system.
    FEBase* side_fe = NULL;
    ctxt.get_side_fe(aux_c_var, side_fe );

    // Element Jacobian * quadrature weights for interior integration
    const std::vector<Real> &JxW = side_fe->get_JxW();

    // Side basis functions
    const std::vector<std::vector<Real> > &phi = side_fe->get_phi();

    // Side Quadrature points
    const std::vector<Point > &qside_point = side_fe->get_xyz();

    //normal vector
    const std::vector<Point> &face_normals = side_fe->get_normals();

    // The number of local degrees of freedom in each variable
    const unsigned int n_c_dofs = ctxt.get_dof_indices(aux_c_var).size();

    // The subvectors and submatrices we need to fill:
    DenseSubMatrix<Number> &J_c_auxz = ctxt.get_elem_jacobian(aux_c_var, aux_zc_var);

    DenseSubVector<Number> &Rc = ctxt.get_elem_residual( aux_c_var );
    //Rf gets no contribution from sides

    unsigned int n_qpoints = ctxt.get_side_qrule().n_points();

    bool isEast = false;
    if (dim == 3) {
        isEast = ctxt.has_side_boundary_id(2);
    }
    else if (dim == 2) {
        isEast = ctxt.has_side_boundary_id(1);
    }

    //set (in)flux boundary condition on west side
    //homogeneous neumann (Danckwerts) outflow boundary condition on east side
    //no-flux (equivalently, homoegenous neumann) boundary conditions on north, south, top, bottom sides
    //"strong" enforcement of boundary conditions
    for (unsigned int qp=0; qp != n_qpoints; qp++)
    {
        Number aux_z = ctxt.side_value(aux_zc_var, qp);

        //velocity vector
        NumberVectorValue U(porosity*vx, 0., 0.);

        for (unsigned int i=0; i != n_c_dofs; i++)
        {
            if(isEast)
                Rc(i) += JxW[qp]*(-U*face_normals[qp]*aux_z)*phi[i][qp];

            if(request_jacobian && context.get_elem_solution_derivative())
            {
                for (unsigned int j=0; j != n_c_dofs; j++)
                {
                    if(isEast)
                        J_c_auxz(i,j) += JxW[qp]*(-U*face_normals[qp]*phi[j][qp])*phi[i][qp];
                }
            } // end - if (request_jacobian && context.get_elem_solution_derivative())
        } //end of outer dof (i) loop
    }

    return request_jacobian;
}

bool HeatSystem::element_time_derivative (bool request_jacobian,
                                          DiffContext & context)
{
  FEMContext & c = libmesh_cast_ref<FEMContext &>(context);

  const unsigned int mesh_dim =
    c.get_system().get_mesh().mesh_dimension();

  // First we get some references to cell-specific data that
  // will be used to assemble the linear system.
  const unsigned int dim = c.get_elem().dim();
  FEBase * fe = libmesh_nullptr;
  c.get_element_fe(T_var, fe, dim);

  // Element Jacobian * quadrature weights for interior integration
  const std::vector<Real> & JxW = fe->get_JxW();

  const std::vector<Point> & xyz = fe->get_xyz();

  const std::vector<std::vector<Real> > & phi = fe->get_phi();

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

  // The number of local degrees of freedom in each variable
  const unsigned int n_T_dofs = c.get_dof_indices(T_var).size();

  // The subvectors and submatrices we need to fill:
  DenseSubMatrix<Number> & K = c.get_elem_jacobian(T_var, T_var);
  DenseSubVector<Number> & F = c.get_elem_residual(T_var);

  // Now we will build the element Jacobian and residual.
  // Constructing the residual requires the solution and its
  // gradient from the previous timestep.  This must be
  // calculated at each quadrature point by summing the
  // solution degree-of-freedom values by the appropriate
  // weight functions.
  unsigned int n_qpoints = c.get_element_qrule().n_points();

  for (unsigned int qp=0; qp != n_qpoints; qp++)
    {
      // Compute the solution gradient at the Newton iterate
      Gradient grad_T = c.interior_gradient(T_var, qp);

      const Number k = _k[dim];

      const Point & p = xyz[qp];

      // solution + laplacian depend on problem dimension
      const Number u_exact = (mesh_dim == 2) ?
        std::sin(libMesh::pi*p(0)) * std::sin(libMesh::pi*p(1)) :
        std::sin(libMesh::pi*p(0)) * std::sin(libMesh::pi*p(1)) *
        std::sin(libMesh::pi*p(2));

      // Only apply forcing to interior elements
      const Number forcing = (dim == mesh_dim) ?
        -k * u_exact * (dim * libMesh::pi * libMesh::pi) : 0;

      const Number JxWxNK = JxW[qp] * -k;

      for (unsigned int i=0; i != n_T_dofs; i++)
        F(i) += JxWxNK * (grad_T * dphi[i][qp] + forcing * phi[i][qp]);
      if (request_jacobian)
        {
          const Number JxWxNKxD = JxWxNK *
            context.get_elem_solution_derivative();

          for (unsigned int i=0; i != n_T_dofs; i++)
            for (unsigned int j=0; j != n_T_dofs; ++j)
              K(i,j) += JxWxNKxD * (dphi[i][qp] * dphi[j][qp]);
        }
    } // end of the quadrature point qp-loop

  return request_jacobian;
}