C++ GenericMatrix::setrow方法代码示例

// Base case for all divergence computations. 
// Compute divergence of vector field u.
void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A,
                       const FunctionSpace& DGscalar, 
                       const FunctionSpace& CRvector)
{
  std::shared_ptr<const GenericDofMap>
    CR1_dofmap = CRvector.dofmap(),
    DG0_dofmap = DGscalar.dofmap();

  // Figure out about the local dofs of DG0 
  std::pair<std::size_t, std::size_t>
  first_last_dof = DG0_dofmap->ownership_range();
  std::size_t first_dof = first_last_dof.first;
  std::size_t last_dof = first_last_dof.second;
  std::size_t n_local_dofs = last_dof - first_dof;

  // Get topological dimension so that we know what Facet is
  const Mesh mesh = *DGscalar.mesh();
  std::size_t tdim = mesh.topology().dim(); 
  std::size_t gdim = mesh.geometry().dim();
  
  std::vector<std::size_t> columns;
  std::vector<double> values;
  
  // Fill the values
  for(CellIterator cell(mesh); !cell.end(); ++cell)
  {
    auto dg_dofs = DG0_dofmap->cell_dofs(cell->index());
    // There is only one DG0 dof per cell
    dolfin::la_index cell_dof = dg_dofs[0];
    
    Point cell_mp = cell->midpoint();
    double cell_volume = cell->volume();
    std::size_t local_facet_index = 0;      
    
    auto cr_dofs = CR1_dofmap->cell_dofs(cell->index());
    
    for(FacetIterator facet(*cell); !facet.end(); ++facet)
    {
      double facet_measure=0;
      if(tdim == 2)
        facet_measure = Edge(mesh, facet->index()).length();
      else if(tdim == 3)
        facet_measure = Face(mesh, facet->index()).area();
      // Tdim 1 will not happen because CR is not defined there
      
      Point facet_normal = facet->normal();

      // Flip the normal if it is not outer already
      Point facet_mp = facet->midpoint();
      double sign = (facet_normal.dot(facet_mp - cell_mp) > 0.0) ? 1.0 : -1.0;
      facet_normal *= (sign*facet_measure/cell_volume);
      
      // Dofs of CR on the facet, local order
      std::vector<std::size_t> facet_dofs;
      CR1_dofmap->tabulate_facet_dofs(facet_dofs, local_facet_index);
      
      for (std::size_t j = 0 ; j < facet_dofs.size(); j++)
      {   
        columns.push_back(cr_dofs[facet_dofs[j]]);
        values.push_back(facet_normal[j]);
      }        
      local_facet_index += 1;
    }
    M.setrow(cell_dof, columns, values);
    columns.clear();
    values.clear();
  }
  M.apply("insert");
  //std::shared_ptr<GenericMatrix> Cp = MatMatMult(M, A);
  //return Cp;
}

//.........这里部分代码省略.........
      
      for (std::size_t i = 0; i < columns.size(); i++)
      {
        std::size_t dof = columns[i];
        if (dof < weight_range.first || dof >= weight_range.second)
        {
          std::size_t owner = dof_owner(all_ranges, dof);
          dofs_needed[owner].push_back(dof);
        }
      }
    }

    // Communicate to all which weights are needed by the process
    std::vector<std::vector<std::size_t> > dofs_needed_recv;
    MPI::all_to_all(mpi_comm, dofs_needed, dofs_needed_recv);
    
    // Fetch the weights that must be communicated
    std::vector<std::vector<double> > weights_to_send(num_processes);    
    for (std::size_t p = 0; p < num_processes; p++)
    {
      if (p == MPI::rank(mpi_comm))  
        continue;
      
      std::vector<std::size_t> dofs = dofs_needed_recv[p];
      std::map<std::size_t, double> send_weights;
      for (std::size_t k = 0; k < dofs.size(); k++)
      {
        weights_to_send[p].push_back(weight[dofs[k]-weight_range.first]);
      }
    }
    std::vector<std::vector<double> > weights_to_send_recv;
    MPI::all_to_all(mpi_comm, weights_to_send, weights_to_send_recv);
    
    // Create a map for looking up received weights
    std::map<std::size_t, double> received_weights;
    for (std::size_t p = 0; p < num_processes; p++)
    {
      if (p == MPI::rank(mpi_comm))
        continue;
      
      for (std::size_t k = 0; k < dofs_needed[p].size(); k++)
      {
        received_weights[dofs_needed[p][k]] = weights_to_send_recv[p][k];         
      }
    }
    
    for (std::size_t row = 0; row < m; row++)
    {   
      // Get global row number
      const std::size_t global_row = row + row_range.first;
      
      A.getrow(global_row, columns, values);
      for (std::size_t i = 0; i < values.size(); i++)
      {
        std::size_t dof = columns[i];
        if (dof < weight_range.first || dof >= weight_range.second)
        {
          values[i] = received_weights[dof];
        }
        else
        {
          values[i] = weight[columns[i]-weight_range.first];  
        }
//        values[i] = 1./values[i];
      }
      
      double s = std::accumulate(values.begin(), values.end(), 0.0);
      std::transform(values.begin(), values.end(), values.begin(),
                     std::bind2nd(std::multiplies<double>(), 1./s));      

      for (std::size_t i=0; i<values.size(); i++)
      {
        double w;
        std::size_t dof = columns[i];
        if (dof < weight_range.first || dof >= weight_range.second)
        {
          w = received_weights[dof];
        }
        else
        {
          w = weight[dof-weight_range.first];  
        }        
        values[i] = values[i]*w;
//        values[i] = values[i]*values[i];
        
      }     
      
      allvalues.push_back(values);
      allcolumns.push_back(columns);
    }

    for (std::size_t row = 0; row < m; row++)
    {       
      // Get global row number
      const std::size_t global_row = row + row_range.first;
      
      A.setrow(global_row, allcolumns[row], allvalues[row]);
    }
    A.apply("insert");  
  }  

//-----------------------------------------------------------------------------
void DirichletBC::zero_columns(GenericMatrix& A,
                               GenericVector& b,
                               double diag_val) const
{
  // Check arguments
  check_arguments(&A, &b, NULL, 1);

  // A map to hold the mapping from boundary dofs to boundary values
  Map bv_map;
  get_boundary_values(bv_map);

  // Create lookup table of dofs
  //const std::size_t nrows = A.size(0); // should be equal to b.size()
  const std::size_t ncols = A.size(1); // should be equal to max possible dof+1

  std::pair<std::size_t, std::size_t> rows = A.local_range(0);

  std::vector<char> is_bc_dof(ncols);
  std::vector<double> bc_dof_val(ncols);
  for (Map::const_iterator bv = bv_map.begin(); bv != bv_map.end(); ++bv)
  {
    is_bc_dof[bv->first] = 1;
    bc_dof_val[bv->first] = bv->second;
  }

  // Scan through all columns of all rows, setting to zero if
  // is_bc_dof[column]. At the same time, we collect corrections to
  // the RHS

  std::vector<std::size_t> cols;
  std::vector<double> vals;
  std::vector<double> b_vals;
  std::vector<dolfin::la_index> b_rows;

  for (std::size_t row = rows.first; row < rows.second; row++)
  {
    // If diag_val is nonzero, the matrix is a diagonal block
    // (nrows==ncols), and we can set the whole BC row
    if (diag_val != 0.0 && is_bc_dof[row])
    {
      A.getrow(row, cols, vals);
      for (std::size_t j = 0; j < cols.size(); j++)
        vals[j] = (cols[j] == row)*diag_val;
      A.setrow(row, cols, vals);
      A.apply("insert");
      b.setitem(row, bc_dof_val[row]*diag_val);
    }
    else // Otherwise, we scan the row for BC columns
    {
      A.getrow(row, cols, vals);
      bool row_changed = false;
      for (std::size_t j = 0; j < cols.size(); j++)
      {
        const std::size_t col = cols[j];

        // Skip columns that aren't BC, and entries that are zero
        if (!is_bc_dof[col] || vals[j] == 0.0)
          continue;

        // We're going to change the row, so make room for it
        if (!row_changed)
        {
          row_changed = true;
          b_rows.push_back(row);
          b_vals.push_back(0.0);
        }

        b_vals.back() -= bc_dof_val[col]*vals[j];
        vals[j] = 0.0;
      }
      if (row_changed)
      {
        A.setrow(row, cols, vals);
        A.apply("insert");
      }
    }
  }

  b.add_local(&b_vals.front(), b_rows.size(), &b_rows.front());
  b.apply("add");
}