assemblysystembase.cpp

Go to the documentation of this file.

#include "assemblysystembase.h"
#include <numerics/matrices.h>
void AssemblySystemBase::applyBoundaryValue(MapIntDouble &boundary_values,SparseMatrix<double> &matrix,double factor)
  {
    // Require that diagonals are first
    // in each row

    if (matrix.get_sparsity_pattern().optimize_diagonal())
    {
      if (boundary_values.size() == 0)
        return;


      std::map<unsigned int,double>::const_iterator dof  = boundary_values.begin(),
        endd = boundary_values.end();
      const SparsityPattern    &sparsity    = matrix.get_sparsity_pattern();
      const unsigned int *sparsity_rowstart = (const unsigned int*) sparsity.get_rowstart_indices();
      for (; dof != endd; ++dof)
      {
        Assert (dof->first < matrix.m(), ExcInternalError());
        const unsigned int dof_number = dof->first;
        const unsigned int last = sparsity_rowstart[dof_number+1];
        for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
          matrix.global_entry(j) = 0.;

        matrix.set (dof_number, dof_number,factor);
      }
    }
   else
    {
      std::map<unsigned int,double>::const_iterator dof  = boundary_values.begin(),
        endd = boundary_values.end();
      for (; dof != endd; ++dof)
      {
        SparseMatrix<double>::iterator it     = matrix.begin(dof->first);
        SparseMatrix<double>::iterator itend  = matrix.end(dof->first);
        for (;it!=itend;it++)
        {
          it->value()=0.0;
        }
        if (factor != 0.0) matrix.set (dof->first,dof->first, factor);
      }

    }
  }
  void AssemblySystemBase::applyBoundaryValue(MapIntDouble &boundary_values,VecDouble &right_hand_side,double factor)
  {

    if (boundary_values.size() == 0)
      return;
    std::map<unsigned int,double>::const_iterator dof  = boundary_values.begin(),
      endd = boundary_values.end();
    for (; dof != endd; ++dof)
    {
      right_hand_side(dof->first) = dof->second*factor;
    }
  }
 void AssemblySystemBase::applyBoundaryValue(MapIntDouble &boundary_values,BlockVecDouble &right_hand_side,double factor) 
 {
    if (boundary_values.size() == 0)
      return;
    std::map<unsigned int,double>::const_iterator dof  = boundary_values.begin(),
      endd = boundary_values.end();
    for (; dof != endd; ++dof)
    {
      right_hand_side(dof->first) = dof->second*factor;
    }

 }

/***************************************************
Objetivo:  Aplicar condicoes de contorno mantendo a simetria do
           sistema. Para isso a linha e a coluna correspondente
           a todo grau de liberdade em que prescrevemos algum
           valor deve ser considerada nula exceto o elemento diagonal,
           de forma que as condicoes de fronteira seja satisfeita.

********************************************************/
void AssemblySystemBase::applySymetricBoundaryValue(MapIntDouble &boundary_values,SparseMatrix<double> &matrix, VecDouble &right_hand_side)
{

  // Require that diagonals are first
  // in each row
  assert (matrix.get_sparsity_pattern().optimize_diagonal());
  assert (matrix.n() == right_hand_side.size());
  // if no boundary values are to be applied
  // simply return
  if (boundary_values.size() == 0)
    return;


  const unsigned int n_dofs = matrix.m();

  // if a diagonal entry is zero
  // later, then we use another
  // number instead. take it to be
  // the first nonzero diagonal
  // element of the matrix, or 1 if
  // there is no such thing
  double first_nonzero_diagonal_entry = 1;
  for (unsigned int i=0; i<n_dofs; ++i)
    if (matrix.diag_element(i) != 0)
    {
      first_nonzero_diagonal_entry = matrix.diag_element(i);
      break;
    };


  std::map<unsigned int,double>::const_iterator dof  = boundary_values.begin(),
    endd = boundary_values.end();
  const SparsityPattern    &sparsity    = matrix.get_sparsity_pattern();
  const unsigned int *sparsity_rowstart = (const unsigned int*) (sparsity.get_rowstart_indices());
  const unsigned int *sparsity_colnums  = (const unsigned int*) (sparsity.get_column_numbers());
  for (; dof != endd; ++dof)
  {
    Assert (dof->first < n_dofs, ExcInternalError());

    const unsigned int dof_number = dof->first;
    // for each boundary dof:

    // set entries of this line
    // to zero except for the diagonal
    // entry. Note that the diagonal
    // entry is always the first one
    // for square matrices, i.e.
    // we shall not set
    // matrix.global_entry(
    //     sparsity_rowstart[dof.first])
    const unsigned int last = sparsity_rowstart[dof_number+1];
    for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
      matrix.global_entry(j) = 0.;


    // set right hand side to
    // wanted value: if main diagonal
    // entry nonzero, don't touch it
    // and scale rhs accordingly. If
    // zero, take the first main
    // diagonal entry we can find, or
    // one if no nonzero main diagonal
    // element exists. Normally, however,
    // the main diagonal entry should
    // not be zero.
    //
    // store the new rhs entry to make
    // the gauss step more efficient
    double new_rhs;
    if (matrix.diag_element(dof_number) != 0.0)
    {
      new_rhs = dof->second * matrix.diag_element(dof_number);
      right_hand_side(dof_number) = new_rhs;
    }
    else
    {
      matrix.set (dof_number, dof_number,
                  first_nonzero_diagonal_entry);
      new_rhs = dof->second * first_nonzero_diagonal_entry;
      right_hand_side(dof_number) = new_rhs;
    }


      // store the only nonzero entry
      // of this line for the Gauss
      // elimination step
      const double diagonal_entry = matrix.diag_element(dof_number);

      // we have to loop over all
      // rows of the matrix which
      // have a nonzero entry in
      // the column which we work
      // in presently. if the
      // sparsity pattern is
      // symmetric, then we can
      // get the positions of
      // these rows cheaply by
      // looking at the nonzero
      // column numbers of the
      // present row. we need not
      // look at the first entry,
      // since that is the
      // diagonal element and
      // thus the present row
      for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
      {
        const unsigned int row = sparsity_colnums[j];

        // find the position of
        // element
        // (row,dof_number)
        const unsigned int *
          p = std::lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
                               &sparsity_colnums[sparsity_rowstart[row+1]],
                               dof_number);

        // check whether this line has
        // an entry in the regarding column
        // (check for ==dof_number and
        // != next_row, since if
        // row==dof_number-1, *p is a
        // past-the-end pointer but points
        // to dof_number anyway...)
        //
        // there should be such an entry!
        Assert ((*p == dof_number) &&
                (p != &sparsity_colnums[sparsity_rowstart[row+1]]),
                ExcInternalError());

        const unsigned int global_entry
          = (p - &sparsity_colnums[sparsity_rowstart[0]]);

        // correct right hand side
        right_hand_side(row) -= matrix.global_entry(global_entry) /
          diagonal_entry * new_rhs;

        // set matrix entry to zero
        matrix.global_entry(global_entry) = 0.;
      };

  };

}


void AssemblySystemBase::applySymetricBoundaryValue(MapIntDouble &boundary_values,BlockMatrix &matrix, BlockVecDouble &sol,BlockVecDouble &right_hand_side) 
{
  MatrixTools::apply_boundary_values(boundary_values,matrix,sol,right_hand_side,true);
}

---