Integration on skeleton for 2d.
Testing integration on skeleton and checking of continuity of hcurl space on edges.
static char help[] =
"...\n\n";
};
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();
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_tangent(
i) * t_hcurl_base(
i);
++t_hcurl_base;
}
}
}
}
};
}
if (type == MBEDGE) {
const size_t nb_dofs = data.
getN().size2() / 3;
const size_t nb_integration_pts = data.
getN().size1();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
for (size_t bb = 0; bb != nb_dofs; ++bb) {
++t_hcurl_base;
}
}
auto check_continuity_of_base = [&](auto &vol_dot_data) {
"Inconsistent number of integration points");
"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 =
"Inconsistency (%d, %d) %3.4e != %3.4e", gg, bb,
else
}
};
}
}
};
int main(
int argc,
char *argv[]) {
try {
moab::Core moab_core;
moab::Interface &moab = moab_core;
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmSelf(), "ATOM"));
LogManager::setLog("ATOM");
DMType dm_name = "DMMOFEM";
{
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();
auto dm = simple_interface->
getDM();
boost::shared_ptr<EdgeEle> skeleton_fe =
boost::shared_ptr<EdgeEle>(
new EdgeEle(m_field));
skeleton_fe->getOpPtrVector().push_back(
skeleton_fe);
}
}
return 0;
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MoFEM::EdgeElementForcesAndSourcesCore 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 DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Add operators pushing bases from local to physical configuration.
int nInTheLoop
number currently of processed method
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
FTensor::Tensor1< double, 3 > getFTensor1Direction()
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.
default operator for Face element
int getEdgeSideNumber() const
get face side number in respect to volume
MatrixDouble & getEdgeCoordsAtGaussPts()
get face coordinates at Gauss pts.
auto getFTensor1Direction()
get normal as tensor
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.
Simple interface for fast problem set-up.
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
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.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
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.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
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.