elasticfem.cpp

Go to the documentation of this file.

#include "elasticfem.h"
#include "assemblyorthomeshsystem.h"
#include "hdf5orthowriter.h"
#include "gnuplotanim.h"


ElasticFEM::ElasticFEM(OrthoMesh &mesh,  VecDouble &cK,double &fYoung, double &fPoisson,Function3D &fD, Function3D &fN) 
  :m_mesh(mesh),m_cK(cK),m_eqvar(mesh,cK,fYoung,fPoisson,fD,fN),m_fN(fN)
 {
  AssemblyOrthoMeshSystem<IsotropicElasticEqVar> 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);

}



 void ElasticFEM::iterate(TransportBase &trans) 
{
  AssemblyOrthoMeshSystem<IsotropicElasticEqVar> assembler;
  assembler.buildGlobalMatrix(m_M,m_mesh, m_eqvar);
  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 ElasticFEM::printOutput() 
{
  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 ElasticFEM::getVelocitiesAtFaces(Matrix& M) 
{
  throw new Exception("ElasticFEM::getVelocitiesAtFaces Not implemented for this class\n");
}

---