hcurl_curl_operator.cpp File Reference

Testich curl-curl operator by applying Stokes theorem. More...

#include <MoFEM.hpp>

Go to the source code of this file.

Classes
struct	OpVolCurl
struct	OpFacesRot

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

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

Detailed Description

Testich curl-curl operator by applying Stokes theorem.

Definition in file hcurl_curl_operator.cpp.

Function Documentation

main()

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

< database

< interface

Definition at line 45 of file hcurl_curl_operator.cpp.

                                 {
 
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    enum bases { AINSWORTH, DEMKOWICZ, LASTOP };
 
    const char *list[] = {"ainsworth", "demkowicz"};
 
    PetscBool flg;
    PetscInt choise_value = AINSWORTH;
    CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list, LASTOP,
                                &choise_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);
 
    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 (choise_value) {
    case AINSWORTH:
      CHKERR simple_interface->addDomainField("HCURL", HCURL,
                                              AINSWORTH_LEGENDRE_BASE, 1);
      CHKERR simple_interface->addBoundaryField("HCURL", HCURL,
                                                AINSWORTH_LEGENDRE_BASE, 1);
      break;
    case DEMKOWICZ:
      CHKERR simple_interface->addDomainField("HCURL", HCURL,
                                              DEMKOWICZ_JACOBI_BASE, 1);
      CHKERR simple_interface->addBoundaryField("HCURL", HCURL,
                                                DEMKOWICZ_JACOBI_BASE, 1);
      break;
    }
 
    if (ho_geometry == PETSC_TRUE)
      CHKERR simple_interface->addDataField("MESH_NODE_POSITIONS", H1,
                                            AINSWORTH_LEGENDRE_BASE, 3);
 
    constexpr int order = 3;
    CHKERR simple_interface->setFieldOrder("HCURL", order);
    // Range ents;
    // CHKERR moab.get_entities_by_dimension(0, 2, ents, true);
    // CHKERR simple_interface->setFieldOrder("HCURL", 1, &ents);
 
    if (ho_geometry == PETSC_TRUE)
      CHKERR simple_interface->setFieldOrder("MESH_NODE_POSITIONS", 2);
    CHKERR simple_interface->setUp();
 
    auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
    CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
    CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
 
    FTensor::Tensor1<double, 3> t_curl_vol;
    FTensor::Tensor1<double, 3> t_curl_skin;
 
    auto material_grad_mat = boost::make_shared<MatrixDouble>();
    auto material_det_vec = boost::make_shared<VectorDouble>();
    auto material_inv_grad_mat = boost::make_shared<MatrixDouble>();
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    auto det_ptr = boost::make_shared<VectorDouble>();
 
    boost::dynamic_pointer_cast<VolumeElementForcesAndSourcesCore>(
        pipeline_mng->getDomainRhsFE())
        ->meshPositionsFieldName = "none";
    boost::dynamic_pointer_cast<PipelineManager::FaceEle>(
        pipeline_mng->getBoundaryRhsFE())
        ->meshPositionsFieldName = "none";
 
    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 OpSetHOWeights(det_ptr));
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpSetHOCovariantPiolaTransform(HCURL, inv_jac_ptr));
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpSetHOInvJacVectorBase(HCURL, inv_jac_ptr));
 
    if (ho_geometry) {
      pipeline_mng->getOpDomainRhsPipeline().push_back(
          new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
                                                   material_grad_mat));
      pipeline_mng->getOpDomainRhsPipeline().push_back(new OpInvertMatrix<3>(
          material_grad_mat, material_det_vec, material_inv_grad_mat));
      pipeline_mng->getOpDomainRhsPipeline().push_back(
          new OpSetHOWeights(material_det_vec));
      pipeline_mng->getOpDomainRhsPipeline().push_back(
          new OpSetHOCovariantPiolaTransform(HCURL, material_inv_grad_mat));
      pipeline_mng->getOpDomainRhsPipeline().push_back(
          new OpSetHOInvJacVectorBase(HCURL, material_inv_grad_mat));
    } 
    pipeline_mng->getOpDomainRhsPipeline().push_back(new OpVolCurl(t_curl_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 OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
    pipeline_mng->getOpBoundaryRhsPipeline().push_back(
        new OpFacesRot(t_curl_skin));
 
    FTensor::Index<'i', 3> i;
    t_curl_vol(i) = 0;
    t_curl_skin(i) = 0;
 
    // 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);
    }
 
    // Run pipelines on mesh
    CHKERR pipeline_mng->loopFiniteElements();
 
    std::cout.precision(12);
 
    t_curl_vol(i) -= t_curl_skin(i);
    double nrm2 = sqrt(t_curl_vol(i) * t_curl_vol(i));
 
    constexpr double eps = 1e-8;
    if (fabs(nrm2) > eps)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Curl operator not passed test\n");
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
}

Variable Documentation

help

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

static

Definition at line 19 of file hcurl_curl_operator.cpp.