v0.14.0
Loading...
Searching...
No Matches
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
*
*/
/**
* @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) {
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) {}
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;
// 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();
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(
}
//! [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 = 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]
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(),
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:311
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
@ NOISY
Definition: definitions.h:211
@ 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
constexpr int order
FTensor::Index< 'm', SPACE_DIM > m
@ F
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
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_deque< 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 addBitRefLevel(const Range &ents, const BitRefLevel &bit, int verb=QUIET) const
add bit ref level to ref entity
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 updateRangeByChildren(const Range &parent, Range &child, MoFEMTypes bh=MF_ZERO)
Update range by childrens.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
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.
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
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1767
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
int r
Definition: sdf.py:8
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]
Add operators pushing bases from local to physical configuration.
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 to child.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
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:27
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.