poisson_prob.cpp

Go to the documentation of this file.

#include <stdlib.h>
#include <BasicFiniteElements.hpp>
#include <PoissonOperators.hpp>
 
using namespace MoFEM;
using namespace PoissonOps;
 
static char help[] = "...\n\n";
 
struct Poisson {
public:
  Poisson(moab::Core &mb_instance, MoFEM::Core &core, int order);
 
  MoFEMErrorCode runAnalysis();
 
private:
  MoFEMErrorCode readMesh();
  MoFEMErrorCode setupSystem();
  MoFEMErrorCode setIntegrationRule();
  MoFEMErrorCode assembleSystem();
  MoFEMErrorCode applyBc();
  MoFEMErrorCode solveSystem();
  MoFEMErrorCode setOutput();
 
  boost::shared_ptr<VolEle> volPipelineLhs;
  boost::shared_ptr<VolEle> volPipelineRhs;
  boost::shared_ptr<VolEle> postFaceVolume;
 
  MoFEM::Interface &mField;
  moab::Interface &moab;
  Simple *simpleInterface;
 
  DM dm;
  KSP ksp;
 
  MatrixDouble invJac;
 
      MPI_Comm mpiComm; // mpi parallel communicator
  const int mpiRank; // number of processors involved
 
  int order;
 
  std::string domainField;
};
 
Poisson::Poisson(moab::Core &mb_instance, MoFEM::Core &core, const int order)
    : domainField("U"), moab(mb_instance), mField(core),
      mpiComm(mField.get_comm()), mpiRank(mField.get_comm_rank()),
      order(order) {
  volPipelineLhs = boost::shared_ptr<VolEle>(new VolEle(mField));
  volPipelineRhs = boost::shared_ptr<VolEle>(new VolEle(mField));
}
 
MoFEMErrorCode Poisson::readMesh() {
  MoFEMFunctionBegin;
 
  CHKERR mField.getInterface(simpleInterface);
  CHKERR simpleInterface->getOptions();
  CHKERR simpleInterface->loadFile();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode Poisson::setupSystem() {
  MoFEMFunctionBegin;
 
  CHKERR simpleInterface->addDomainField(domainField, H1,
                                         AINSWORTH_LEGENDRE_BASE, 1);
  CHKERR simpleInterface->setFieldOrder(domainField, order);
  CHKERR simpleInterface->setUp();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode Poisson::setIntegrationRule() {
  MoFEMFunctionBegin;
  auto vol_ruleLhs = [](int, int, int p) -> int { return 2 * (p - 1); };
  auto vol_ruleRhs = [](int, int, int p) -> int { return 2 * (p - 1); };
  volPipelineLhs->getRuleHook = vol_ruleLhs;
  volPipelineRhs->getRuleHook = vol_ruleRhs;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode Poisson::assembleSystem() {
  MoFEMFunctionBegin;
  auto det_ptr = boost::make_shared<VectorDouble>();
  auto jac_ptr = boost::make_shared<MatrixDouble>();
  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
  volPipelineLhs->getOpPtrVector().push_back(
      new OpCalculateHOJacForFace(jac_ptr));
  volPipelineLhs->getOpPtrVector().push_back(
      new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
  volPipelineLhs->getOpPtrVector().push_back(
      new OpLhsUU(domainField, domainField));
  volPipelineRhs->getOpPtrVector().push_back(new OpRhsU(domainField));
 
  // dm = simpleInterface->getDM();
  CHKERR simpleInterface->getDM(&dm);
  boost::shared_ptr<VolEle> null;
 
  CHKERR setIntegrationRule();
 
  CHKERR DMMoFEMKSPSetComputeOperators(dm, simpleInterface->getDomainFEName(), volPipelineLhs, null, null);
  CHKERR DMMoFEMKSPSetComputeRHS(dm, simpleInterface->getDomainFEName(), volPipelineRhs, null, null);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode Poisson::solveSystem() {
  MoFEMFunctionBegin;
 
  Vec global_rhs, global_solution;
  CHKERR DMCreateGlobalVector(dm, &global_rhs);
  CHKERR VecDuplicate(global_rhs, &global_solution);
 
 
  // ksp = createKSP(mField.get_comm());
 
  CHKERR KSPCreate(PETSC_COMM_WORLD, &ksp);
 
  // CHKERR KSPSetOperators(ksp, volPipelineLhs->getFEMethod()->ksp_B,
  //                        volPipelineLhs->getFEMethod()->ksp_B);
 
  CHKERR KSPSetFromOptions(ksp);
 
  CHKERR KSPSetDM(ksp, dm);
 
  CHKERR KSPSetUp(ksp);
 
  CHKERR KSPSolve(ksp, global_rhs, global_solution);
 
  CHKERR DMoFEMMeshToGlobalVector(dm, global_solution, INSERT_VALUES,
                                  SCATTER_REVERSE);
 
  
  
  CHKERR KSPDestroy(&ksp);
  CHKERR VecDestroy(&global_solution);
  CHKERR VecDestroy(&global_rhs);
 
  CHKERR DMoFEMLoopFiniteElements(dm, simpleInterface->getDomainFEName(),
                                  postFaceVolume);
 
  CHKERR
  boost::static_pointer_cast<PostProcFaceOnRefinedMesh>(postFaceVolume)
      ->writeFile("out_vol.h5m");
 
  CHKERR DMDestroy(&dm);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode Poisson::setOutput() {
  MoFEMFunctionBegin;
 
  postFaceVolume =
      boost::shared_ptr<VolEle>(new PostProcFaceOnRefinedMesh(mField));
 
  CHKERR
  boost::static_pointer_cast<PostProcFaceOnRefinedMesh>(postFaceVolume)
      ->generateReferenceElementMesh();
  CHKERR
  boost::static_pointer_cast<PostProcFaceOnRefinedMesh>(postFaceVolume)
      ->addFieldValuesPostProc(domainField);
  // CHKERR
  // boost::static_pointer_cast<PostProcFaceOnRefinedMesh>(postFaceVolume)
  //     ->addFieldValuesGradientPostProc(domainField);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode Poisson::applyBc() {
  MoFEMFunctionBegin;
  Range bdry_entities;
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
    string name = it->getName();
    if (name.compare(0, 14, "ESSENTIAL") == 0) {
      CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
                                                 bdry_entities, true);
    }
  }
 
  cerr << bdry_entities;
  Range bdr_verts;
  CHKERR mField.get_moab().get_connectivity(bdry_entities, bdr_verts, true);
  bdry_entities.merge(bdr_verts);
 
  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
      simpleInterface->getProblemName(), domainField, bdry_entities);
 
  
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode Poisson::runAnalysis() {
  MoFEMFunctionBegin;
 
  CHKERR readMesh();
  CHKERR setupSystem();
  CHKERR applyBc();
  
  CHKERR assembleSystem();
  
  CHKERR setOutput();
  
  CHKERR solveSystem();
 
  MoFEMFunctionReturn(0);
}
 
int main(int argc, char *argv[]) {
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
  try {
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
    MoFEM::Core core(moab);
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    int order = 1;
    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
 
    Poisson poisson_problem(mb_instance, core, order);
 
    CHKERR poisson_problem.runAnalysis();
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}