solverAlgorithms.h

// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
//
// Contributor(s):
//   Eric Bechet
//

#ifndef _SOLVERALGORITHMS_H_
#define _SOLVERALGORITHMS_H_


#include "dofManager.h"
#include "terms.h"
#include "quadratureRules.h"
#include "MVertex.h"



template<class Iterator, class Assembler> void Assemble(BilinearTermBase &term, FunctionSpaceBase &space,
                                                        Iterator itbegin, Iterator itend,
                                                        QuadratureBase &integrator, Assembler &assembler)
  // symmetric
{
  fullMatrix<typename Assembler::dataMat> localMatrix;
  std::vector<Dof> R;
  for (Iterator it = itbegin; it != itend; ++it){
    MElement *e = *it;
    R.clear();
    IntPt *GP;
    int npts = integrator.getIntPoints(e, &GP);
    term.get(e, npts, GP, localMatrix); //localMatrix.print();
    space.getKeys(e, R);
    assembler.assemble(R, localMatrix);
  }
}

template<class Assembler> void Assemble(BilinearTermBase &term, FunctionSpaceBase &space, MElement *e,
                                        QuadratureBase &integrator, Assembler &assembler) // symmetric
{
  fullMatrix<typename Assembler::dataMat> localMatrix;
  std::vector<Dof> R;
  IntPt *GP;
  int npts = integrator.getIntPoints(e, &GP);
  term.get(e, npts, GP, localMatrix);
  space.getKeys(e, R);
  assembler.assemble(R, localMatrix);
}

template<class Iterator, class Assembler> void Assemble(BilinearTermBase &term,
                                                        FunctionSpaceBase &shapeFcts,
                                                        FunctionSpaceBase &testFcts,
                                                        Iterator itbegin, Iterator itend,
                                                        QuadratureBase &integrator,
                                                        Assembler &assembler) // non symmetric
{
  fullMatrix<typename Assembler::dataMat> localMatrix;
  std::vector<Dof> R;
  std::vector<Dof> C;
  for (Iterator it = itbegin; it != itend; ++it){
    MElement *e = *it;
    R.clear();
    C.clear();
    IntPt *GP;
    int npts = integrator.getIntPoints(e, &GP);
    term.get(e, npts, GP, localMatrix); //localMatrix.print();
    printf("local matrix size = %d %d\n", localMatrix.size1(), localMatrix.size2());
    shapeFcts.getKeys(e, R);
    testFcts.getKeys(e, C);