Integration on skeleton for 2dTeting integration on skeleton and checking of continuity of hcurl space on edges.
static char help[] =
"...\n\n";
FaceElementForcesAndSourcesCore::NO_CONTRAVARIANT_TRANSFORM_HDIV |
FaceElementForcesAndSourcesCore::NO_COVARIANT_TRANSFORM_HCURL>;
};
if (
type == MBEDGE && side == getEdgeSideNumber()) {
MatrixDouble diff = getCoordsAtGaussPts() - getEdgeCoordsAtGaussPts();
const double eps = 1e-12;
if (norm_inf(diff) >
eps)
"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();
auto t_hdiv_base = data.getFTensor1N<3>();
if (getFEMethod()->nInTheLoop == 0)
else
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;
}
}
}
}
};
faceSideFe(m_field), elemData(elem_data) {
}
const size_t nb_dofs = data.getN().size2() / 3;
const size_t nb_integration_pts = data.getN().size1();
auto &dir = getDirection();
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;
}
}
CHKERR loopSideFaces(
"dFE", faceSideFe);
auto check_continuity_of_base = [&](auto &vol_dot_data) {
if (vol_dot_data.size1() != elemData.
dotEdge.size1())
"Inconsistent number of integration points");
if (vol_dot_data.size2() != elemData.
dotEdge.size2())
"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));
"Inconsistency (%d, %d) %3.4e != %3.4e", gg, bb,
vol_dot_data(gg, bb), elemData.
dotEdge(gg, bb));
}
};
}
}
};
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
Simple *simple_interface;
{
CHKERR simple_interface->getOptions();
CHKERR simple_interface->loadFile(
"");
enum bases { AINSWORTH, DEMKOWICZ, LASTBASEOP };
const char *list_bases[] = {"ainsworth", "demkowicz"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
LASTBASEOP, &choice_base_value, &flg);
if (flg == PETSC_TRUE) {
if (choice_base_value == AINSWORTH)
else if (choice_base_value == DEMKOWICZ)
return base;
}
};
auto base = get_base();
CHKERR simple_interface->addDomainField(
"FIELD",
HCURL, base, 1);
CHKERR simple_interface->addSkeletonField(
"FIELD",
HCURL, base, 1);
CHKERR simple_interface->setFieldOrder(
"FIELD", 2);
CHKERR simple_interface->setUp();
auto dm = simple_interface->getDM();
boost::shared_ptr<EdgeEle> skeleton_fe =
boost::shared_ptr<EdgeEle>(
new EdgeEle(m_field));
skeleton_fe->getOpPtrVector().push_back(
new OpSetContravariantPiolaTransformOnEdge());
skeleton_fe->getOpPtrVector().push_back(
skeleton_fe);
}
}
return 0;
}
int main(int argc, char *argv[])
MoFEM::EdgeElementForcesAndSourcesCoreSwitch< EdgeElementForcesAndSourcesCore::NO_COVARIANT_TRANSFORM_HCURL > EdgeEle
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ HCURL
field with continuous tangents
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
FTensor::Index< 'i', 3 > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
implementation of Data Operators for Forces and Sources
OpCalculateJacForFaceImpl< 2 > OpCalculateJacForFace
OpSetContravariantPiolaTransformFaceImpl< 2 > OpSetContravariantPiolaTransformFace
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
DataForcesAndSourcesCore::EntData EntData
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
default operator for EDGE element
@ NO_COVARIANT_TRANSFORM_HCURL
default operator for Face element
FaceElementForcesAndSourcesCoreOnSideBase::UserDataOperator UserDataOperator
friend class UserDataOperator
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.