An Object Finite Element Library

Assembly.h

Go to the documentation of this file.

/*==============================================================================

                                    O  F  E  L  I

                            Object  Finite  Element  Library

  ==============================================================================

   Copyright (C) 1998 - 2008 Rachid Touzani

   This program is free software; you can redistribute it and/or modify it under
   the terms of the GNU General Public License as published by the Free
   Software Foundation; Version 2 of the License.

   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 General Public License for more
   details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the :

   Free Software Foundation
   Inc., 59 Temple Place - Suite 330
   Boston, MA  02111-1307, USA

  ==============================================================================

                       Functions for Finite Element Assembly

  ==============================================================================*/

#ifndef __ASSEMBLY_H
#define __ASSEMBLY_H

#include "OFELI_Config.h"
#include "Vect.h"
#include "LocalVect.h"
#include "LocalMatrix.h"
#include "SkMatrix.h"
#include "SkSMatrix.h"
#include "SpMatrix.h"

namespace OFELI {

template<class T_,size_t N_,class E_>
void Assembly(const E_ &e, const LocalVect<T_,N_> &be, Vect<T_> &b)
{
   size_t i=1;
   for (size_t n=1; n<=e.getNbNodes(); n++) {
      Node *nd = e.getPtrNode(n);
      for (size_t k=1; k<=nd->getNbDOF(); k++)
         b(nd->getDOF(k)) += be(i++);
   }
}


template<class T_,size_t N_,class E_>
void Assembly(const E_ &e, const LocalMatrix<T_,N_,N_> &ae, Vect<T_> &b)
{
   size_t i = 0;
   for (size_t n=1; n<=e.getNbNodes(); n++) {
      Node *nd = e.getPtrNode(n);
      for (size_t k=1; k<=nd->getNbDOF(); k++) {
         i++;
         b(nd->getDOF(k)) += ae(i,i);
      }
   }
}


template<class T_,size_t N_,class E_>
void Assembly(const E_ &e, const LocalMatrix<T_,N_,N_> &ae, SkMatrix<T_> &a)
{
   size_t i = 1;
   for (size_t in=1; in<=e.getNbNodes(); in++) {
      Node *nd1 = e.getPtrNode(in);
      for (size_t k=1; k<=nd1->getNbDOF(); k++) {
         size_t ii = nd1->getDOF(k);
         size_t j = 1;
         for (size_t jn=1; jn<=e.getNbNodes(); jn++) {
            Node *nd2 = e.getPtrNode(jn);
            for (size_t l=1; l<=nd2->getNbDOF(); l++) {
               a.Add(ii,nd2->getDOF(l),ae(i,j));
               j++;
            }
         }
         i++;
      }
   }
}


template<class T_,size_t N_,class E_>
void Assembly(const E_ &e, const LocalMatrix<T_,N_,N_> &ae, SkSMatrix<T_> &a)
{
   size_t i=1, ii, jj;
   for (size_t in=1; in<=e.getNbNodes(); in++) {
      Node *nd1 = e.getPtrNode(in);
      for (size_t k=1; k<=nd1->getNbDOF(); k++) {
         if ((ii=nd1->getDOF(k)) > 0) {
           size_t j = 1;
           for (size_t jn=1; jn<=e.getNbNodes(); jn++) {
              Node *nd2 = e.getPtrNode(jn);
              for (size_t l=1; l<=nd2->getNbDOF(); l++) {
                 if ((jj=nd2->getDOF(l)) > 0)
                   a.Add(ii,jj,ae(i,j));
                 j++;
              }
           }
         }
         i++;
      }
   }
}


template<class T_,size_t N_,class E_>
void Assembly(const E_ &e, const LocalMatrix<T_,N_,N_> &ae, SpMatrix<T_> &a)
{
   size_t i = 1;
   for (size_t in=1; in<=e.getNbNodes(); in++) {
      Node *nd1 = e.getPtrNode(in);
      for (size_t k=1; k<=nd1->getNbDOF(); k++) {
         size_t ii=nd1->getDOF(k), j=1;
         for (size_t jn=1; jn<=e.getNbNodes(); jn++) {
            Node *nd2 = e.getPtrNode(jn);
            for (size_t l=1; l<=nd2->getNbDOF(); l++) {
               size_t jj = nd2->getDOF(l);
               a.Add(ii,jj,ae(i,j));
               j++;
            }
         }
         i++;
      }
   }
}


#ifndef DOXYGEN_SHOULD_SKIP_THIS

template<class T_,size_t N_,class E_>
void SAssembly(const E_ &e, const LocalVect<T_,N_> &be, Vect<T_> &b)
{
   size_t i = 1;
   for (size_t n=1; n<=N_; n++)
      b(e.getNodeLabel(n)) += be(i++);
}


template<class T_,size_t N_,class E_>
void SAssembly(const E_ &e, const LocalMatrix<T_,N_,N_> &ae, Vect<T_> &b)
{
   for (size_t n=1; n<=N_; n++)
      b(e.getNodeLabel(n)) += ae(n,n);
}


template<class T_,size_t N_,class E_>
void SAssembly(const E_ &e, const LocalMatrix<T_,N_,N_> &ae, SkSMatrix<T_> &a)
{
   size_t i = 1;
   for (size_t in=1; in<=N_; in++) {
      size_t ii = e.getPtrNode(in)->getLabel();
      size_t j = 1;
      for (size_t jn=1; jn<=N_; jn++) {
         size_t jj = e.getNodeLabel(jn);
         if (ii >= jj)
           a.Add(ii,jj,ae(i,j));
         j++;
      }
      i++;
   }
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


template<class T_,size_t N_>
void SideAssembly(const Element &e, const LocalMatrix<T_,N_,N_> &ae, SpMatrix<T_> &a)
{
   size_t i = 1;
   for (size_t n=1; n<=e.getNbSides(); n++) {
      Side *s1 = e.getPtrSide(n);
      for (size_t k=1; k<=s1->getNbDOF(); k++) {
         if (s1->getDOF(k)!=0) {
            size_t j = 1;
            for (size_t m=1; m<=e.getNbSides(); m++) {
               Side *s2 = e.getPtrSide(m);
               for (size_t l=1; l<=s2->getNbDOF(); l++) {
                  if (s2->getDOF(l)!=0)
                     a.Add(s1->getDOF(k),s2->getDOF(l),ae(i,j));
                  j++;
               }
            }
         }
         i++;
      }
   }
}


template<class T_,size_t N_>
void SideAssembly(const Element &e, const LocalMatrix<T_,N_,N_> &ae, SkSMatrix<T_> &a)
{
   size_t i=1, j, ii, jj;
   for (size_t n=1; n<=e.getNbSides(); n++) {
      Side *s1 = e.getPtrSide(n);
      for (size_t k=1; k<=s1->getNbDOF(); k++) {
         if ((ii=s1->getDOF(k)) > 0) {
            j = 1;
            for (size_t m=1; m<=e.getNbSides(); m++) {
               Side *s2 = e.getPtrSide(m);
               for (size_t l=1; l<=s2->getNbDOF(); l++) {
                  if ((jj=s2->getDOF(l)) > 0)
                     a.Add(ii,jj,ae(i,j));
                  j++;
               }
            }
         }
         i++;
      }
   }
}


template<class T_,size_t N_>
void SideAssembly(const Element &e, const LocalMatrix<T_,N_,N_> &ae, SkMatrix<T_> &a)
{
   size_t i = 1;
   for (size_t n=1; n<=e.getNbSides(); n++) {
      Side *s1 = e.getPtrSide(n);
      for (size_t k=1; k<=s1->getNbDOF(); k++) {
         size_t ii = s1->getDOF(k);
         if (ii != 0) {
            size_t j = 1;
            for (size_t m=1; m<=e.getNbSides(); m++) {
               Side *s2 = e.getPtrSide(m);
               for (size_t l=1; l<=s2->getNbDOF(); l++) {
                  size_t jj = s2->getDOF(l);
                  if (jj)
                     a.Add(ii,jj,ae(i,j));
                  j++;
               }
            }
         }
         i++;
      }
   }
}


template<class T_,size_t N_>
void SideAssembly(const Element &e, const LocalVect<T_,N_> &be, Vect<T_> &b)
{
   size_t i = 1;
   for (size_t n=1; n<=e.getNbSides(); n++) {
      Side *sd = e.getPtrSide(n);
      for (size_t k=1; k<=sd->getNbDOF(); k++) {
         if (sd->getDOF(k) != 0)
            b(sd->getDOF(k)) += be(i);
         i++;
      }
   }
}


/*template<class T_,size_t N_>
void SideAssembly(const Side &s, const LocalVect<T_,N_> &be, Vect<T_> &b)
{
   size_t i = 1;
   for (size_t n=1; n<=s.getNbNodes(); n++) {
      Side *sd = s.getPtrSide(n);
      for (size_t k=1; k<=sd->getNbDOF(); k++) {
         if (sd->getDOF(k) != 0)
            b(sd->getDOF(k)) += be(i);
         i++;
      }
   }
}*/

/*  @} */

} /* namespace OFELI */

#endif

Copyright © 1998-2007 Rachid Touzani