v0.13.1
Loading...
Searching...
No Matches
hanging_node_approx.cpp

Tetsing approximation with hanging nodes.

/**
* \example hanging_node_approx.cpp
*
* Tetsing approximation with hanging nodes.
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
constexpr char FIELD_NAME[] = "U";
constexpr int FIELD_DIM = 1;
constexpr int SPACE_DIM = 2;
constexpr int nb_ref_levels = 3; ///< Three levels of refinement
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
template <> struct ElementsAndOps<3> {
};
template <int FIELD_DIM> struct ApproxFieldFunction;
template <int FIELD_DIM> struct ApproxFieldFunctionDerivative;
/**
* @brief third order polynomial used for testing
*
*/
template <> struct ApproxFieldFunction<1> {
auto operator()(const double x, const double y, const double z) {
return x * x + y * y + x * y * y + x * x * y;
}
};
/**
* @brief third order polynomial used for testing
*
*/
template <> struct ApproxFieldFunctionDerivative<1> {
auto operator()(const double x, const double y, const double z) {
// x * x + y * y + x * y * y + x * x * y
return FTensor::Tensor1<double, SPACE_DIM>{2 * x + y * y + 2 * x * y,
2 * y + 2 * x * y + x * x};
}
};
/**
* @brief evaluate mass matrix
*
*/
/**
* @brief evaluate source, i.e. rhs vector
*
*/
/**
* @brief set bit
*
*/
auto bit = [](auto l) { return BitRefLevel().set(l); };
/**
* @brief set bit to marker
*
* Marker is used to mark field entities on skin on which we have hanging nodes
*/
auto marker = [](auto l) {
return BitRefLevel().set(BITREFLEVEL_SIZE - 1 - l);
};
/**
* @brief set levels of projection operators, which project field data from
* parent entities, to child, up to to level, i.e. last mesh refinement.
*
*/
template <typename PARENT_FE>
boost::shared_ptr<FEMethod> &fe_top,
ForcesAndSourcesCore::UserDataOperator::OpType op,
int verbosity, LogManager::SeverityLevel sev) {
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
BitRefLevel bit_marker;
for (auto l = 1; l <= nb_ref_levels; ++l)
bit_marker |= marker(l);
boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
add_parent_level =
[&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
if (level > 0) {
auto fe_ptr_current = boost::shared_ptr<ForcesAndSourcesCore>(
new PARENT_FE(m_field));
if (op == DomainEleOp::OPSPACE) {
if (typeid(PARENT_FE) == typeid(DomainParentEle)) {
fe_ptr_current->getOpPtrVector().push_back(
new OpCalculateHOJac<2>(jac_ptr));
fe_ptr_current->getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
fe_ptr_current->getOpPtrVector().push_back(
new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
}
}
add_parent_level(
boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
fe_ptr_current),
level - 1);
if (op == DomainEleOp::OPSPACE) {
parent_fe_pt->getOpPtrVector().push_back(
H1, op, fe_ptr_current,
BitRefLevel().set(), bit(0).flip(),
bit_marker, BitRefLevel().set(),
verbosity, sev));
} else {
parent_fe_pt->getOpPtrVector().push_back(
FIELD_NAME, op, fe_ptr_current,
BitRefLevel().set(), bit(0).flip(),
bit_marker, BitRefLevel().set(),
verbosity, sev));
}
}
};
add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
};
/**
* @brief lambda function used to select elements on which finite element
* pipelines are executed.
*
* @note childs elements on pipeline, retrive data from parents using operators
* pushed by \ref set_parent_dofs
*
*/
auto test_bit_child = [](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
struct AtomTest {
AtomTest(MoFEM::Interface &m_field) : mField(m_field) {}
private:
/**
* @brief red mesh and randomly refine three times
*
* @return MoFEMErrorCode
*/
/**
* @brief add field, and set up problem
*
* @return MoFEMErrorCode
*/
struct CommonData {
boost::shared_ptr<VectorDouble> approxVals;
boost::shared_ptr<MatrixDouble> divApproxVals;
};
template <int FIELD_DIM> struct OpError;
template <int FIELD_DIM> struct OpErrorSkel;
};
template <> struct AtomTest::OpError<1> : public DomainEleOp {
boost::shared_ptr<CommonData> commonDataPtr;
OpError(boost::shared_ptr<CommonData> &common_data_ptr)
: DomainEleOp(FIELD_NAME, OPROW), commonDataPtr(common_data_ptr) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
if (const size_t nb_dofs = data.getIndices().size()) {
const int nb_integration_pts = getGaussPts().size2();
auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
auto t_grad_val =
getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->divApproxVals));
auto t_coords = getFTensor1CoordsAtGaussPts();
VectorDouble nf(nb_dofs, false);
nf.clear();
const double volume = getMeasure();
auto t_row_base = data.getFTensor0N();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * volume;
double diff = t_val - AtomTest::approxFunction(t_coords(0), t_coords(1),
t_coords(2));
auto t_grad_diff =
AtomTest::divApproxFunction(t_coords(0), t_coords(1), t_coords(2));
t_grad_diff(i) -= t_grad_val(i);
error += alpha * (pow(diff, 2) + t_grad_diff(i) * t_grad_diff(i));
for (size_t r = 0; r != nb_dofs; ++r) {
nf[r] += alpha * t_row_base * diff;
++t_row_base;
}
++t_w;
++t_val;
++t_grad_val;
++t_coords;
}
const int index = 0;
CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
}
}
};
template <> struct AtomTest::OpErrorSkel<1> : public BoundaryEleOp {
boost::shared_ptr<CommonData> commonDataPtr;
OpErrorSkel(boost::shared_ptr<CommonData> &common_data_ptr)
: BoundaryEleOp(H1, OPSPACE), commonDataPtr(common_data_ptr) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
auto t_coords = getFTensor1CoordsAtGaussPts();
const double volume = getMeasure();
double error2 = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * volume;
double diff = t_val - AtomTest::approxFunction(t_coords(0), t_coords(1),
t_coords(2));
error2 += alpha * (pow(diff, 2));
++t_w;
++t_val;
++t_coords;
}
MOFEM_LOG("SELF", Sev::verbose) << "Boundary error " << sqrt(error2);
constexpr double eps = 1e-8;
if (sqrt(error2) > eps)
SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Error on boundary = %6.4e", sqrt(error2));
}
};
//! [Run programme]
}
//! [Run programme]
//! [Read mesh]
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
Skinner skin(&mField.get_moab());
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
MOFEM_LOG("WORLD", Sev::verbose) << "Dim " << simpleInterface->getDim();
auto &moab = mField.get_moab();
Range level0_ents;
CHKERR mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
bit(0), BitRefLevel().set(), SPACE_DIM, level0_ents);
Range level0_skin;
CHKERR skin.find_skin(0, level0_ents, false, level0_skin);
CHKERR pcomm->filter_pstatus(level0_skin,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr);
auto refine_mesh = [&](auto l) {
auto meshset_level0_ptr = get_temp_meshset_ptr(moab);
SPACE_DIM, *meshset_level0_ptr);
// random mesh refinement
auto meshset_ref_edges_ptr = get_temp_meshset_ptr(moab);
Range els;
CHKERR moab.get_entities_by_dimension(*meshset_level0_ptr, SPACE_DIM, els);
CHKERR bit_mng->filterEntitiesByRefLevel(bit(l - 1), bit(l - 1), els);
Range ele_to_refine;
if (l == 1) {
int ii = 0;
for (auto t : els) {
if ((ii % 2)) {
ele_to_refine.insert(t);
std::vector<EntityHandle> adj_edges;
CHKERR mField.get_moab().get_adjacencies(&t, 1, SPACE_DIM - 1, false,
adj_edges);
CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*adj_edges.begin(),
adj_edges.size());
}
++ii;
}
} else {
Range level_skin;
CHKERR skin.find_skin(0, els, false, level_skin);
CHKERR pcomm->filter_pstatus(level_skin,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr);
level_skin = subtract(level_skin, level0_skin);
Range adj;
CHKERR mField.get_moab().get_adjacencies(level_skin, SPACE_DIM, false,
adj, moab::Interface::UNION);
els = subtract(els, adj);
ele_to_refine.merge(els);
Range adj_edges;
CHKERR mField.get_moab().get_adjacencies(
els, SPACE_DIM - 1, false, adj_edges, moab::Interface::UNION);
CHKERR moab.add_entities(*meshset_ref_edges_ptr, adj_edges);
}
CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr, bit(l),
false, VERBOSE);
CHKERR refine->refineTrisHangingNodes(*meshset_level0_ptr, bit(l), VERBOSE);
CHKERR bit_mng->updateRangeByChildren(level0_skin, level0_skin);
simpleInterface->getBoundaryMeshSet(), bit(l - 1), BitRefLevel(),
bit(l), BitRefLevel(), simpleInterface->getBoundaryMeshSet(), MBMAXTYPE,
true);
bit(l), BitRefLevel().set(), SPACE_DIM,
(boost::lexical_cast<std::string>(l) + "_ref_mesh.vtk").c_str(), "VTK",
"");
bit(l), bit(l), MBTRI,
(boost::lexical_cast<std::string>(l) + "_only_ref_mesh.vtk").c_str(),
"VTK", "");
};
auto mark_skins = [&](auto l, auto m) {
Range ents;
ents);
Range level_skin;
CHKERR skin.find_skin(0, ents, false, level_skin);
CHKERR pcomm->filter_pstatus(level_skin,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr);
level_skin = subtract(level_skin, level0_skin);
CHKERR mField.get_moab().get_adjacencies(level_skin, 0, false, level_skin,
moab::Interface::UNION);
CHKERR bit_mng->addBitRefLevel(level_skin, marker(m));
};
BitRefLevel bit_sum;
for (auto l = 0; l != nb_ref_levels; ++l) {
CHKERR refine_mesh(l + 1);
CHKERR mark_skins(l, l + 1);
CHKERR mark_skins(l + 1, l + 1);
bit_sum |= bit(l);
}
bit_sum |= bit(nb_ref_levels);
simpleInterface->getBitRefLevel() = bit_sum;
simpleInterface->getBitRefLevelMask() = BitRefLevel().set();
}
//! [Read mesh]
//! [Set up problem]
// Add field
CHKERR simpleInterface->addDomainField(FIELD_NAME, H1,
CHKERR simpleInterface->addBoundaryField(FIELD_NAME, H1,
constexpr int order = 3;
// Simple interface will resolve adjacency to DOFs of parent of the element.
// Using that information MAtrixManager allocate appropriately size of
// matrix.
simpleInterface->getParentAdjacencies() = true;
BitRefLevel bit_marker;
for (auto l = 1; l <= nb_ref_levels; ++l)
bit_marker |= marker(l);
simpleInterface->getBitAdjEnt() = bit_marker;
// remove obsolete DOFs from problem
for (int l = 0; l != nb_ref_levels; ++l) {
CHKERR prb_mng->removeDofsOnEntities(simpleInterface->getProblemName(),
CHKERR prb_mng->removeDofsOnEntities(simpleInterface->getProblemName(),
bit(l).flip());
}
}
//! [Set up problem]
//! [Push operators to pipeline]
auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
pipeline_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
auto beta = [](const double, const double, const double) { return 1; };
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
DomainEleOp::OPSPACE, QUIET, Sev::noisy);
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
DomainEleOp::OPROW, QUIET, Sev::noisy);
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
DomainEleOp::OPCOL, QUIET, Sev::noisy);
pipeline_mng->getOpDomainLhsPipeline().push_back(
auto field_op_row = new ForcesAndSourcesCore::UserDataOperator(
field_op_row->doWorkRhsHook = [](DataOperator *op_ptr, int side,
if (type == MBENTITYSET) {
MOFEM_LOG("SELF", Sev::verbose)
<< "ROW: side/type: " << side << "/" << CN::EntityTypeName(type)
<< " op space/base: " << FieldSpaceNames[data.getSpace()] << "/"
<< ApproximationBaseNames[data.getBase()] << " DOFs "
<< data.getIndices() << " nb base functions " << data.getN().size2()
<< " nb base functions integration points " << data.getN().size1();
auto get_bool = [](auto fe, auto bit) {
return (bit & fe->getBitRefLevel()).any();
};
for (auto &field_ent : data.getFieldEntities()) {
MOFEM_LOG("SELF", Sev::verbose)
<< "\t" << CN::EntityTypeName(field_ent->getEntType());
}
}
};
set_parent_dofs<DomainParentEle>(
Sev::verbose);
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
DomainEleOp::OPROW, VERBOSE, Sev::noisy);
pipeline_mng->getOpDomainRhsPipeline().push_back(field_op_row);
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
//! [Push operators to pipeline]
//! [Solve]
MOFEM_LOG("WORLD", Sev::inform) << "Solve problem";
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto dm = simpleInterface->getDM();
auto D = smartCreateDMVector(dm);
CHKERR KSPSolve(solver, F, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
//! [Check results]
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
pipeline_mng->getBoundaryRhsFE()->exeTestHook = test_bit_child;
auto common_data_ptr = boost::make_shared<CommonData>();
common_data_ptr->resVec = smartCreateDMVector(simpleInterface->getDM());
common_data_ptr->L2Vec = createSmartVectorMPI(
mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
common_data_ptr->divApproxVals = boost::make_shared<MatrixDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHOJac<2>(jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
DomainEleOp::OPSPACE, QUIET, Sev::noisy);
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
DomainEleOp::OPROW, VERBOSE, Sev::noisy);
pipeline_mng->getOpDomainRhsPipeline().push_back(
FIELD_NAME, common_data_ptr->divApproxVals));
pipeline_mng->getOpDomainRhsPipeline().push_back(
common_data_ptr->approxVals));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpError<FIELD_DIM>(common_data_ptr));
set_parent_dofs<BoundaryParentEle>(mField, pipeline_mng->getBoundaryRhsFE(),
set_parent_dofs<BoundaryParentEle>(mField, pipeline_mng->getBoundaryRhsFE(),
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
common_data_ptr->approxVals));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpErrorSkel<FIELD_DIM>(common_data_ptr));
CHKERR VecZeroEntries(common_data_ptr->L2Vec);
CHKERR VecZeroEntries(common_data_ptr->resVec);
CHKERR pipeline_mng->loopFiniteElements();
CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
CHKERR VecAssemblyBegin(common_data_ptr->resVec);
CHKERR VecAssemblyEnd(common_data_ptr->resVec);
double nrm2;
CHKERR VecNorm(common_data_ptr->resVec, NORM_2, &nrm2);
const double *array;
CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
MOFEM_LOG_C("WORLD", Sev::inform, "Error %6.4e Vec norm %6.4e\n",
std::sqrt(array[0]), nrm2);
constexpr double eps = 1e-8;
if (nrm2 > eps)
SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Not converged solution err = %6.4e", nrm2);
if (std::sqrt(array[0]) > eps)
SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Error in approximation err = %6.4e", std::sqrt(array[0]));
CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
}
//! [Check results]
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
auto rule = [](int, int, int p) -> int { return -1; };
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
static_cast<ForcesAndSourcesCore *>(pipeline_mng->getDomainRhsFE().get())
->setRuleHook = [](ForcesAndSourcesCore *fe_raw_ptr, int order_row,
int order_col, int order_data) -> MoFEMErrorCode {
fe_raw_ptr->gaussPts.resize(3, 3);
fe_raw_ptr->gaussPts(0, 0) = 0;
fe_raw_ptr->gaussPts(1, 0) = 0;
fe_raw_ptr->gaussPts(2, 0) = 0;
fe_raw_ptr->gaussPts(0, 1) = 1;
fe_raw_ptr->gaussPts(1, 1) = 0;
fe_raw_ptr->gaussPts(2, 1) = 0;
fe_raw_ptr->gaussPts(0, 2) = 0;
fe_raw_ptr->gaussPts(1, 2) = 1;
fe_raw_ptr->gaussPts(2, 2) = 0;
};
auto field_op_row = new ForcesAndSourcesCore::UserDataOperator(
auto approx_vals = boost::make_shared<VectorDouble>();
auto &moab = mField.get_moab();
Tag th;
double def_val[] = {0};
CHKERR moab.tag_get_handle("FIELD", 1, MB_TYPE_DOUBLE, th,
MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
field_op_row->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
if (type == MBVERTEX) {
auto op_ptr =
base_op_ptr);
auto t_field = getFTensor0FromVec(*approx_vals);
auto nb_gauss_pts = op_ptr->getGaussPts().size2();
if (nb_gauss_pts != 3)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Should be three guass pts.");
auto conn = op_ptr->getConn();
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
const double v = t_field;
CHKERR moab.tag_set_data(th, &conn[gg], 1, &v);
++t_field;
}
}
};
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
DomainEleOp::OPROW, VERBOSE, Sev::noisy);
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(field_op_row);
CHKERR pipeline_mng->loopFiniteElements();
CHKERR mField.getInterface<BitRefManager>()->writeBitLevelByType(
bit(nb_ref_levels), BitRefLevel().set(), MBTRI, "out.vtk", "VTK", "");
}
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
MoFEM::Core::Initialize(&argc, &argv, NULL, help);
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [AtomTest]
AtomTest ex(m_field);
CHKERR ex.runProblem();
//! [AtomTest]
}
}
static Index< 'p', 3 > p
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
static char help[]
int main()
Definition: adol-c_atom.cpp:46
static const double eps
ElementsAndOps< SPACE_DIM >::DomianParentEle DomainParentEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
constexpr char FIELD_NAME[]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
constexpr int SPACE_DIM
constexpr int FIELD_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
@ QUIET
Definition: definitions.h:208
@ VERBOSE
Definition: definitions.h:209
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:219
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
static const char *const FieldSpaceNames[]
Definition: definitions.h:92
#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
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
static const char *const ApproximationBaseNames[]
Definition: definitions.h:72
#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
FTensor::Index< 'm', SPACE_DIM > m
ElementsAndOps< SPACE_DIM >::DomianParentEle DomianParentEle
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
MoFEMErrorCode updateMeshsetByEntitiesChildren(const EntityHandle parent, const BitRefLevel &parent_bit, const BitRefLevel &parent_mask, const BitRefLevel &child_bit, const BitRefLevel &child_mask, const EntityHandle child, EntityType child_type, const bool recursive=false, int verb=0)
Get child entities form meshset containing parent entities.
MoFEMErrorCode getEntitiesByDimAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const int dim, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
MoFEMErrorCode addBitRefLevel(const Range &ents, const BitRefLevel bit, int verb=QUIET) const
add bit ref level to ref entity
MoFEMErrorCode updateRangeByChildren(const Range &parent, Range &child, MoFEMTypes bh=MF_ZERO)
Update range by childresn.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
auto marker
set bit to marker
constexpr int nb_ref_levels
Three levels of refinement.
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
auto bit
set bit
auto set_parent_dofs(MoFEM::Interface &m_field, boost::shared_ptr< FEMethod > &fe_top, ForcesAndSourcesCore::UserDataOperator::OpType op, int verbosity, LogManager::SeverityLevel sev)
set levels of projection operators, which project field data from parent entities,...
FTensor::Index< 'i', SPACE_DIM > i
constexpr char FIELD_NAME[]
constexpr int SPACE_DIM
constexpr int FIELD_DIM
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, FIELD_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainSource
Operator to integrate the right hand side matrix for the problem.
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temprary meshset.
Definition: Templates.hpp:1584
const double r
rate factor
constexpr double t
plate stiffness
Definition: plate.cpp:59
boost::shared_ptr< VectorDouble > approxVals
SmartPetscObj< Vec > L2Vec
SmartPetscObj< Vec > resVec
boost::shared_ptr< MatrixDouble > divApproxVals
Operator to evaluate errors.
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Simple * simpleInterface
MoFEMErrorCode checkResults()
[Check results]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
static ApproxFieldFunction< FIELD_DIM > approxFunction
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode readMesh()
[Run programme]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
static ApproxFieldFunctionDerivative< FIELD_DIM > divApproxFunction
MoFEM::Interface & mField
MoFEMErrorCode printResults()
[Check results]
MoFEMErrorCode runProblem()
[Run programme]
Managing BitRefLevels.
MoFEMErrorCode writeBitLevelByDim(const BitRefLevel bit, const BitRefLevel mask, const int dim, const char *file_name, const char *file_type, const char *options, const bool check_for_empty=true) const
Write bit ref level to file.
MoFEMErrorCode filterEntitiesByRefLevel(const BitRefLevel bit, const BitRefLevel mask, Range &ents, int verb=QUIET) const
filter entities by bit ref level
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =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
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
FieldApproximationBase & getBase()
Get approximation base.
const VectorFieldEntities & getFieldEntities() const
get field entities
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FieldSpace & getSpace()
Get field space.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPCOL
operator doWork function is executed on FE columns
@ OPROW
operator doWork function is executed on FE rows
@ OPSPACE
operator do Work is execute on space data
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure to get information form mofem into EntitiesFieldData
MatrixDouble gaussPts
Matrix of integration points.
Mesh refinement interface.
Operator to project base functions from parent entity.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.