v0.14.0
hanging_node_approx.cpp

Testing approximation with hanging nodes.

/**
* \example hanging_node_approx.cpp
*
* Testing 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
using BoundaryEle =
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
*
*/
PETSC>::LinearForm<GAUSS>::OpSource<1, FIELD_DIM>;
/**
* @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,
int verbosity, LogManager::SeverityLevel sev) {
BitRefLevel bit_marker;
for (auto l = 1; l <= nb_ref_levels; ++l)
bit_marker |= marker(l);
/**
* @brief Collect data from parent elements to child
*/
boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
add_parent_level =
[&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
// Evaluate if not last parent element
if (level > 0) {
// Create domain parent FE
auto fe_ptr_current = boost::shared_ptr<ForcesAndSourcesCore>(
new PARENT_FE(m_field));
if (op == DomainEleOp::OPSPACE) {
// Push base function
if (typeid(PARENT_FE) == typeid(DomainParentEle))
fe_ptr_current->getOpPtrVector(), {H1});
}
// Call next level
add_parent_level(
boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
fe_ptr_current),
level - 1);
// Add data to curent fe level
if (op == DomainEleOp::OPSPACE) {
// Only base
parent_fe_pt->getOpPtrVector().push_back(
H1, op, fe_ptr_current,
BitRefLevel().set(), bit(0).flip(),
bit_marker, BitRefLevel().set(),
verbosity, sev));
} else {
// Filed data
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, retrieve 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) {}
MoFEMErrorCode runProblem();
private:
Simple *simpleInterface;
static ApproxFieldFunction<FIELD_DIM> approxFunction;
static ApproxFieldFunctionDerivative<FIELD_DIM> divApproxFunction;
/**
* @brief red mesh and randomly refine three times
*
* @return MoFEMErrorCode
*/
MoFEMErrorCode readMesh();
/**
* @brief add field, and set up problem
*
* @return MoFEMErrorCode
*/
MoFEMErrorCode setupProblem();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode checkResults();
MoFEMErrorCode printResults();
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_w = getFTensor0IntegrationWeight();
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]
CHKERR readMesh();
CHKERR setupProblem();
CHKERR assembleSystem();
CHKERR solveSystem();
CHKERR checkResults();
CHKERR printResults();
}
//! [Run programme]
//! [Read mesh]
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
Skinner skin(&mField.get_moab());
CHKERR mField.getInterface(simpleInterface);
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 refine = mField.getInterface<MeshRefinement>();
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;
CHKERR simpleInterface->setFieldOrder(FIELD_NAME, order);
// Simple interface will resolve adjacency to DOFs of parent of the element.
// Using that information MAtrixManager allocate appropriately size of
// matrix.
simpleInterface->getParentAdjacencies() = true;
CHKERR simpleInterface->setUp();
// 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]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
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_NAME, DomainEleOp::OPROW);
field_op_row->doWorkRhsHook = [](DataOperator *op_ptr, int side,
EntityType type,
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();
for (auto &field_ent : data.getFieldEntities()) {
MOFEM_LOG("SELF", Sev::verbose)
<< "\t" << CN::EntityTypeName(field_ent->getEntType());
}
}
};
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
DomainEleOp::OPSPACE, NOISY, Sev::verbose);
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
DomainEleOp::OPROW, NOISY, Sev::noisy);
pipeline_mng->getOpDomainRhsPipeline().push_back(field_op_row);
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpDomainSource(FIELD_NAME, approxFunction));
}
//! [Push operators to pipeline]
//! [Solve]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
MOFEM_LOG("WORLD", Sev::inform) << "Solve problem";
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto dm = simpleInterface->getDM();
auto D = createDMVector(dm);
auto F = vectorDuplicate(D);
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]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
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 = createDMVector(simpleInterface->getDM());
common_data_ptr->L2Vec = createVectorMPI(
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>();
pipeline_mng->getOpDomainRhsPipeline(), {H1});
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(),
BoundaryEleOp::OPSPACE, QUIET, Sev::noisy);
set_parent_dofs<BoundaryParentEle>(mField, pipeline_mng->getBoundaryRhsFE(),
BoundaryEleOp::OPROW, VERBOSE, Sev::noisy);
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]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
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(
FIELD_NAME, DomainEleOp::OPROW);
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,
EntityType type,
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(),
DomainEleOp::OPSPACE, VERBOSE, Sev::noisy);
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]
}
}
AtomTest
Definition: child_and_parent.cpp:57
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
evaluate source, i.e. rhs vector
Definition: hanging_node_approx.cpp:69
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::EntitiesFieldData::EntData::getFieldEntities
const VectorFieldEntities & getFieldEntities() const
get field entities
Definition: EntitiesFieldData.hpp:1277
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MoFEM::PipelineManager::getDomainRhsFE
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Definition: PipelineManager.hpp:405
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
OpError
Definition: initial_diffusion.cpp:61
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition: DataOperators.hpp:24
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
NOISY
@ NOISY
Definition: definitions.h:224
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
ApproximationBaseNames
const static char *const ApproximationBaseNames[]
Definition: definitions.h:72
MoFEM::PipelineManager::loopFiniteElements
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
Definition: PipelineManager.cpp:19
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
EntData
EntitiesFieldData::EntData EntData
Definition: hanging_node_approx.cpp:29
MoFEM::PipelineManager::getOpDomainRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Definition: PipelineManager.hpp:773
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::ForcesAndSourcesCore::UserDataOperator::OpType
OpType
Controls loop over entities on element.
Definition: ForcesAndSourcesCore.hpp:566
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
MoFEM::DMoFEMMeshToLocalVector
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:523
MoFEM::PipelineManager::getBoundaryRhsFE
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Definition: PipelineManager.hpp:413
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
AtomTest::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition: child_and_parent.cpp:290
nb_ref_levels
constexpr int nb_ref_levels
Three levels of refinement.
Definition: hanging_node_approx.cpp:17
test_bit_child
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
Definition: hanging_node_approx.cpp:173
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
Definition: child_and_parent.cpp:53
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
sdf.r
int r
Definition: sdf.py:8
SPACE_DIM
constexpr int SPACE_DIM
Definition: hanging_node_approx.cpp:16
MoFEM::BitRefManager::filterEntitiesByRefLevel
MoFEMErrorCode filterEntitiesByRefLevel(const BitRefLevel bit, const BitRefLevel mask, Range &ents, int verb=QUIET) const
filter entities by bit ref level
Definition: BitRefManager.cpp:746
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
help
static char help[]
Definition: hanging_node_approx.cpp:12
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
AtomTest::OpError
Operator to evaluate errors.
Definition: child_and_parent.cpp:82
VERBOSE
@ VERBOSE
Definition: definitions.h:222
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
BoundaryEleOp
MoFEM::BitRefManager::getEntitiesByDimAndRefLevel
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
Definition: BitRefManager.cpp:822
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MoFEM::MeshRefinement
Mesh refinement interface.
Definition: MeshRefinement.hpp:26
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1886
FIELD_DIM
constexpr int FIELD_DIM
Definition: hanging_node_approx.cpp:15
double
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MoFEM::ProblemsManager::removeDofsOnEntities
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
Definition: ProblemsManager.cpp:2954
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
DomainParentEle
ElementsAndOps< SPACE_DIM >::DomianParentEle DomainParentEle
Definition: child_and_parent.cpp:35
AtomTest::divApproxFunction
static ApproxFieldFunctionDerivative< FIELD_DIM > divApproxFunction
Definition: hanging_node_approx.cpp:189
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::PipelineManager::createKSP
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
Definition: PipelineManager.cpp:49
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
MoFEM::BitRefManager::updateMeshsetByEntitiesChildren
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.
Definition: BitRefManager.cpp:1029
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::BitRefManager::addBitRefLevel
MoFEMErrorCode addBitRefLevel(const Range &ents, const BitRefLevel &bit, int verb=QUIET) const
add bit ref level to ref entity
Definition: BitRefManager.cpp:554
MoFEM::EntitiesFieldData::EntData::getBase
FieldApproximationBase & getBase()
Get approximation base.
Definition: EntitiesFieldData.hpp:1313
AtomTest::OpErrorSkel
Definition: hanging_node_approx.cpp:217
MoFEM::BitRefManager::writeBitLevelByDim
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.
Definition: BitRefManager.cpp:660
MoFEM::PipelineManager::setDomainRhsIntegrationRule
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:530
MoFEM::LogManager::SeverityLevel
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
AtomTest::approxFunction
static ApproxFieldFunction< FIELD_DIM > approxFunction
Definition: child_and_parent.cpp:67
MoFEM::PipelineManager::getOpDomainLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition: PipelineManager.hpp:749
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
MoFEM::ForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: ForcesAndSourcesCore.hpp:482
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
BiLinearForm
FTensor::Index< 'i', SPACE_DIM >
AtomTest::solveSystem
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
Definition: child_and_parent.cpp:377
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
MoFEM::PipelineManager::setDomainLhsIntegrationRule
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:503
MoFEM::BitRefManager::updateRangeByChildren
MoFEMErrorCode updateRangeByChildren(const Range &parent, Range &child, MoFEMTypes bh=MF_ZERO)
Update range by childrens.
Definition: BitRefManager.cpp:1144
main
int main(int argc, char *argv[])
Definition: hanging_node_approx.cpp:730
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
DomainEleOp
MoFEM::CoreTmp< 0 >::Initialize
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
FieldSpaceNames
const static char *const FieldSpaceNames[]
Definition: definitions.h:92
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1318
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
ApproxFieldFunctionDerivative
Definition: hanging_node_approx.cpp:32
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
BITREFLEVEL_SIZE
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:232
CommonData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::EntitiesFieldData::EntData::getSpace
FieldSpace & getSpace()
Get field space.
Definition: EntitiesFieldData.hpp:1315
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::ForcesAndSourcesCore::gaussPts
MatrixDouble gaussPts
Matrix of integration points.
Definition: ForcesAndSourcesCore.hpp:109
AtomTest::checkResults
MoFEMErrorCode checkResults()
[Check results]
Definition: dg_projection.cpp:213
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
evaluate mass matrix
Definition: hanging_node_approx.cpp:62
MoFEM::PipelineManager::getOpBoundaryRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
Definition: PipelineManager.hpp:821
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
AtomTest::readMesh
MoFEMErrorCode readMesh()
[Run programme]
Definition: child_and_parent.cpp:208
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::PipelineManager::setBoundaryRhsIntegrationRule
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:584
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
marker
auto marker
set bit to marker
Definition: hanging_node_approx.cpp:82
BoundaryParentEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
Definition: child_and_parent.cpp:41
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
QUIET
@ QUIET
Definition: definitions.h:221
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::SmartPetscObj< Vec >
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
Definition: child_and_parent.cpp:55
set_parent_dofs
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,...
Definition: hanging_node_approx.cpp:92
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::OpAddParentEntData
Operator to project base functions from parent entity to child.
Definition: MeshProjectionDataOperators.hpp:66
AtomTest::printResults
MoFEMErrorCode printResults()
[Check results]
Definition: hanging_node_approx.cpp:655
AtomTest::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: child_and_parent.cpp:268
AtomTest::runProblem
MoFEMErrorCode runProblem()
[Run programme]
Definition: child_and_parent.cpp:188
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
FIELD_NAME
constexpr char FIELD_NAME[]
Definition: hanging_node_approx.cpp:14
ApproxFieldFunction
Definition: child_and_parent.cpp:43
F
@ F
Definition: free_surface.cpp:394