v0.14.0
Loading...
Searching...
No Matches
continuity_check_on_skeleton_with_simple_2d_for_hdiv.cpp

Integration on skeleton for 2d.

Integration on skeleton for 2d

Teting integration on skeleton and checking of continuity of hcurl space on edges.

/**
* \file continuity_check_on_skeleton_with_simple_2d_for_hdiv.cpp
* \ingroup mofem_simple_interface
* \example continuity_check_on_skeleton_with_simple_2d_for_hdiv.cpp
*
* \brief Integration on skeleton for 2d
*
* Teting integration on skeleton and checking of continuity of hcurl space on
* edges.
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
struct CommonData {
};
struct SkeletonFE : public EdgeEleOp {
struct OpFaceSide : public FaceEleOnSideOp {
OpFaceSide(CommonData &elem_data)
: FaceEleOnSideOp("FIELD", UserDataOperator::OPROW), elemData(elem_data) {}
if (type == MBEDGE && side == getEdgeSideNumber()) {
const double eps = 1e-12;
if (norm_inf(diff) > eps)
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
"coordinates at integration pts are different");
const size_t nb_dofs = data.getN().size2() / 3;
const size_t nb_integration_pts = data.getN().size1();
auto &dir = getDirection();
FTensor::Tensor1<double, 2> t_normal{-dir(1), dir(0)};
auto t_hdiv_base = data.getFTensor1N<3>();
MatrixDouble *ptr_dot_elem_data = nullptr;
if (getFEMethod()->nInTheLoop == 0)
ptr_dot_elem_data = &elemData.dotEleLeft;
else
ptr_dot_elem_data = &elemData.dotEleRight;
MatrixDouble &dot_elem_data = *ptr_dot_elem_data;
dot_elem_data.resize(nb_integration_pts, nb_dofs, false);
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
for (size_t bb = 0; bb != nb_dofs; ++bb) {
dot_elem_data(gg, bb) = t_normal(i) * t_hdiv_base(i);
++t_hdiv_base;
}
}
}
}
};
SkeletonFE(MoFEM::Interface &m_field, CommonData &elem_data)
faceSideFe(m_field), elemData(elem_data) {
auto jac_ptr = boost::make_shared<MatrixDouble>();
new OpCalculateHOJacForFace(jac_ptr));
}
if (type == MBEDGE) {
const size_t nb_dofs = data.getN().size2() / 3;
const size_t nb_integration_pts = data.getN().size1();
auto &dir = getDirection();
FTensor::Tensor1<double, 2> t_normal{-dir(1), dir(0)};
auto t_hdiv_base = data.getFTensor1N<3>();
elemData.dotEdge.resize(nb_integration_pts, nb_dofs, false);
elemData.dotEleLeft.resize(nb_integration_pts, 0, false);
elemData.dotEleRight.resize(nb_integration_pts, 0, false);
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
for (size_t bb = 0; bb != nb_dofs; ++bb) {
elemData.dotEdge(gg, bb) = t_normal(i) * t_hdiv_base(i);
++t_hdiv_base;
}
}
auto check_continuity_of_base = [&](auto &vol_dot_data) {
if (vol_dot_data.size1() != elemData.dotEdge.size1())
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Inconsistent number of integration points");
if (vol_dot_data.size2() != elemData.dotEdge.size2())
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Inconsistent number of base functions");
const double eps = 1e-12;
for (size_t gg = 0; gg != vol_dot_data.size1(); ++gg)
for (size_t bb = 0; bb != vol_dot_data.size2(); ++bb) {
const double error =
std::abs(vol_dot_data(gg, bb) - elemData.dotEdge(gg, bb));
if (error > eps)
SETERRQ4(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Inconsistency (%d, %d) %3.4e != %3.4e", gg, bb,
vol_dot_data(gg, bb), elemData.dotEdge(gg, bb));
}
};
if (elemData.dotEleLeft.size2() != 0)
CHKERR check_continuity_of_base(elemData.dotEleLeft);
else if (elemData.dotEleRight.size2() != 0)
CHKERR check_continuity_of_base(elemData.dotEleRight);
}
}
};
int main(int argc, char *argv[]) {
// initialize petsc
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
// Create MoAB database
moab::Core moab_core;
moab::Interface &moab = moab_core;
// Create MoFEM database and link it to MoAB
MoFEM::Core mofem_core(moab);
MoFEM::Interface &m_field = mofem_core;
// Register DM Manager
DMType dm_name = "DMMOFEM";
// Simple interface
Simple *simple_interface;
CHKERR m_field.getInterface(simple_interface);
{
// get options from command line
CHKERR simple_interface->getOptions();
// load mesh file
CHKERR simple_interface->loadFile("");
auto get_base = []() -> FieldApproximationBase {
enum bases { AINSWORTH, DEMKOWICZ, LASTBASEOP };
const char *list_bases[] = {"ainsworth", "demkowicz"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
LASTBASEOP, &choice_base_value, &flg);
if (flg == PETSC_TRUE) {
if (choice_base_value == AINSWORTH)
else if (choice_base_value == DEMKOWICZ)
return base;
}
return LASTBASE;
};
// add fields
auto base = get_base();
CHKERR simple_interface->addDomainField("FIELD", HCURL, base, 1);
CHKERR simple_interface->addSkeletonField("FIELD", HCURL, base, 1);
// set fields order
CHKERR simple_interface->setFieldOrder("FIELD", 2);
// setup problem
CHKERR simple_interface->setUp();
// get dm
auto dm = simple_interface->getDM();
// create elements
CommonData elem_data;
boost::shared_ptr<EdgeEle> skeleton_fe =
boost::shared_ptr<EdgeEle>(new EdgeEle(m_field));
skeleton_fe->getOpPtrVector().push_back(
skeleton_fe->getOpPtrVector().push_back(
new SkeletonFE(m_field, elem_data));
// iterate skeleton finite elements
skeleton_fe);
}
}
// finish work cleaning memory, getting statistics, etc.
return 0;
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static char help[]
int main()
Definition: adol-c_atom.cpp:46
static const double eps
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ 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
@ 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_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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
FTensor::Index< 'i', SPACE_DIM > i
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: Common.hpp:10
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
int nInTheLoop
number currently of processed method
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.
MoFEMErrorCode loopSideFaces(const string fe_name, FaceElementForcesAndSourcesCoreOnSide &fe_side)
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Base face element used to integrate on skeleton.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
@ OPROW
operator doWork function is executed on FE rows
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Make Hdiv space from Hcurl space in 2d.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
Definition: Simple.hpp:341
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:264
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEMErrorCode addSkeletonField(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 skeleton.
Definition: Simple.cpp:300
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.