h_adaptive_transport.cpp File Reference

Example implementation of transport problem using mixed formulation. More...

#include <BasicFiniteElements.hpp>
#include <MixTransportElement.hpp>

Go to the source code of this file.

Classes
struct	BcFluxData
struct	MyTransport
	Application of mix transport data structure. More...

Typedefs
typedef map< int, BcFluxData >	BcFluxMap

Functions
int	main (int argc, char *argv[])

Variables
static char	help []

Detailed Description

Example implementation of transport problem using mixed formulation.

Todo:

Should be implemented and tested problem from this article Demkowicz, Leszek, and Jayadeep Gopalakrishnan. "Analysis of the DPG method for the Poisson equation." SIAM Journal on Numerical Analysis 49.5 (2011): 1788-1809.

Definition in file h_adaptive_transport.cpp.

Typedef Documentation

BcFluxMap

typedef map<int, BcFluxData> BcFluxMap

Definition at line 46 of file h_adaptive_transport.cpp.

Function Documentation

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 375 of file h_adaptive_transport.cpp.

                                 {

  const string default_options = "-ksp_type fgmres \n"
                                 "-pc_type lu \n"
                                 "-pc_factor_mat_solver_package mumps \n"
                                 "-ksp_monitor\n";

  string param_file = "param_file.petsc";
  if (!static_cast<bool>(ifstream(param_file))) {
    std::ofstream file(param_file.c_str(), std::ios::ate);
    if (file.is_open()) {
      file << default_options;
      file.close();
    }
  }

  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);

  try {

    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
    int rank;
    MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

    // get file name form command line
    PetscBool flg = PETSC_TRUE;
    char mesh_file_name[255];
    CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
                                 mesh_file_name, 255, &flg);
    if (flg != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "*** ERROR -my_file (MESH FILE NEEDED)");
    }

    // create MOAB communicator
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    if (pcomm == NULL)
      pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);

    const char *option;
    option = "";
    CHKERR moab.load_file(mesh_file_name, 0, option);

    // Create mofem interface
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;

    // Add meshsets with material and boundary conditions
    MeshsetsManager *meshsets_manager_ptr;
    CHKERR m_field.getInterface(meshsets_manager_ptr);
    CHKERR meshsets_manager_ptr->setMeshsetFromFile();

    PetscPrintf(PETSC_COMM_WORLD,
                "Read meshsets add added meshsets for bc.cfg\n");
    for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
      PetscPrintf(PETSC_COMM_WORLD, "%s",
                  static_cast<std::ostringstream &>(
                      std::ostringstream().seekp(0) << *it << endl)
                      .str()
                      .c_str());
      cerr << *it << endl;
    }

    // set entities bit level
    BitRefLevel ref_level;
    ref_level.set(0);
    CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
        0, 3, ref_level);

    // set app. order
    // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
    // (Mark Ainsworth & Joe Coyle)
    PetscInt order;
    CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &order,
                              &flg);
    if (flg != PETSC_TRUE) {
      order = 0;
    }

    // finite elements

    BcFluxMap bc_flux_map;
    MyTransport ufe(m_field, bc_flux_map);

    // Initially calculate problem on coarse mesh

    CHKERR ufe.addFields("VALUES", "FLUXES", order);
    CHKERR ufe.addFiniteElements("FLUXES", "VALUES");
    // Set boundary conditions
    CHKERR ufe.addBoundaryElements(ref_level);
    CHKERR ufe.buildProblem(ref_level);
    CHKERR ufe.createMatrices();
    CHKERR ufe.solveLinearProblem();
    CHKERR ufe.calculateResidual();
    CHKERR ufe.evaluateError();
    CHKERR ufe.destroyMatrices();
    CHKERR ufe.postProc("out_0.h5m");

    int nb_levels = 5; // default number of refinement levels
    // get number of refinement levels form command line
    CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-nb_levels", &nb_levels,
                              PETSC_NULL);

    // refine mesh, solve problem and do it again until number of refinement
    // levels are exceeded.
    for (int ll = 1; ll != nb_levels; ll++) {
      const int nb_levels = ll;
      CHKERR ufe.squashBits();
      CHKERR ufe.refineMesh(ufe, nb_levels, order);
      ref_level = BitRefLevel().set(nb_levels);
      bc_flux_map.clear();
      CHKERR ufe.addBoundaryElements(ref_level);
      CHKERR ufe.buildProblem(ref_level);
      CHKERR ufe.createMatrices();
      CHKERR ufe.solveLinearProblem();
      CHKERR ufe.calculateResidual();
      CHKERR ufe.evaluateError();
      CHKERR ufe.destroyMatrices();
      CHKERR ufe.postProc(
          static_cast<std::ostringstream &>(std::ostringstream().seekp(0)
                                            << "out_" << nb_levels << ".h5m")
              .str());
    }
  }
  CATCH_ERRORS;

  MoFEM::Core::Finalize();

  return 0;
}

Variable Documentation

help

char help[]

static

Initial value:

= "-my_file input file"
                     "-my_order order of approximation"
                     "-nb_levels number of refinements levels"
                     "\n\n"

Definition at line 32 of file h_adaptive_transport.cpp.