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

Unit test for:

  1. Integration linear forms, and consistency of boundary integrals
  2. Integration integration on high-order geometry (2d and 3d)
  3. Integration for axi-symmetric case
  4. Consistency of Lhs operators, and Rhs operators (using OperatorsTester)

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 \]

/**
* \file tensor_divergence_operator.cpp
* \example tensor_divergence_operator.cpp
*
* Unit test for:
* 1. Integration linear forms, and consistency of boundary integrals
* 2. Integration integration on high-order geometry (2d and 3d)
* 3. Integration for axi-symmetric case
* 4. Consistency of Lhs operators, and Rhs operators (using OperatorsTester)
*
* Note: Two sets are tested, when Sigma or u is variation.
*
* \f[
* \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
* \f]
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType I =
IntegrationType::GAUSS; //< selected integration type
constexpr CoordinateTypes COORD_TYPE = EXECUTABLE_COORD_TYPE;
template <int DIM> struct ElementsAndOps {};
static constexpr FieldSpace HDIV_SPACE = HCURL;
};
static constexpr FieldSpace HDIV_SPACE = HDIV;
};
int main(int argc, char *argv[]) {
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
// Create MoFEM instance
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
// Declare elements
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
const char *list_bases[] = {"ainsworth", "demkowicz"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
LASBASETOP, &choice_base_value, &flg);
if (flg != PETSC_TRUE)
SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
if (choice_base_value == AINSWORTH)
else if (choice_base_value == DEMKOWICZ)
int order = 4;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
// Register DM Manager
DMType dm_name = "DMMOFEM";
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "ATOM"));
LogManager::setLog("ATOM");
MOFEM_LOG_TAG("ATOM", "ATOM");
// Simple interface
auto simple = m_field.getInterface<Simple>();
// get options from command line
CHKERR simple->getOptions();
// load mesh file
CHKERR simple->loadFile();
CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
CHKERR simple->addDomainField("SIGMA", HDIV_SPACE, DEMKOWICZ_JACOBI_BASE,
CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
CHKERR simple->addBoundaryField("SIGMA", HDIV_SPACE, DEMKOWICZ_JACOBI_BASE,
CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
// set fields order, i.e. for most first cases order is sufficient.
CHKERR simple->setFieldOrder("U", order);
// CHKERR simple->setFieldOrder("SIGMA", order);
CHKERR simple->setFieldOrder("GEOMETRY", 2);
auto get_skin = [&]() {
Range body_ents;
CHKERR m_field.get_moab().get_entities_by_dimension(0, SPACE_DIM,
body_ents);
Skinner skin(&m_field.get_moab());
Range skin_ents;
CHKERR skin.find_skin(0, body_ents, false, skin_ents);
return skin_ents;
};
auto filter_true_skin = [&](auto skin) {
Range boundary_ents;
ParallelComm *pcomm =
ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
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 simple->setFieldOrder("SIGMA", 0);
CHKERR simple->setFieldOrder("SIGMA", order, &boundary_ents);
// setup problem
CHKERR simple->setUp();
auto bc_mng = m_field.getInterface<BcManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM",
"U", 0, MAX_DOFS_ON_ENTITY, true);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(), "SYM", "SIGMA", 0, MAX_DOFS_ON_ENTITY, true);
auto project_ho_geometry = [&]() {
Projection10NodeCoordsOnField ent_method(m_field, "GEOMETRY");
return m_field.loop_dofs("GEOMETRY", ent_method);
};
CHKERR project_ho_geometry();
// get operators tester
auto opt = m_field.getInterface<OperatorsTester>(); // get interface to
// OperatorsTester
auto pip_mng = m_field.getInterface<PipelineManager>(); // get interface to
// pipeline manager
pip_mng->getOpDomainLhsPipeline().clear();
pip_mng->getOpDomainRhsPipeline().clear();
// integration rule
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);
auto x = opt->setRandomFields(simple->getDM(),
{{"U", {-1, 1}}, {"SIGMA", {-1, 1}}});
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), x, INSERT_VALUES,
SCATTER_REVERSE);
// set integration rules
auto jacobian = [&](double r) {
if constexpr (COORD_TYPE == CYLINDRICAL)
return 2 * M_PI * r;
else
return 1.;
};
auto beta_domain = [jacobian](double r, double, double) {
return jacobian(r);
};
auto beta_bdy = [jacobian](double r, double, double) {
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(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{}, {{"U", u_ptr}}, {{"SIGMA", sigma_ptr}}, {})
);
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), f_res, INSERT_VALUES,
SCATTER_REVERSE);
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
post_proc_fe);
post_proc_fe->writeFile(out_name);
};
auto test_consistency_of_domain_and_bdy_integrals = [&]() {
using OpMixDivURhs =
GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM, COORD_TYPE>;
using OpMixLambdaGradURhs = FormsIntegrators<DomainEleOp>::Assembly<
using OpMixLambdaGradURhs = FormsIntegrators<DomainEleOp>::Assembly<
using OpMixNormalLambdaURhs = FormsIntegrators<BoundaryEleOp>::Assembly<
using OpUTimeTractionRhs = FormsIntegrators<BoundaryEleOp>::Assembly<
auto ops_rhs_interior = [&](auto &pip) {
"GEOMETRY");
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
"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(
new OpMixUTimesDivLambdaRhs("U", sigma_div_ptr, beta_domain));
pip.push_back(new OpMixUTimesLambdaRhs("U", sigma_ptr, beta_domain));
};
auto ops_rhs_boundary = [&](auto &pip) {
"GEOMETRY");
auto u_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
auto traction_ptr = boost::make_shared<MatrixDouble>();
"SIGMA", traction_ptr));
// We have to integrate on curved face geometry, thus integration weight
// have to adjusted.
pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
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());
auto f = createDMVector(simple->getDM());
pip_mng->getDomainRhsFE()->f = f;
pip_mng->getBoundaryRhsFE()->f = f;
CHKERR VecZeroEntries(f);
simple->getDomainFEName(),
pip_mng->getDomainRhsFE());
simple->getBoundaryFEName(),
pip_mng->getBoundaryRhsFE());
CHKERR VecAssemblyBegin(f);
CHKERR VecAssemblyEnd(f);
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) {
CHKERR post_proc(simple->getDM(), f,
"tensor_divergence_operator_res_vec.h5m");
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Test failed");
}
};
// Testing Lhs domain operators
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(
new OpLambdaGraULhs("SIGMA", "U", unity, beta_domain, true, false));
};
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(),
pip_mng->getDomainLhsFE(), x, SmartPetscObj<Vec>(),
SmartPetscObj<Vec>(), diff_x, 0, 1, eps);
// Calculate norm of difference between directive calculated from finite
// difference, and tangent matrix.
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)
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Test failed");
};
CHKERR test_consistency_of_domain_and_bdy_integrals();
CHKERR test_lhs_ops();
}
}
static Index< 'p', 3 > p
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
static char help[]
int main()
Definition: adol-c_atom.cpp:46
static const double eps
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
FieldSpace
approximation spaces
Definition: definitions.h:82
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
CoordinateTypes
Coodinate system.
Definition: definitions.h:114
@ CYLINDRICAL
Definition: definitions.h:117
#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
constexpr int order
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
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
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
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
Definition: Common.hpp:10
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)
int r
Definition: sdf.py:8
constexpr IntegrationType I
constexpr AssemblyType A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
virtual moab::Interface & get_moab()=0
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)
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.
Definition: Simple.hpp:27
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