Note: Two sets are tested, when Sigma or u is variation.
\[
\int_\Gamma n_i \Sigma_{ij} u_j \textrm{d}\Gamma
=
\int_\Omega \left( \Sigma_{ij} u_j \right)_i \textrm{d}\Omega
=
\int_\Omega \Sigma_{ij,i} u_j \textrm{d}\Omega
+
\int_\Omega \Sigma_{ij} u_{j,i} \textrm{d}\Omega
\]
static char help[] =
"...\n\n";
IntegrationType::GAUSS;
};
};
int main(
int argc,
char *argv[]) {
try {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
const char *list_bases[] = {"ainsworth", "demkowicz"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
LASBASETOP, &choice_base_value, &flg);
if (flg != PETSC_TRUE)
if (choice_base_value == AINSWORTH)
else if (choice_base_value == DEMKOWICZ)
DMType dm_name = "DMMOFEM";
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "ATOM"));
LogManager::setLog("ATOM");
auto get_skin = [&]() {
body_ents);
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
auto filter_true_skin = [&](auto skin) {
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
auto boundary_ents = filter_true_skin(get_skin());
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM",
CHKERR bc_mng->removeBlockDOFsOnEntities(
auto project_ho_geometry = [&]() {
return m_field.
loop_dofs(
"GEOMETRY", ent_method);
};
pip_mng->getOpDomainLhsPipeline().clear();
pip_mng->getOpDomainRhsPipeline().clear();
auto rule = [&](
int,
int,
int p) {
return 2 *
p + 2; };
CHKERR pip_mng->setDomainRhsIntegrationRule(rule);
CHKERR pip_mng->setDomainLhsIntegrationRule(rule);
CHKERR pip_mng->setBoundaryRhsIntegrationRule(rule);
{{"U", {-1, 1}}, {"SIGMA", {-1, 1}}});
SCATTER_REVERSE);
auto jacobian = [&](
double r) {
return 2 * M_PI * r;
else
return 1.;
};
return jacobian(r);
};
return -jacobian(r);
};
auto post_proc = [&](auto dm, auto f_res, auto out_name) {
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
auto sigma_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
sigma_ptr));
auto u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{}, {{"U", u_ptr}}, {{"SIGMA", sigma_ptr}}, {})
);
SCATTER_REVERSE);
post_proc_fe);
post_proc_fe->writeFile(out_name);
};
auto test_consistency_of_domain_and_bdy_integrals = [&]() {
GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM, COORD_TYPE>;
auto ops_rhs_interior = [&](auto &pip) {
"GEOMETRY");
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
"U", grad_u_ptr));
pip.push_back(
new OpMixDivURhs(
"SIGMA", u_ptr, beta_domain));
pip.push_back(
new OpMixLambdaGradURhs("SIGMA", grad_u_ptr, beta_domain));
auto sigma_ptr = boost::make_shared<MatrixDouble>();
auto sigma_div_ptr = boost::make_shared<MatrixDouble>();
"SIGMA", sigma_div_ptr));
"SIGMA", sigma_ptr));
pip.push_back(
};
auto ops_rhs_boundary = [&](auto &pip) {
"GEOMETRY");
auto u_ptr = boost::make_shared<MatrixDouble>();
auto traction_ptr = boost::make_shared<MatrixDouble>();
"SIGMA", traction_ptr));
pip.push_back(new OpMixNormalLambdaURhs("SIGMA", u_ptr, beta_bdy));
pip.push_back(new OpUTimeTractionRhs("U", traction_ptr, beta_bdy));
};
CHKERR ops_rhs_interior(pip_mng->getOpDomainRhsPipeline());
CHKERR ops_rhs_boundary(pip_mng->getOpBoundaryRhsPipeline());
pip_mng->getDomainRhsFE()->f =
f;
pip_mng->getBoundaryRhsFE()->f =
f;
pip_mng->getDomainRhsFE());
pip_mng->getBoundaryRhsFE());
double f_nrm2;
CHKERR VecNorm(f, NORM_2, &f_nrm2);
MOFEM_LOG(
"ATOM", Sev::inform) <<
"f_norm2 = " << f_nrm2;
if (std::fabs(f_nrm2) > 1e-10) {
"tensor_divergence_operator_res_vec.h5m");
}
};
auto test_lhs_ops = [&]() {
auto op_lhs_domain = [&](auto &pip) {
"GEOMETRY");
auto unity = []() { return 1; };
pip.push_back(
new OpMixDivULhs(
"SIGMA",
"U", unity, beta_domain,
true,
false));
pip.push_back(
};
CHKERR op_lhs_domain(pip_mng->getOpDomainLhsPipeline());
auto diff_x = opt->setRandomFields(
simple->getDM(),
{{"U", {-1, 1}}, {"SIGMA", {-1, 1}}});
constexpr double eps = 1e-5;
auto diff_res = opt->checkCentralFiniteDifference(
simple->getDM(),
simple->getDomainFEName(), pip_mng->getDomainRhsFE(),
double fnorm_res;
CHKERR VecNorm(diff_res, NORM_2, &fnorm_res);
MOFEM_LOG_C(
"ATOM", Sev::inform,
"Test Lhs OPs %3.4e", fnorm_res);
if (std::fabs(fnorm_res) > 1e-8)
};
CHKERR test_consistency_of_domain_and_bdy_integrals();
}
}
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
FieldSpace
approximation spaces
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
CoordinateTypes
Coodinate system.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
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.
auto createDMVector(DM dm)
Get smart vector from DM.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
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.
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
constexpr IntegrationType I
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
virtual moab::Interface & get_moab()=0
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate trace of vector (Hdiv/Hcurl) space.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
SmartPetscObj< Vec > setRandomFields(SmartPetscObj< DM > dm, std::vector< RandomFieldData > random_fields, boost::shared_ptr< Range > ents=nullptr)
Generate random fileds.
PipelineManager interface.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
constexpr FieldSpace HDIV_SPACE
constexpr CoordinateTypes COORD_TYPE