local_algebra.h

Go to the documentation of this file.

/*
 * Copyright (c) 2010-2015:  G-CSC, Goethe University Frankfurt
 * Author: Andreas Vogel
 * 
 * This file is part of UG4.
 * 
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 * 
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 */
 
#ifndef __H__UG__LIB_DISC__COMMON__LOCAL_ALGEBRA__
#define __H__UG__LIB_DISC__COMMON__LOCAL_ALGEBRA__
 
//#define UG_LOCALALGEBRA_ASSERT(cond, exp)
// include define below to assert arrays used in stabilization
#define UG_LOCALALGEBRA_ASSERT(cond, exp) UG_ASSERT((cond), exp)
 
#include <vector>
 
#include "./multi_index.h"
#include "./function_group.h"
#include "lib_algebra/small_algebra/small_algebra.h"
 
namespace ug{
 
 
class LocalIndices
{
  public:
    typedef size_t index_type;
 
    typedef index_type comp_type;
 
  public:
    LocalIndices() {};
 
    void resize_fct(size_t numFct)
    {
      m_vvIndex.resize(numFct);
      m_vLFEID.resize(numFct);
    }
 
    void set_lfeID(size_t fct, const LFEID& lfeID)
    {
      UG_ASSERT(fct < m_vLFEID.size(), "Invalid index: "<<fct);
      m_vLFEID[fct] = lfeID;
    }
 
    const LFEID& local_finite_element_id(size_t fct) const
    {
      UG_ASSERT(fct < m_vLFEID.size(), "Invalid index: "<<fct);
      return m_vLFEID[fct];
    }
 
    void resize_dof(size_t fct, size_t numDoF)
    {
      check_fct(fct);
      m_vvIndex[fct].resize(numDoF);
    }
 
    void clear_dof(size_t fct) {resize_dof(fct, 0);}
 
    void reserve_dof(size_t fct, size_t numDoF)
    {
      check_fct(fct);
      m_vvIndex[fct].reserve(numDoF);
    }
 
    void push_back_index(size_t fct, size_t index) {push_back_multi_index(fct,index,0);}
 
    void push_back_multi_index(size_t fct, size_t index, size_t comp)
    {
      check_fct(fct);
      m_vvIndex[fct].push_back(DoFIndex(index,comp));
    }
 
    void clear() {m_vvIndex.clear();}
 
    size_t num_fct() const {return m_vvIndex.size();}
 
    size_t num_dof(size_t fct) const
    {
      check_fct(fct);
      return m_vvIndex[fct].size();
    }
 
    size_t num_dof() const
    {
      size_t num = 0;
      for(size_t fct = 0; fct < num_fct(); ++fct) num += num_dof(fct);
      return num;
    }
 
    const DoFIndex& multi_index(size_t fct, size_t dof) const
    {
      check_dof(fct, dof);
      return m_vvIndex[fct][dof];
    }
 
    index_type index(size_t fct, size_t dof) const
    {
      check_dof(fct, dof);
      return m_vvIndex[fct][dof][0];
    }
 
    index_type& index(size_t fct, size_t dof)
    {
      check_dof(fct, dof);
      return m_vvIndex[fct][dof][0];
    }
 
    comp_type comp(size_t fct, size_t dof) const
    {
      check_dof(fct, dof);
      return m_vvIndex[fct][dof][1];
    }
 
    comp_type& comp(size_t fct, size_t dof)
    {
      check_dof(fct, dof);
      return m_vvIndex[fct][dof][1];
    }
  
    bool contains_index(index_type ind)
    {
      for(size_t fct = 0; fct < num_fct(); fct++)
        for(size_t dof = 0; dof < num_dof(fct); dof++)
          if(m_vvIndex[fct][dof][0] == ind)
            return true;
      return false;
    }
 
  protected:
    inline void check_fct(size_t fct) const
    {
      UG_LOCALALGEBRA_ASSERT(fct < num_fct(), "Wrong index.");
    }
    inline void check_dof(size_t fct, size_t dof) const
    {
      check_fct(fct);
      UG_LOCALALGEBRA_ASSERT(dof < num_dof(fct), "Wrong index.");
    }
 
  protected:
  //  Mapping (fct, dof) -> local index
    std::vector<std::vector<DoFIndex> > m_vvIndex;
 
  //  Local finite element ids
    std::vector<LFEID> m_vLFEID;
};
 
class LocalVector
{
  public:
    typedef LocalVector this_type;
 
    typedef number value_type;
 
  public:
    LocalVector() : m_pIndex(NULL) {m_vvValue.clear();}
 
    LocalVector(const LocalIndices& ind) {resize(ind);}
 
    void resize(const LocalIndices& ind)
    {
      m_pIndex = &ind;
      m_vvValue.resize(ind.num_fct());
      m_vvValueAcc.resize(m_vvValue.size());
      for(size_t fct = 0; fct < m_vvValue.size(); ++fct)
        m_vvValue[fct].resize(ind.num_dof(fct));
      access_all();
    }
 
    const LocalIndices& get_indices() const {return *m_pIndex;}
 
    // vector functions
 
  this_type& operator=(const this_type& other)
  {
    m_pIndex = other.m_pIndex;
    const size_t numFcts = m_pIndex->num_fct();
    m_vvValue.resize(numFcts);
    for (size_t fct = 0; fct < numFcts; ++fct)
      m_vvValue[fct].resize(m_pIndex->num_dof(fct));
    m_vvValueAcc.resize(numFcts);
    if (other.m_pFuncMap)
      access_by_map(*other.m_pFuncMap);
    else
      access_all();
 
    return *this;
  }
 
    this_type& operator=(number val)
    {
      for(size_t fct = 0; fct < m_vvValue.size(); ++fct)
        for(size_t dof = 0; dof < m_vvValue[fct].size(); ++dof)
          m_vvValue[fct][dof] = val;
      return *this;
    }
 
    this_type& operator*(number val)
    {
      return this->operator*=(val);
    }
 
    this_type& operator*=(number val)
    {
      for(size_t fct = 0; fct < m_vvValue.size(); ++fct)
        for(size_t dof = 0; dof < m_vvValue[fct].size(); ++dof)
          m_vvValue[fct][dof] *= val;
      return *this;
    }
 
    this_type& operator+=(const this_type& rhs)
    {
      UG_LOCALALGEBRA_ASSERT(m_pIndex==rhs.m_pIndex, "Not same indices.");
      for(size_t fct = 0; fct < m_vvValue.size(); ++fct)
        for(size_t dof = 0; dof < m_vvValue[fct].size(); ++dof)
          m_vvValue[fct][dof] += rhs.m_vvValue[fct][dof];
      return *this;
    }
 
    this_type& operator-=(const this_type& rhs)
    {
      UG_LOCALALGEBRA_ASSERT(m_pIndex==rhs.m_pIndex, "Not same indices.");
      for(size_t fct = 0; fct < m_vvValue.size(); ++fct)
        for(size_t dof = 0; dof < m_vvValue[fct].size(); ++dof)
          m_vvValue[fct][dof] -= rhs.m_vvValue[fct][dof];
      return *this;
    }
 
    this_type& scale_append(number s, const this_type& rhs)
    {
      UG_LOCALALGEBRA_ASSERT(m_pIndex==rhs.m_pIndex, "Not same indices.");
      for(size_t fct = 0; fct < m_vvValue.size(); ++fct)
        for(size_t dof = 0; dof < m_vvValue[fct].size(); ++dof)
          m_vvValue[fct][dof] += s * rhs.m_vvValue[fct][dof];
      return *this;
    }
 
    // restricted DoF access
 
    void access_by_map(const FunctionIndexMapping& funcMap)
    {
      m_pFuncMap = &funcMap;
      for(size_t i = 0; i < funcMap.num_fct(); ++i)
      {
        const size_t mapFct = funcMap[i];
        m_vvValueAcc[i] = &(m_vvValue[mapFct][0]);
      }
    }
 
    void access_all()
    {
      m_pFuncMap = NULL;
 
      if(m_pIndex==NULL) {m_vvValueAcc.clear(); return;}
 
      for(size_t i = 0; i < m_pIndex->num_fct(); ++i)
        m_vvValueAcc[i] = &(m_vvValue[i][0]);
    }
 
    size_t num_fct() const
    {
      if(m_pFuncMap == NULL) return m_vvValue.size();
      return m_pFuncMap->num_fct();
    }
 
    const LFEID& local_finite_element_id(size_t fct) const
    {
      UG_ASSERT(m_pIndex != NULL, "No indices present");
      check_fct(fct);
      if(m_pFuncMap == NULL) return m_pIndex->local_finite_element_id(fct);
      else return m_pIndex->local_finite_element_id((*m_pFuncMap)[fct]);
    }
 
    size_t num_dof(size_t fct) const
    {
      check_fct(fct);
      if(m_pFuncMap == NULL) return m_vvValue[fct].size();
      else return m_vvValue[ (*m_pFuncMap)[fct] ].size();
    }
 
    number& operator()(size_t fct, size_t dof)
    {
      check_dof(fct,dof);
      return m_vvValueAcc[fct][dof];
    }
 
    number operator()(size_t fct, size_t dof) const
    {
      check_dof(fct,dof);
      return m_vvValueAcc[fct][dof];
    }
 
    // all DoF access
 
    size_t num_all_fct() const {return m_vvValue.size();}
 
    size_t num_all_dof(size_t fct) const {check_all_fct(fct); return m_vvValue[fct].size();}
 
    number& value(size_t fct, size_t dof){check_all_dof(fct,dof);return m_vvValue[fct][dof];}
 
    const number& value(size_t fct, size_t dof) const{check_all_dof(fct,dof);return m_vvValue[fct][dof];}
 
  protected:
    inline void check_fct(size_t fct) const
    {
      UG_LOCALALGEBRA_ASSERT(fct < num_fct(), "Wrong index: fct: "<<fct<<
                  " of num_fct: "<<num_fct());
    }
    inline void check_dof(size_t fct, size_t dof) const
    {
      check_fct(fct);
      UG_LOCALALGEBRA_ASSERT(dof < num_dof(fct), "Wrong index: dof: "
          <<dof<<", num_dof: "<<num_dof(fct)<<" of fct: "<<fct);
    }
    inline void check_all_fct(size_t fct) const
    {
      UG_LOCALALGEBRA_ASSERT(fct < num_all_fct(), "Wrong index.");
    }
    inline void check_all_dof(size_t fct, size_t dof) const
    {
      check_all_fct(fct);
      UG_LOCALALGEBRA_ASSERT(dof < num_all_dof(fct), "Wrong index.");
    }
 
  protected:
    const LocalIndices* m_pIndex;
 
    const FunctionIndexMapping* m_pFuncMap;
 
    std::vector<value_type*> m_vvValueAcc;
 
    std::vector<std::vector<value_type> > m_vvValue;
};
 
class LocalMatrix
{
  public:
    typedef LocalMatrix this_type;
 
    typedef number value_type;
 
  public:
    LocalMatrix() :
      m_pRowIndex(NULL), m_pColIndex(NULL) ,
      m_pRowFuncMap(NULL), m_pColFuncMap(NULL)
    {}
 
    LocalMatrix(const LocalIndices& rowInd, const LocalIndices& colInd)
      : m_pRowFuncMap(NULL), m_pColFuncMap(NULL)
    {
      resize(rowInd, colInd);
    }
 
    void resize(const LocalIndices& ind) {resize(ind, ind);}
 
    void resize(const LocalIndices& rowInd, const LocalIndices& colInd)
    {
      m_pRowIndex = &rowInd;
      m_pColIndex = &colInd;
 
      m_CplMat.resize(rowInd.num_fct(), colInd.num_fct());
      m_CplMatAcc.resize(rowInd.num_fct(), colInd.num_fct());
 
      for(size_t fct1 = 0; fct1 < m_CplMat.num_rows(); ++fct1)
        for(size_t fct2 = 0; fct2 < m_CplMat.num_cols(); ++fct2)
        {
          m_CplMat(fct1, fct2).resize(colInd.num_dof(fct1),
                        rowInd.num_dof(fct2));
        }
 
      access_all();
    }
 
    const LocalIndices& get_row_indices() const {return *m_pRowIndex;}
 
    const LocalIndices& get_col_indices() const {return *m_pColIndex;}
 
    // Matrix functions
 
    this_type& operator=(number val)
    {
      for(size_t fct1=0; fct1 < m_CplMat.num_rows(); ++fct1)
        for(size_t fct2=0; fct2 < m_CplMat.num_cols(); ++fct2)
          m_CplMat(fct1,fct2) = val;
      return *this;
    }
 
    this_type& operator*(number val)
    {
      return this->operator*=(val);
    }
 
    this_type& operator*=(number val)
    {
      for(size_t fct1=0; fct1 < m_CplMat.num_rows(); ++fct1)
        for(size_t fct2=0; fct2 < m_CplMat.num_cols(); ++fct2)
          m_CplMat(fct1,fct2) *= val;
      return *this;
    }
 
    this_type& operator+=(const this_type& rhs)
    {
      UG_LOCALALGEBRA_ASSERT(m_pRowIndex==rhs.m_pRowIndex &&
                m_pColIndex==rhs.m_pColIndex, "Not same indices.");
      for(size_t fct1=0; fct1 < m_CplMat.num_rows(); ++fct1)
        for(size_t fct2=0; fct2 < m_CplMat.num_cols(); ++fct2)
          m_CplMat(fct1,fct2) += rhs.m_CplMat(fct1,fct2);
      return *this;
    }
 
    this_type& operator-=(const this_type& rhs)
    {
      UG_LOCALALGEBRA_ASSERT(m_pRowIndex==rhs.m_pRowIndex &&
                m_pColIndex==rhs.m_pColIndex, "Not same indices.");
      for(size_t fct1=0; fct1 < m_CplMat.num_rows(); ++fct1)
        for(size_t fct2=0; fct2 < m_CplMat.num_cols(); ++fct2)
          m_CplMat(fct1,fct2) -= rhs.m_CplMat(fct1,fct2);
      return *this;
    }
 
    this_type& scale_append(number s, const this_type& rhs)
    {
      UG_LOCALALGEBRA_ASSERT(m_pRowIndex==rhs.m_pRowIndex &&
            m_pColIndex==rhs.m_pColIndex, "Not same indices.");
      for(size_t fct1=0; fct1 < m_CplMat.num_rows(); ++fct1)
        for(size_t fct2=0; fct2 < m_CplMat.num_cols(); ++fct2)
          MatScaleAppend(m_CplMat(fct1,fct2), s, rhs.m_CplMat(fct1,fct2));
      return *this;
    }
 
    // restricted DoF access
 
    void access_by_map(const FunctionIndexMapping& funcMap)
    {
      access_by_map(funcMap, funcMap);
    }
 
    void access_by_map(const FunctionIndexMapping& rowFuncMap,
                       const FunctionIndexMapping& colFuncMap)
    {
      m_pRowFuncMap = &rowFuncMap;
      m_pColFuncMap = &colFuncMap;
 
      for(size_t i = 0; i < m_pRowFuncMap->num_fct(); ++i)
        for(size_t j = 0; j < m_pColFuncMap->num_fct(); ++j)
        {
          const size_t rowMapFct = rowFuncMap[i];
          const size_t colMapFct = colFuncMap[j];
 
          m_CplMatAcc(i,j) = &(m_CplMat(rowMapFct, colMapFct));
        }
    }
 
    void access_all()
    {
      m_pRowFuncMap = NULL;
      m_pColFuncMap = NULL;
 
      if(m_pRowIndex==NULL) {m_CplMatAcc.resize(0,0); return;}
 
      for(size_t i = 0; i < m_pRowIndex->num_fct(); ++i)
        for(size_t j = 0; j < m_pColIndex->num_fct(); ++j)
          m_CplMatAcc(i,j) = &(m_CplMat(i,j));
    }
 
    size_t num_row_fct() const
    {
      if(m_pRowFuncMap != NULL) return m_pRowFuncMap->num_fct();
      return m_CplMatAcc.num_rows();
    }
 
    size_t num_col_fct() const
    {
      if(m_pColFuncMap != NULL) return m_pColFuncMap->num_fct();
      return m_CplMatAcc.num_cols();
    }
 
    size_t num_row_dof(size_t fct) const
    {
      if(m_CplMat.num_rows()==0) return 0;
      if(m_pRowFuncMap == NULL) return m_CplMat(fct, 0).num_rows();
      else return m_CplMat( (*m_pRowFuncMap)[fct], 0).num_rows();
    }
 
    size_t num_col_dof(size_t fct) const
    {
      if(m_CplMat.num_cols()==0) return 0;
      if(m_pRowFuncMap == NULL) return m_CplMat(0, fct).num_cols();
      else return m_CplMat( 0, (*m_pColFuncMap)[fct]).num_cols();
    }
 
    number& operator()(size_t rowFct, size_t rowDoF,
                       size_t colFct, size_t colDoF)
    {
      check_dof(rowFct, rowDoF, colFct, colDoF);
      return (*m_CplMatAcc(rowFct,colFct))(rowDoF, colDoF);
    }
 
    number operator()(size_t rowFct, size_t rowDoF,
                            size_t colFct, size_t colDoF) const
    {
      check_dof(rowFct, rowDoF, colFct, colDoF);
      return (*m_CplMatAcc(rowFct,colFct))(rowDoF, colDoF);
    }
 
    // all DoF access
 
    size_t num_all_row_fct() const{ return m_CplMat.num_rows();}
 
    size_t num_all_col_fct() const{return m_CplMat.num_cols();}
 
    size_t num_all_row_dof(size_t fct) const {return m_CplMat(fct, 0).num_rows();}
 
    size_t num_all_col_dof(size_t fct) const {return m_CplMat(0, fct).num_cols();}
 
    number& value(size_t rowFct, size_t rowDoF,
                  size_t colFct, size_t colDoF)
    {
      check_all_dof(rowFct, rowDoF, colFct, colDoF);
      return (m_CplMat(rowFct,colFct))(rowDoF, colDoF);
    }
 
    number value(size_t rowFct, size_t rowDoF,
                       size_t colFct, size_t colDoF) const
    {
      check_all_dof(rowFct, rowDoF, colFct, colDoF);
      return (m_CplMat(rowFct,colFct))(rowDoF, colDoF);
    }
 
  protected:
    inline void check_fct(size_t rowFct, size_t colFct) const
    {
      UG_LOCALALGEBRA_ASSERT(rowFct < num_row_fct(), "Wrong index.");
      UG_LOCALALGEBRA_ASSERT(colFct < num_col_fct(), "Wrong index.");
    }
    inline void check_dof(size_t rowFct, size_t rowDoF,
                          size_t colFct, size_t colDoF) const
    {
      check_fct(rowFct, colFct);
      UG_LOCALALGEBRA_ASSERT(rowDoF < num_row_dof(rowFct), "Wrong index.");
      UG_LOCALALGEBRA_ASSERT(colDoF < num_col_dof(colFct), "Wrong index.");
    }
    inline void check_all_fct(size_t rowFct, size_t colFct) const
    {
      UG_LOCALALGEBRA_ASSERT(rowFct < num_all_row_fct(), "Wrong index.");
      UG_LOCALALGEBRA_ASSERT(colFct < num_all_col_fct(), "Wrong index.");
    }
    inline void check_all_dof(size_t rowFct, size_t rowDoF,
                size_t colFct, size_t colDoF) const
    {
      check_all_fct(rowFct, colFct);
      UG_LOCALALGEBRA_ASSERT(rowDoF < num_all_row_dof(rowFct), "Wrong index.");
      UG_LOCALALGEBRA_ASSERT(colDoF < num_all_col_dof(colFct), "Wrong index.");
    }
 
  protected:
  //  Row indices
    const LocalIndices* m_pRowIndex;
 
  //  Column indices
    const LocalIndices* m_pColIndex;
 
    const FunctionIndexMapping* m_pRowFuncMap;
 
    const FunctionIndexMapping* m_pColFuncMap;
 
  //  \todo: Think of a better (faster) storage (one plain array?)
 
  //  type of cpl matrices between two functions
    typedef DenseMatrix<VariableArray2<value_type> > LocalCplMatrix;
 
  //  type of Func-Coupling matrices
    typedef DenseMatrix<VariableArray2<LocalCplMatrix> > FctCplMatrix;
 
  //  type of Func-Coupling pointer matrices
    typedef DenseMatrix<VariableArray2<LocalCplMatrix*> > FctCplAccMatrix;
 
  //  Entries (fct1, fct2, dof1, dof2)
    FctCplMatrix m_CplMat;
 
  //  Entries (fct1, fct2, dof1, dof2)
    FctCplAccMatrix m_CplMatAcc;
};
 
inline
std::ostream& operator<< (std::ostream& outStream, const ug::LocalMatrix& mat)
{
  for(size_t fct1 = 0; fct1 < mat.num_row_fct(); ++fct1)
    for(size_t fct2 = 0; fct2 < mat.num_col_fct(); ++fct2)
      for(size_t dof1 = 0; dof1 < mat.num_row_dof(fct1); ++dof1)
        for(size_t dof2 = 0; dof2 < mat.num_col_dof(fct2); ++dof2)
        {
          outStream << "[("<< fct1 << "," << dof1 << ") x ("
                       << fct2 << "," << dof2 << ")]: "
                       << mat(fct1, dof1, fct2, dof2) << "\n";
        }
  return outStream;
}
 
inline
std::ostream& operator<< (std::ostream& outStream, const ug::LocalVector& vec)
{
  for(size_t fct = 0; fct < vec.num_fct(); ++fct)
    for(size_t dof = 0; dof < vec.num_dof(fct); ++dof)
    {
      outStream << "["<<fct<<","<<dof<<"]:" << vec(fct, dof) << "\n";
    }
  return outStream;
}
 
template <typename TVector>
void GetLocalVector(LocalVector& lvec, const TVector& vec)
{
  const LocalIndices& ind = lvec.get_indices();
 
  for(size_t fct=0; fct < lvec.num_all_fct(); ++fct)
    for(size_t dof=0; dof < lvec.num_all_dof(fct); ++dof)
    {
      const size_t index = ind.index(fct,dof);
      const size_t comp = ind.comp(fct,dof);
      lvec.value(fct,dof) = BlockRef(vec[index], comp);
    }
}
 
template <typename TVector>
void AddLocalVector(TVector& vec, const LocalVector& lvec)
{
  const LocalIndices& ind = lvec.get_indices();
 
  for(size_t fct=0; fct < lvec.num_all_fct(); ++fct)
    for(size_t dof=0; dof < lvec.num_all_dof(fct); ++dof)
    {
      const size_t index = ind.index(fct,dof);
      const size_t comp = ind.comp(fct,dof);
      BlockRef(vec[index], comp) += lvec.value(fct,dof);
    }
}
 
template <typename TMatrix>
void AddLocalMatrixToGlobal(TMatrix& mat, const LocalMatrix& lmat)
{
  //PROFILE_FUNC_GROUP("algebra")
  const LocalIndices& rowInd = lmat.get_row_indices();
  const LocalIndices& colInd = lmat.get_col_indices();
 
  for(size_t fct1=0; fct1 < lmat.num_all_row_fct(); ++fct1)
    for(size_t dof1=0; dof1 < lmat.num_all_row_dof(fct1); ++dof1)
    {
      const size_t rowIndex = rowInd.index(fct1,dof1);
      const size_t rowComp = rowInd.comp(fct1,dof1);
 
      for(size_t fct2=0; fct2 < lmat.num_all_col_fct(); ++fct2)
        for(size_t dof2=0; dof2 < lmat.num_all_col_dof(fct2); ++dof2)
        {
          const size_t colIndex = colInd.index(fct2,dof2);
          const size_t colComp = colInd.comp(fct2,dof2);
 
          BlockRef(mat(rowIndex, colIndex), rowComp, colComp)
                += lmat.value(fct1,dof1,fct2,dof2);
        }
    }
}
 
} // end namespace ug
 
#endif /* __H__UG__LIB_DISC__COMMON__LOCAL_ALGEBRA__ */