v0.13.1
hcurl_curl_operator.cpp

Using PipelineManager interface calculate the curl of base functions, and integral of the vector tangent vector with normal on the boundary. Since the h-curl space is used, volume integral and boundary integral should give the same result, as a result, as we are applying Stokes theorem on h-curl space.

/**
* \file hcurl_curl_operator.cpp
* \brief Testich curl-curl operator by applying Stokes theorem
* \example hcurl_curl_operator.cpp
*
* Using PipelineManager interface calculate the curl of base functions, and
* integral of the vector tangent vector with normal on the boundary. Since the
* h-curl space is used, volume integral and boundary integral should give the
* same result, as a result, as we are applying Stokes theorem on h-curl space.
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
"HCURL", UserDataOperator::OPROW),
tCurl(t_curl) {}
};
"HCURL", UserDataOperator::OPROW),
tCurl(t_curl) {}
};
int main(int argc, char *argv[]) {
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_IMPOSIBLE_CASE, "base not set");
}
PetscBool ho_geometry = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-ho_geometry", &ho_geometry,
PETSC_NULL);
DMType dm_name = "DMMOFEM";
// 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,
CHKERR simple_interface->addBoundaryField("HCURL", HCURL,
break;
case DEMKOWICZ:
CHKERR simple_interface->addDomainField("HCURL", HCURL,
CHKERR simple_interface->addBoundaryField("HCURL", HCURL,
break;
}
if (ho_geometry == PETSC_TRUE)
CHKERR simple_interface->addDataField("MESH_NODE_POSITIONS", H1,
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);
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(
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(
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpFacesRot(t_curl_skin));
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");
}
}
if (data.getFieldData().size() == 0)
const unsigned int nb_gauss_pts = data.getDiffN().size1();
const unsigned int nb_dofs = data.getFieldData().size();
MatrixDouble curl_mat;
auto t_curl_base = data.getFTensor2DiffN<3, 3>();
unsigned int gg = 0;
for (; gg < nb_gauss_pts; gg++) {
double w = getGaussPts()(3, gg) * getVolume();
for (unsigned int dd = 0; dd != nb_dofs; dd++) {
t_curl(i) = levi_civita(j, i, k) * t_curl_base(j, k);
tCurl(i) += w * t_curl(i);
++t_curl_base;
}
}
}
int nb_dofs = data.getFieldData().size();
if (nb_dofs == 0)
int nb_gauss_pts = data.getN().size1();
auto t_curl_base = data.getFTensor1N<3>();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
for (int dd = 0; dd < nb_dofs; dd++) {
double w = getGaussPts()(2, gg);
const double n0 = getNormalsAtGaussPts(gg)[0];
const double n1 = getNormalsAtGaussPts(gg)[1];
const double n2 = getNormalsAtGaussPts(gg)[2];
if (getFEType() == MBTRI) {
w *= 0.5;
}
tCurl(0) += (n1 * t_curl_base(2) - n2 * t_curl_base(1)) * w;
tCurl(1) += (n2 * t_curl_base(0) - n0 * t_curl_base(2)) * w;
tCurl(2) += (n0 * t_curl_base(1) - n1 * t_curl_base(0)) * w;
++t_curl_base;
}
}
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static const double eps
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
auto integration_rule
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
virtual bool check_field(const std::string &name) const =0
check if field is in database
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
int main(int argc, char *argv[])
static char help[]
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1086
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
double w(const double g, const double t)
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Calculate normals at Gauss points of triangle element.
transform Hcurl base fluxes from reference element to physical triangle
Apply covariant (Piola) transfer to Hcurl space for HO geometry.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Set inverse jacobian to base functions.
PipelineManager interface.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEMErrorCode addDataField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:430
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:374
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
MoFEMErrorCode addBoundaryField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on boundary.
Definition: Simple.cpp:392
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
OpFacesRot(FTensor::Tensor1< double, 3 > &t_curl)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
FTensor::Tensor1< double, 3 > & tCurl
OpVolCurl(FTensor::Tensor1< double, 3 > &t_curl)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
FTensor::Tensor1< double, 3 > & tCurl