naviersystem.c

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

  DenseSubMatrix<Number> &Kpu = *elem_subjacobians[p_var][u_var];
  DenseSubMatrix<Number> &Kpv = *elem_subjacobians[p_var][v_var];
  DenseSubMatrix<Number> &Kpw = *elem_subjacobians[p_var][w_var];
  DenseSubVector<Number> &Fp = *elem_subresiduals[p_var];

  // Add the constraint given by the continuity equation
  unsigned int n_qpoints = element_qrule->n_points();
  for (unsigned int qp=0; qp != n_qpoints; qp++)
    {
      // Compute the velocity gradient at the old Newton iterate
      Gradient grad_u = interior_gradient(u_var, qp),
               grad_v = interior_gradient(v_var, qp),
               grad_w = interior_gradient(w_var, qp);

      // Now a loop over the pressure degrees of freedom.  This
      // computes the contributions of the continuity equation.
      for (unsigned int i=0; i != n_p_dofs; i++)
        {
          Fp(i) += JxW[qp] * psi[i][qp] *
                   (grad_u(0) + grad_v(1));
          if (dim == 3)
            Fp(i) += JxW[qp] * psi[i][qp] *
                     (grad_w(2));

          if (request_jacobian)
            for (unsigned int j=0; j != n_u_dofs; j++)
              {
                Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
                Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
                if (dim == 3)
                  Kpw(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](2);
              }
        }
    } // end of the quadrature point qp-loop

  return request_jacobian;
}

      

bool NavierSystem::side_constraint (bool request_jacobian)
{
  // Here we define some references to cell-specific data that
  // will be used to assemble the linear system.

  // Element Jacobian * quadrature weight for side integration
  const std::vector<Real> &JxW_side = fe_side_vel->get_JxW();

  // The velocity shape functions at side quadrature points.
  const std::vector<std::vector<Real> >& phi_side =
    fe_side_vel->get_phi();

  // The number of local degrees of freedom in u and v
  const unsigned int n_u_dofs = dof_indices_var[u_var].size(); 

  // The subvectors and submatrices we need to fill:
  const unsigned int dim = this->get_mesh().mesh_dimension();
  DenseSubMatrix<Number> &Kuu = *elem_subjacobians[u_var][u_var];
  DenseSubMatrix<Number> &Kvv = *elem_subjacobians[v_var][v_var];
  DenseSubMatrix<Number> &Kww = *elem_subjacobians[w_var][w_var];

  DenseSubVector<Number> &Fu = *elem_subresiduals[u_var];
  DenseSubVector<Number> &Fv = *elem_subresiduals[v_var];
  DenseSubVector<Number> &Fw = *elem_subresiduals[w_var];

  // For this example we will use Dirichlet velocity boundary
  // conditions imposed at each timestep via the penalty method.

  // The penalty value.  \f$ \frac{1}{\epsilon} \f$
  const Real penalty = 1.e10;

  unsigned int n_sidepoints = side_qrule->n_points();
  for (unsigned int qp=0; qp != n_sidepoints; qp++)
    {
      // Compute the solution at the old Newton iterate
      Number u = side_value(u_var, qp),
             v = side_value(v_var, qp),
             w = side_value(w_var, qp);

      // Set u = 1 on the top boundary, (which build_square() has
      // given boundary id 2, and build_cube() has called 5),
      // u = 0 everywhere else
      short int boundary_id =
        this->get_mesh().boundary_info->boundary_id(elem, side);
      libmesh_assert (boundary_id != BoundaryInfo::invalid_id);
      const short int top_id = (dim==3) ? 5 : 2;

      Real u_value = 0.;

      // For lid-driven cavity, set u=1 on the lid.
      if ((application == 0) && (boundary_id == top_id))
	u_value = 1.;
              
      // Set v, w = 0 everywhere
      const Real v_value = 0.;
      const Real w_value = 0.;

      for (unsigned int i=0; i != n_u_dofs; i++)
        {
          Fu(i) += JxW_side[qp] * penalty *
                   (u - u_value) * phi_side[i][qp];
          Fv(i) += JxW_side[qp] * penalty *
                   (v - v_value) * phi_side[i][qp];
          if (dim == 3)
            Fw(i) += JxW_side[qp] * penalty *
                     (w - w_value) * phi_side[i][qp];

          if (request_jacobian)
            for (unsigned int j=0; j != n_u_dofs; j++)
              {
                Kuu(i,j) += JxW_side[qp] * penalty *
                            phi_side[i][qp] * phi_side[j][qp];
                Kvv(i,j) += JxW_side[qp] * penalty *
                            phi_side[i][qp] * phi_side[j][qp];
                if (dim == 3)
                  Kww(i,j) += JxW_side[qp] * penalty *
                              phi_side[i][qp] * phi_side[j][qp];
              }
        }
    }

  // Pin p = 0 at the origin
  if (elem->contains_point(Point(0.,0.)))
    {
      // The pressure penalty value.  \f$ \frac{1}{\epsilon} \f$
      const Real penalty = 1.e9;

      DenseSubMatrix<Number> &Kpp = *elem_subjacobians[p_var][p_var];
      DenseSubVector<Number> &Fp = *elem_subresiduals[p_var];
      const unsigned int n_p_dofs = dof_indices_var[p_var].size(); 

      Point zero(0.);
      Number p = point_value(p_var, zero);
      Number p_value = 0.;

      unsigned int dim = get_mesh().mesh_dimension();
      FEType fe_type = element_fe_var[p_var]->get_fe_type();
      Point p_master = FEInterface::inverse_map(dim, fe_type, elem, zero);

      std::vector<Real> point_phi(n_p_dofs);
      for (unsigned int i=0; i != n_p_dofs; i++)
        {
          point_phi[i] = FEInterface::shape(dim, fe_type, elem, i, p_master);
        }

      for (unsigned int i=0; i != n_p_dofs; i++)
        {
          Fp(i) += penalty * (p - p_value) * point_phi[i];
          if (request_jacobian)
            for (unsigned int j=0; j != n_p_dofs; j++)
              {
		Kpp(i,j) += penalty * point_phi[i] * point_phi[j];
              }
        }
    }

  return request_jacobian;
}


// We override the default mass_residual function,
// because in the non-dimensionalized Navier-Stokes equations
// the time derivative of velocity has a Reynolds number coefficient.
// Alternatively we could divide the whole equation by
// Reynolds number (or choose a more complicated non-dimensionalization
// of time), but this gives us an opportunity to demonstrate overriding
// FEMSystem::mass_residual()
bool NavierSystem::mass_residual (bool request_jacobian)
{
  // The subvectors and submatrices we need to fill:
  const unsigned int dim = this->get_mesh().mesh_dimension();

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

  // The velocity shape functions at interior quadrature points.
  const std::vector<std::vector<Real> >& phi = fe_velocity->get_phi();

  // The subvectors and submatrices we need to fill:
  DenseSubVector<Number> &Fu = *elem_subresiduals[u_var];
  DenseSubVector<Number> &Fv = *elem_subresiduals[v_var];
  DenseSubVector<Number> &Fw = *elem_subresiduals[w_var];
  DenseSubMatrix<Number> &Kuu = *elem_subjacobians[u_var][u_var];
  DenseSubMatrix<Number> &Kvv = *elem_subjacobians[v_var][v_var];
  DenseSubMatrix<Number> &Kww = *elem_subjacobians[w_var][w_var];

  // The number of local degrees of freedom in velocity
  const unsigned int n_u_dofs = dof_indices_var[u_var].size();

  unsigned int n_qpoints = element_qrule->n_points();

  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
    {
      Number u = interior_value(u_var, qp),
             v = interior_value(v_var, qp),
             w = interior_value(w_var, qp);

      // We pull as many calculations as possible outside of loops
      Number JxWxRe   = JxW[qp] * Reynolds;
      Number JxWxRexU = JxWxRe * u;
      Number JxWxRexV = JxWxRe * v;
      Number JxWxRexW = JxWxRe * w;

      for (unsigned int i = 0; i != n_u_dofs; ++i)
        {
          Fu(i) += JxWxRexU * phi[i][qp];
          Fv(i) += JxWxRexV * phi[i][qp];
          if (dim == 3)
            Fw(i) += JxWxRexW * phi[i][qp];

          if (request_jacobian && elem_solution_derivative)
            {
              libmesh_assert (elem_solution_derivative == 1.0);

              Number JxWxRexPhiI = JxWxRe * phi[i][qp];
              Number JxWxRexPhiII = JxWxRexPhiI * phi[i][qp];
              Kuu(i,i) += JxWxRexPhiII;
              Kvv(i,i) += JxWxRexPhiII;
              if (dim == 3)
                Kww(i,i) += JxWxRexPhiII;

              // The mass matrix is symmetric, so we calculate
              // one triangle and add it to both upper and lower
              // triangles
              for (unsigned int j = i+1; j != n_u_dofs; ++j)
                {
                  Number Kij = JxWxRexPhiI * phi[j][qp];
                  Kuu(i,j) += Kij;
                  Kuu(j,i) += Kij;
                  Kvv(i,j) += Kij;
                  Kvv(j,i) += Kij;
                  if (dim == 3)
                    {
                      Kww(i,j) += Kij;
                      Kww(j,i) += Kij;
                    }
                }
            }
        }
    }

  return request_jacobian;
}



Point NavierSystem::forcing(const Point& p)
{
  switch (application)
    {
      // lid driven cavity
    case 0:
      {
	// No forcing
	return Point(0.,0.,0.);
      }

      // Homogeneous Dirichlet BCs + sinusoidal forcing
    case 1:
      {
	const unsigned int dim = this->get_mesh().mesh_dimension();

	// This assumes your domain is defined on [0,1]^dim.
	Point f;

	// Counter-Clockwise vortex in x-y plane
	if (dim==2)
	  {
	    f(0) =  std::sin(2.*libMesh::pi*p(1));
	    f(1) = -std::sin(2.*libMesh::pi*p(0));
	  }

	// Counter-Clockwise vortex in x-z plane
	else if (dim==3)
	  {
	    f(0) =  std::sin(2.*libMesh::pi*p(1));
	    f(1) = 0.;
	    f(2) = -std::sin(2.*libMesh::pi*p(0));
	  }

	return f;
      }

    default:
      {
	libmesh_error();
      }
    }
}