C++ SparseMatrix::add方法代码示例

void AssembleOptimization::hessian (const NumericVector<Number> & /*soln*/,
                                    SparseMatrix<Number> & H_f,
                                    OptimizationSystem & sys)
{
  H_f.zero();
  H_f.add(1., *A_matrix);

  // We also need to add the Hessian of the inequality and equality constraints,
  //  \sum_i^n_eq lambda_eq_i H_eq_i + \sum_i^n_ineq lambda_ineq_i H_ineq_i
  // where lambda_eq and lambda_ineq denote Lagrange multipliers associated
  // with the equality and inequality constraints, respectively.
  // In this example, our equality constraints are linear, hence H_eq_i = 0.
  // However, our inequality constraint is nonlinear, so it will contribute
  // to the Hessian matrix.
  sys.optimization_solver->get_dual_variables();
  dof_id_type ineq_index = 0;
  Number lambda_ineq_0 = 0.;
  unsigned int lambda_rank = 0;
  if ((sys.lambda_ineq->first_local_index() <= ineq_index) &&
      (ineq_index < sys.lambda_ineq->last_local_index()))
    {
      lambda_ineq_0 = (*sys.lambda_ineq)(0);
      lambda_rank = sys.comm().rank();
    }

  // Sync lambda_rank across all processors.
  sys.comm().sum(lambda_rank);
  sys.comm().broadcast(lambda_rank, lambda_rank);

  if ((sys.get_dof_map().first_dof() <= 200) && (200 < sys.get_dof_map().end_dof()))
    H_f.add(200, 200, 2. * lambda_ineq_0);
}

void CSFSolver::compute_system_matrices(const Vector<double> &last_relax_step,
                                        const Vector<double> &previous_time_step,
                                        SparseMatrix<double> &mass_output,
                                        SparseMatrix<double> &stiffness_output) {
    mass_output = 0;
    stiffness_output = 0;
    const QGauss<2> quadrature_formula(3);
    FEValues<2> fe_values(fe, quadrature_formula,
        update_values | update_gradients | update_JxW_values | update_quadrature_points);
    const unsigned int dofs_per_cell = fe.dofs_per_cell;
    const unsigned int n_q_points = quadrature_formula.size();

    FullMatrix<double> local_mass_matrix(dofs_per_cell, dofs_per_cell);
    FullMatrix<double> local_stiffness_matrix(dofs_per_cell, dofs_per_cell);
    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

    std::vector<Tensor<1, 2, double>> last_relax_grads(n_q_points);
    std::vector<Tensor<1, 2, double>> previous_step_grads(n_q_points);

    for(auto cell : dof_handler.active_cell_iterators()) {
        local_mass_matrix = 0;
        local_stiffness_matrix = 0;
        fe_values.reinit(cell);
        fe_values.get_function_gradients(last_relax_step, last_relax_grads);
        fe_values.get_function_gradients(previous_time_step, previous_step_grads);

        for(unsigned int q_pt_idx = 0; q_pt_idx < n_q_points; q_pt_idx++) {
            for(unsigned int i = 0; i < dofs_per_cell; i++) {
                for(unsigned int j = 0; j < dofs_per_cell; j++) {
                    // Compute the shared denominator - normalizing to avoid infinities
                    double grad_sum = last_relax_grads[q_pt_idx].norm() + previous_step_grads[q_pt_idx].norm();
                    double grad_sum_inv = 1.0 / MAX(grad_sum, grad_epsilon);
                    // If there is a weight function, compute its value at the quad point
                    double weight_fn_value = weight_function(fe_values.quadrature_point(q_pt_idx));
                    // Compute the mass component
                    local_mass_matrix(i, j) += (fe_values.shape_value(i, q_pt_idx) * fe_values.shape_value(j, q_pt_idx) * weight_fn_value * grad_sum_inv * fe_values.JxW(q_pt_idx));
                    // Compute the stiffness component
                    local_stiffness_matrix(i, j) += (fe_values.shape_grad(i, q_pt_idx) * fe_values.shape_grad(j, q_pt_idx) * weight_fn_value * grad_sum_inv * fe_values.JxW(q_pt_idx));
                }
            }
        }

        // Copy the local matrices into the global output
        cell->get_dof_indices(local_dof_indices);
        for(unsigned int i = 0; i < dofs_per_cell; i++) {
            for(unsigned int j = 0; j < dofs_per_cell; j++) {
                mass_output.add(local_dof_indices[i], local_dof_indices[j], local_mass_matrix(i, j));
                stiffness_output.add(local_dof_indices[i], local_dof_indices[j], local_stiffness_matrix(i, j));
            }
        }
    }
}

void
NodalConstraint::computeJacobian(SparseMatrix<Number> & jacobian)
{
  _qp = 0;
  dof_id_type & dof_idx = _var.nodalDofIndex();
  dof_id_type & dof_idx_neighbor = _var.nodalDofIndexNeighbor();

  Real scaling_factor = _var.scalingFactor();
  jacobian.add(dof_idx, dof_idx, scaling_factor * computeQpJacobian(Moose::MasterMaster));
  jacobian.add(dof_idx, dof_idx_neighbor, scaling_factor * computeQpJacobian(Moose::MasterSlave));
  jacobian.add(dof_idx_neighbor, dof_idx_neighbor, scaling_factor * computeQpJacobian(Moose::SlaveSlave));
  jacobian.add(dof_idx_neighbor, dof_idx, scaling_factor * computeQpJacobian(Moose::SlaveMaster));
}

void AssembleOptimization::hessian (const NumericVector<Number> & /*soln*/,
                                    SparseMatrix<Number> & H_f,
                                    OptimizationSystem & /*sys*/)
{
  H_f.zero();
  H_f.add(1., *A_matrix);
}

void
Assembly::addCachedJacobian(SparseMatrix<Number> & jacobian)
{
  mooseAssert(_cached_jacobian_rows.size() == _cached_jacobian_cols.size(),
              "Error: Cached data sizes MUST be the same!");

  for(unsigned int i=0; i<_cached_jacobian_rows.size(); i++)
    jacobian.add(_cached_jacobian_rows[i], _cached_jacobian_cols[i], _cached_jacobian_values[i]);

  if (_max_cached_jacobians < _cached_jacobian_values.size())
    _max_cached_jacobians = _cached_jacobian_values.size();

  // Try to be more efficient from now on
  // The 2 is just a fudge factor to keep us from having to grow the vector during assembly
  _cached_jacobian_values.clear();
  _cached_jacobian_values.reserve(_max_cached_jacobians*2);

  _cached_jacobian_rows.clear();
  _cached_jacobian_rows.reserve(_max_cached_jacobians*2);

  _cached_jacobian_cols.clear();
  _cached_jacobian_cols.reserve(_max_cached_jacobians*2);
}

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh;
  H2DReader mloader;
  // We load the mesh on a (-1, 1)^2 domain.
  mloader.load("ref_square.mesh", &mesh);            

  // Enter boundary markers 
  // (If no markers are entered, default is a natural BC).
  BCTypes bc_types;

  // Create an H1 space with default shapeset,
  // natural BC, and linear elements.
  H1Space space(&mesh, &bc_types, P_INIT);
  // The type of element, mesh_mode = 4 means a rectangle element.
  int mesh_mode = 4;
  int ndof = Space::get_num_dofs(&space);
  printf("ndof = %d\n", ndof);
  if(ndof > FNS_NUM) error("Max number of shape functions exceeded.");

  int *fn_idx = new int [FNS_NUM];
  int m = 0;
  int order = P_INIT;

  // Get the vertex fns index.
  info("Get the vertex fns index.");
  for (int i = 0; i < mesh_mode; i++, m++)
  {
    fn_idx[m] = space.get_shapeset()->get_vertex_index(i);
    printf("m = %d, get_vertex_index(m) = %d\n",
          m, space.get_shapeset()->get_vertex_index(m));
  }
  // Get the edge fns index.
  info("Get the edge fns index.");
  for (int edge_order = 2; edge_order <= order; edge_order++)
  {
    for (int j = 0; j < mesh_mode; j++, m++)
    {
      fn_idx[m] = space.get_shapeset()->get_edge_index(j, 0, edge_order);
      printf("m = %d, get_edge_index(m) = %d\n", m, fn_idx[m]);
    }
  }
  // Get the bubble fns index.
  info("Get the bubble fns index.");
  int number_bubble = space.get_shapeset()->get_num_bubbles(H2D_MAKE_QUAD_ORDER(order, order));
  int *bubble_idx = space.get_shapeset()->get_bubble_indices(H2D_MAKE_QUAD_ORDER(order, order));
  for (int i = 0; i < number_bubble; i++, m++ )
  {
    fn_idx[m] = bubble_idx[i];
    printf("m = %d, get_bubble_index(m) = %d\n", m, fn_idx[m]);
  }
  
  // Initialize the matrix solver.
  SparseMatrix* mat = create_matrix(matrix_solver);
  Vector* rhs = create_vector(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, mat, rhs);

  // precalc structure
  mat->prealloc(ndof);
  for (int i = 0; i < ndof; i++)
    for (int j = 0; j < ndof; j++)
      mat->pre_add_ij(i, j);
        
  mat->alloc();
  rhs->alloc(ndof);

  info("Assembling matrix ...");
  for (int i = 0; i < ndof; i++)
  {
    for (int j = 0; j < order; j++)
    {
      double x1 = -1 + (2.0/order)*j;
      double y1 = -1;
      double value = space.get_shapeset()->get_fn_value(fn_idx[i], x1, y1, 0);
      mat->add(i, j, value);
      printf("get fn[%d] value = %f  ", i, value); 
      printf("x1 = %f, y1 = %f\n", x1, y1);
    }
    for (int j = 0; j < order; j++)
    {
      double y2 = -1 + (2.0/order)*j;;
      double x2 = 1;
      double value = space.get_shapeset()->get_fn_value(fn_idx[i], x2, y2, 0);
      mat->add(i, j+order, value);
      printf("get fn[%d] value = %f  ", i, value); 
      printf("x2 = %f, y2 = %f\n", x2, y2);
    }
    for (int j = 0; j < order; j++)
    {
      double x3 = 1 + (-1.0)*(2.0/order)*j;
      double y3 = 1;
      double value = space.get_shapeset()->get_fn_value(fn_idx[i], x3, y3, 0);
      mat->add(i, j+2*order, value);
      printf("get fn[%d] value = %f  ", i, value); 
      printf("x3 = %f, y3 = %f\n", x3, y3);
    }
    for (int j = 0; j < order; j++)
    {
      double x4 = -1;
//.........这里部分代码省略.........