poissonfem.cpp

Go to the documentation of this file.

#include "poissonfem.h"
#include "assemblyorthomeshsystem.h"
#include "hdf5orthowriter.h"
#include "gnuplotanim.h"
#include "fconst3d.h"
#include "assemblyorthomeshsystemII.h"
#include "varform.h"
#include <fstream>


PoissonFEM::PoissonFEM(OrthoMesh &mesh,  VecDouble &cK, Function3D &fD, Function3D &fN,LinearSolver &solver, unsigned debugLevel) 
  :m_mesh(mesh), m_cK(cK),m_eqvar(ff),m_solver(solver),m_debugLevel(debugLevel), m_fN(fN)
{
  m_eqvar.reinit(mesh);
  m_eqvar.setParameters(cK);
  AssemblyOrthoMeshSystemII<VarFormPoisson> ass;
  std::cout << "Resp 2: " << ass.max_nz_per_row(m_eqvar) << std::endl;
  unsigned rank = ass.global_system_rank(m_eqvar);
  std::cout << "Resp 1: " << rank << std::endl;

  //Build the sparsity patter and initialize the sparse matrix
  ass.buildSparsityPattern(m_sparsity,m_eqvar);
  m_M.reinit(m_sparsity);
  m_vRHS.reinit(rank);
  m_sol.reinit(rank);
  //Print the sparsity
  //std::ofstream file("sparse.dat");
  //m_sparsity.print_gnuplot(file);

  ass.setBoundaryValue(m_mesh,m_eqvar,fD,m_dirichlet_bc);

  /*
  AssemblyOrthoMeshSystem<PoissonEqVar> assembler;
  m_vRHS.reinit(m_eqvar.n_dofs());
  m_sol.reinit(m_eqvar.n_dofs());
  m_sparsity.reinit(m_eqvar.n_dofs(),m_eqvar.n_dofs(),m_eqvar.n_max_dof_coupling());
  assembler.buildSparsityPattern(m_sparsity,m_mesh,m_eqvar);
  m_sparsity.compress();
  m_M.reinit(m_sparsity);
  assembler.setBoundaryValue(m_mesh,m_eqvar,fD,dirichlet_bc_map);
  m_eqvar.setSource(*(new FConst3D(0)));

  */
}


 void PoissonFEM::iterate(TransportBase &trans) 
{
  
  printf("Building system\n");
  AssemblyOrthoMeshSystemII<VarFormPoisson> assembler;
  assembler.buildGlobalMatrix(m_M, m_eqvar);
  assembler.buildGlobalRHS(m_vRHS,m_eqvar);

  assembler.addNeumannCondition(m_vRHS,m_eqvar,m_fN);
  //Apply boundary conditions
  m_vRHS.print(std::cout);
  assembler.applyBoundaryValue(m_dirichlet_bc,m_vRHS);
  assembler.applyBoundaryValue(m_dirichlet_bc,m_M);

  //solve the system
  printf("Solving\n");
  m_solver.solve(m_M,m_sol,m_vRHS);
  printf("done");
  
  //std::ofstream file("matrix.dat");
  //m_M.print_formatted(file,10,true,0," 0 ");
  //m_vRHS.print(std::cout);
  /*
  m_eqvar.setK(m_cK);
  assembler.buildGlobalRhs(m_vRHS,m_mesh,m_eqvar);
  assembler.applyBoundaryValue(dirichlet_bc_map,m_M);
  assembler.applyBoundaryValue(dirichlet_bc_map,m_vRHS);
  assembler.addNeummanCondition(m_mesh,m_vRHS,m_eqvar,m_fN);
  m_sol=m_vRHS;
  m_sol.print();
  solveA.factorize(m_M);
  solveA.solve(m_sol);
  */
}



 void PoissonFEM::printOutput() 
{
  AssemblyOrthoMeshSystemII<VarFormPoisson> assembler;
  VecDouble dofSol(m_sol.size());
  assembler.getDOFSolAtVertices(dofSol,m_sol,m_eqvar);
  HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 
  hdf5.writeScalarField(dofSol,"P");

  //dofSol-=m_sol;

  /*
  if (m_debugLevel > 0)
  {
    HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 
    hdf5.writeScalarField(m_sol,"P");

     GnuPlotAnim &plot = GnuPlotAnim::getGnuPlotAnim();
     
     //char str[100];
     //printf("Mass: %g\n",m_mesh.getIntegralAtCells(m_cValues));
     plot.set("yrange","[0:1]");
     plot.plotVerticeValuesAtFixedZ(m_mesh,m_sol,"Sc","with lines",5);
     plot.NextScene();
  }
  */
}


void PoissonFEM::getVelocitiesAtFaces(Matrix& M) 
{
  throw new Exception("PoissonFEM::getVelocitiesAtFaces Not implemented yet\n");
}

---