fem_system.c

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

#endif
              // Logging of numerical jacobians is done separately
              this->numerical_side_jacobian();

              // If we're in DEBUG mode or if
	      // verify_analytic_jacobians is on, we've moved
              // elem_jacobian's accumulated values into old_jacobian.
              // Now let's add them back.
#ifndef DEBUG
              if (verify_analytic_jacobians != 0.0)
#endif // ifndef DEBUG
                elem_jacobian += old_jacobian;
            }

          // Compute a numeric jacobian if we're asked to verify the
          // analytic jacobian we got
	  if (get_jacobian && jacobian_computed &&
              verify_analytic_jacobians != 0.0)
            {
              DenseMatrix<Number> analytic_jacobian(elem_jacobian);

              elem_jacobian.zero();
              // Logging of numerical jacobians is done separately
              this->numerical_side_jacobian();

              Real analytic_norm = analytic_jacobian.l1_norm();
              Real numerical_norm = elem_jacobian.l1_norm();

              // If we can continue, we'll probably prefer the analytic jacobian
              analytic_jacobian.swap(elem_jacobian);

              // The matrix "analytic_jacobian" will now hold the error matrix
              analytic_jacobian.add(-1.0, elem_jacobian);
              Real error_norm = analytic_jacobian.l1_norm();

              Real relative_error = error_norm /
                                    std::max(analytic_norm, numerical_norm);

              if (relative_error > verify_analytic_jacobians)
                {
                  std::cerr << "Relative error " << relative_error
                            << " detected in analytic jacobian on element "
                            << elem->id()
			    << ", side " << side << '!' << std::endl;

                  unsigned int old_precision = std::cout.precision();
                  std::cout.precision(16);
	          std::cout << "J_analytic " << elem->id() << " = "
                            << elem_jacobian << std::endl;
                  analytic_jacobian.add(1.0, elem_jacobian);
	          std::cout << "J_numeric " << elem->id() << " = "
                            << analytic_jacobian << std::endl;
                  std::cout.precision(old_precision);

                  libmesh_error();
                }
              // Once we've verified a side, we'll want to add back the
              // rest of the accumulated jacobian
              elem_jacobian += old_jacobian;
            }
	  // In DEBUG mode, we've set elem_jacobian == 0, and we
          // may still need to add the old jacobian back
#ifdef DEBUG
	  if (get_jacobian && jacobian_computed &&
              verify_analytic_jacobians == 0.0)
            {
              elem_jacobian += old_jacobian;
            }
#endif // ifdef DEBUG
        }

#ifdef ENABLE_AMR
      // We turn off the asymmetric constraint application;
      // enforce_constraints_exactly() should be called in the solver
      if (get_residual && get_jacobian)
        this->get_dof_map().constrain_element_matrix_and_vector
          (elem_jacobian, elem_residual, dof_indices, false);
      else if (get_residual)
        this->get_dof_map().constrain_element_vector
          (elem_residual, dof_indices, false);
      else if (get_jacobian)
        this->get_dof_map().constrain_element_matrix
          (elem_jacobian, dof_indices, false);
#endif // #ifdef ENABLE_AMR

      if (get_jacobian && print_element_jacobians)
        {
          unsigned int old_precision = std::cout.precision();
          std::cout.precision(16);
	  std::cout << "J_elem " << elem->id() << " = "
                    << elem_jacobian << std::endl;
          std::cout.precision(old_precision);
        }

      if (get_jacobian)
        this->matrix->add_matrix (elem_jacobian, dof_indices);
      if (get_residual)
        this->rhs->add_vector (elem_residual, dof_indices);
    }


  if (get_residual && print_residual_norms)
    {
      this->rhs->close();
      std::cout << "|F| = " << this->rhs->l1_norm() << std::endl;
    }
  if (get_residual && print_residuals)
    {
      this->rhs->close();
      unsigned int old_precision = std::cout.precision();
      std::cout.precision(16);
      std::cout << "F = [" << *(this->rhs) << "];" << std::endl;
      std::cout.precision(old_precision);
    }
  if (get_jacobian && print_jacobian_norms)
    {
      this->matrix->close();
      std::cout << "|J| = " << this->matrix->l1_norm() << std::endl;
    }
  if (get_jacobian && print_jacobians)
    {
      this->matrix->close();
      unsigned int old_precision = std::cout.precision();
      std::cout.precision(16);
      std::cout << "J = [" << *(this->matrix) << "];" << std::endl;
      std::cout.precision(old_precision);
    }
  STOP_LOG(log_name, "FEMSystem");
}



void FEMSystem::postprocess ()
{
  START_LOG("postprocess()", "FEMSystem");

  const MeshBase& mesh = this->get_mesh();

  const DofMap& dof_map = this->get_dof_map();

//  this->get_vector("_nonlinear_solution").localize
//    (*current_local_nonlinear_solution,
//     dof_map.get_send_list());
  this->update();

  // Loop over every active mesh element on this processor
  MeshBase::const_element_iterator el =
    mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el =
    mesh.active_local_elements_end();

  for ( ; el != end_el; ++el)
    {
      elem = *el;

      // Initialize the per-element data for elem.
      dof_map.dof_indices (elem, dof_indices);
      unsigned int n_dofs = dof_indices.size();

      elem_solution.resize(n_dofs);
      // This resize call also zeros out the residual
      elem_residual.resize(n_dofs);

      for (unsigned int i=0; i != n_dofs; ++i)
//        elem_solution(i) = current_nonlinear_solution(dof_indices[i]);
        elem_solution(i) = current_solution(dof_indices[i]);

      // Initialize the per-variable data for elem.
      unsigned int sub_dofs = 0;
      for (unsigned int i=0; i != this->n_vars(); ++i)
        {
          dof_map.dof_indices (elem, dof_indices_var[i], i);

          elem_subsolutions[i]->reposition
            (sub_dofs, dof_indices_var[i].size());

          if (use_fixed_solution)
            elem_fixed_subsolutions[i]->reposition
              (sub_dofs, dof_indices_var[i].size());

          elem_subresiduals[i]->reposition
            (sub_dofs, dof_indices_var[i].size());

          sub_dofs += dof_indices_var[i].size();
        }
      libmesh_assert(sub_dofs == n_dofs);

      // Optionally initialize all the interior FE objects on elem.
      // Logging of FE::reinit is done in the FE functions
      if (fe_reinit_during_postprocess)
        {
          std::map<FEType, FEBase *>::iterator fe_end = element_fe.end();
          for (std::map<FEType, FEBase *>::iterator i = element_fe.begin();
               i != fe_end; ++i)
            {
              i->second->reinit(elem);
            }
        }
      
      this->element_postprocess();

      for (side = 0; side != elem->n_sides(); ++side)
        {
          // Don't compute on non-boundary sides unless requested
          if (!postprocess_sides ||
              (!compute_internal_sides &&
               elem->neighbor(side) != NULL))
            continue;

          // Optionally initialize all the interior FE objects on elem/side.
          // Logging of FE::reinit is done in the FE functions
          if (fe_reinit_during_postprocess)
            {
              std::map<FEType, FEBase *>::iterator fe_end = element_fe.end();
              fe_end = side_fe.end();
              for (std::map<FEType, FEBase *>::iterator i = side_fe.begin();
                   i != fe_end; ++i)
                {
                  i->second->reinit(elem, side);
                }
            }

          this->side_postprocess();
        }
    }

  STOP_LOG("postprocess()", "FEMSystem");
}



void FEMSystem::numerical_elem_jacobian ()
{
  START_LOG("numerical_elem_jacobian()", "FEMSystem");

  DenseVector<Number> original_residual(elem_residual);
  DenseVector<Number> backwards_residual(elem_residual);
  DenseMatrix<Number> numerical_jacobian(elem_jacobian);

  for (unsigned int j = 0; j != dof_indices.size(); ++j)
    {
      Number original_solution = elem_solution(j);
      elem_solution(j) -= numerical_jacobian_h;
      elem_residual.zero();
      time_solver->element_residual(false);

      backwards_residual = elem_residual;
      elem_solution(j) = original_solution + numerical_jacobian_h;
      elem_residual.zero();
      time_solver->element_residual(false);

      for (unsigned int i = 0; i != dof_indices.size(); ++i)
        {
          numerical_jacobian(i,j) =
            (elem_residual(i) - backwards_residual(i)) /
            2. / numerical_jacobian_h;
        }
      elem_solution(j) = original_solution;
    }

  elem_residual = original_residual;
  elem_jacobian = numerical_jacobian;
  STOP_LOG("numerical_elem_jacobian()", "FEMSystem");
}



void FEMSystem::numerical_side_jacobian ()
{
  START_LOG("numerical_side_jacobian()", "FEMSystem");

  DenseVector<Number> original_residual(elem_residual);
  DenseVector<Number> backwards_residual(elem_residual);
  DenseMatrix<Number> numerical_jacobian(elem_jacobian);

  for (unsigned int j = 0; j != dof_indices.size(); ++j)
    {
      Number original_solution = elem_solution(j);
      elem_solution(j) -= numerical_jacobian_h;
      elem_residual.zero();
      time_solver->side_residual(false);

      backwards_residual = elem_residual;
      elem_solution(j) = original_solution + numerical_jacobian_h;
      elem_residual.zero();
      time_solver->side_residual(false);

      for (unsigned int i = 0; i != dof_indices.size(); ++i)
        {
          numerical_jacobian(i,j) +=
            (elem_residual(i) - backwards_residual(i))
            / 2. / numerical_jacobian_h;
        }
      elem_solution(j) = original_solution;
    }

  elem_residual = original_residual;
  elem_jacobian = numerical_jacobian;
  STOP_LOG("numerical_side_jacobian()", "FEMSystem");
}



void FEMSystem::time_evolving (unsigned int var)
{
  // Call the parent function
  Parent::time_evolving(var);

  // Then make sure we're prepared to do mass integration
  element_fe_var[var]->get_JxW();
  element_fe_var[var]->get_phi();
}



bool FEMSystem::mass_residual (bool request_jacobian)
{
  unsigned int n_qpoints = element_qrule->n_points();

  for (unsigned int var = 0; var != this->n_vars(); ++var)
    {
      if (!_time_evolving[var])
        continue;

      const std::vector<Real> &JxW = 
        element_fe_var[var]->get_JxW();

      const std::vector<std::vector<Real> > &phi =
        element_fe_var[var]->get_phi();

      const unsigned int n_dofs = dof_indices_var[var].size();

      DenseSubVector<Number> &Fu = *elem_subresiduals[var];
      DenseSubMatrix<Number> &Kuu = *elem_subjacobians[var][var];

      for (unsigned int qp = 0; qp != n_qpoints; ++qp)
        {
          Number u = interior_value(var, qp);
          Number JxWxU = JxW[qp] * u;
          for (unsigned int i = 0; i != n_dofs; ++i)
            {
              Fu(i) += JxWxU * phi[i][qp];
              if (request_jacobian && elem_solution_derivative)
                {
                  libmesh_assert (elem_solution_derivative == 1.0);

                  Number JxWxPhiI = JxW[qp] * phi[i][qp];
                  Kuu(i,i) += JxWxPhiI * phi[i][qp];
                  for (unsigned int j = i+1; j != n_dofs; ++j)
                    {
                      Number Kij = JxWxPhiI * phi[j][qp];
                      Kuu(i,j) += Kij;
                      Kuu(j,i) += Kij;
                    }
                }
            }
        }
    }

  return request_jacobian;
}