hdiv_divergence_operator.cpp File Reference

#include <MoFEM.hpp>

Go to the source code of this file.

Classes
struct	OpVolDivergence
struct	OpFacesFluxes

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

Variables
static char	help [] = "...\n\n"

Function Documentation

main()

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

< database

< interface

This has no real effect, folling line are only for atom test purpose

Definition at line 44 of file hdiv_divergence_operator.cpp.

                                 {
 
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    enum bases { AINSWORTH, DEMKOWICZ, LASTOP };
 
    const char *list[] = {"ainsworth", "demkowicz"};
 
    PetscBool flg;
    PetscInt choice_value = AINSWORTH;
    CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list, LASTOP,
                                &choice_value, &flg);
    if (flg != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
    }
 
    PetscBool ho_geometry = PETSC_FALSE;
    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-ho_geometry", &ho_geometry,
                               PETSC_NULL);
 
    PetscInt ho_choice_value = AINSWORTH;
    CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-ho_base", list, LASTOP,
                                &ho_choice_value, &flg);
 
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    // Create MoAB
    moab::Core mb_instance;              ///< database
    moab::Interface &moab = mb_instance; ///< interface
 
    // Create MoFEM
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    Simple *simple_interface = m_field.getInterface<Simple>();
    PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
    CHKERR simple_interface->getOptions();
    CHKERR simple_interface->loadFile("");
 
    // fields
    switch (choice_value) {
    case AINSWORTH:
      CHKERR simple_interface->addDomainField("HDIV", HDIV,
                                              AINSWORTH_LEGENDRE_BASE, 1);
      CHKERR simple_interface->addBoundaryField("HDIV", HDIV,
                                                AINSWORTH_LEGENDRE_BASE, 1);
      break;
    case DEMKOWICZ:
      CHKERR simple_interface->addDomainField("HDIV", HDIV,
                                              DEMKOWICZ_JACOBI_BASE, 1);
      CHKERR simple_interface->addBoundaryField("HDIV", HDIV,
                                                DEMKOWICZ_JACOBI_BASE, 1);
      break;
    }
 
    if (ho_geometry == PETSC_TRUE) {
      switch (ho_choice_value) {
      case AINSWORTH:
        CHKERR simple_interface->addDataField("MESH_NODE_POSITIONS", H1,
                                              AINSWORTH_LEGENDRE_BASE, 3);
      case DEMKOWICZ:
        CHKERR simple_interface->addDataField("MESH_NODE_POSITIONS", H1,
                                              DEMKOWICZ_JACOBI_BASE, 3);
      }
    }
 
    constexpr int order = 3;
    CHKERR simple_interface->setFieldOrder("HDIV", order);
    if (ho_geometry == PETSC_TRUE)
      CHKERR simple_interface->setFieldOrder("MESH_NODE_POSITIONS", 2);
    CHKERR simple_interface->setUp();
 
    /// This has no real effect, folling line are only for atom test purpose
    pipeline_mng->getDomainLhsFE().reset();
    pipeline_mng->getDomainRhsFE().reset();
    pipeline_mng->getBoundaryLhsFE().reset();
    pipeline_mng->getBoundaryRhsFE().reset();
    pipeline_mng->getSkeletonLhsFE().reset();
    pipeline_mng->getSkeletonRhsFE().reset();
 
    auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
    CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
    CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
 
    double divergence_vol = 0;
    double divergence_skin = 0;
 
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    auto det_ptr = boost::make_shared<VectorDouble>();
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpCalculateHOJac<3>(jac_ptr));
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpSetHOWeights(det_ptr));
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpVolDivergence(divergence_vol));
 
    if (m_field.check_field("MESH_NODE_POSITIONS"))
      pipeline_mng->getOpBoundaryRhsPipeline().push_back(
          new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
    pipeline_mng->getOpBoundaryRhsPipeline().push_back(
        new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
    pipeline_mng->getOpBoundaryRhsPipeline().push_back(
        new OpFacesFluxes(divergence_skin));
 
    // project geometry form 10 node tets on higher order approx. functions
    if (ho_geometry == PETSC_TRUE) {
      Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
      CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
    }
    CHKERR pipeline_mng->loopFiniteElements();
 
    MOFEM_LOG("WORLD", Sev::inform)
        << "divergence_vol " << std::setprecision(12) << divergence_vol;
    MOFEM_LOG("WORLD", Sev::inform)
        << "divergence_skin " << std::setprecision(12) << divergence_skin;
 
    constexpr double eps = 1e-8;
    const double error = divergence_skin - divergence_vol;
    if (fabs(error) > eps)
      SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
               "invalid surface flux or divergence or both error = %3.4e",
               error);
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
}

Variable Documentation

help

char help[] = "...\n\n"

static

Definition at line 16 of file hdiv_divergence_operator.cpp.