v0.14.0
test_cache_on_entities.cpp

Testing entities cache. Entities cache is data abstraction enabling user to pass data between operators, finite elements, and problems, in convenient way.

It can be used to store indices, or history variables, and any other data, which I can not imagine at that moment of time, and how to make use of it depends on user creativity and imagination.

/**
* \file test_cache_on_entities.cpp
* \example test_cache_on_entities.cpp
*
* Testing entities cache. Entities cache is data abstraction enabling user to
* pass data between operators, finite elements, and problems, in convenient
* way.
*
* It can be used to store indices, or history variables, and any other data,
* which I can not imagine at that moment of time, and how to make use of it
* depends on user creativity and imagination.
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
MyStorage(VectorInt &data) : globalIndices(data) {}
VectorInt globalIndices;
};
/**
* @brief Operator set cache stored data, in this example, global indices, but
* it can be any structure
*
*/
struct OpVolumeSet : public VolOp {
OpVolumeSet(const std::string &field_name) : VolOp(field_name, OPROW) {}
MoFEMErrorCode doWork(int side, EntityType type,
MOFEM_LOG_TAG("SYNC", "OpVolumeSet");
// Clear data when start process element
if (type == MBVERTEX)
entsIndices.clear();
// Get field entity pointer
auto field_ents = data.getFieldEntities();
if (auto e_ptr = field_ents[0]) {
// Add indices to global storage
entsIndices.push_back(boost::make_shared<MyStorage>(data.getIndices()));
// Store pointer to data on entity
e_ptr->getWeakStoragePtr() = entsIndices.back();
MOFEM_LOG("SYNC", Sev::inform) << "Set " << e_ptr->getEntTypeName() << " "
<< side << " : " << entsIndices.size();
// Check if all works
if (auto ptr = e_ptr->getSharedStoragePtr<EntityStorage>()) {
if (auto cast_ptr = boost::dynamic_pointer_cast<MyStorage>(ptr))
MOFEM_LOG("SYNC", Sev::noisy) << "Cast works";
else
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Cast not works");
} else {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Pointer to stored data is expected to be set");
}
if (auto ptr = e_ptr->getSharedStoragePtr<MyStorage>()) {
MOFEM_LOG("SYNC", Sev::verbose)
<< data.getIndices() << " : " << ptr->globalIndices;
} else {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Pointer to stored data is expected to be set");
}
} else {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Pointer to entity should be set by finite element");
}
}
using SharedVecInt = boost::shared_ptr<MyStorage>;
static std::vector<SharedVecInt>
entsIndices; ///< This is global static storage
};
std::vector<OpVolumeSet::SharedVecInt> OpVolumeSet::entsIndices;
/**
* @brief Test if cached data can be accessed, and check consistency of data
*
*/
struct OpVolumeTest : public VolOp {
OpVolumeTest(const std::string &field_name) : VolOp(field_name, OPROW) {}
MoFEMErrorCode doWork(int side, EntityType type,
MOFEM_LOG_TAG("SYNC", "OpVolumeTest");
// Get pointer to field entities
auto field_ents = data.getFieldEntities();
if (auto e_ptr = field_ents[0]) {
MOFEM_LOG("SYNC", Sev::inform)
<< "Test " << e_ptr->getEntTypeName() << " " << side;
// Check if data are cached on entity, and if code is correct, data should
// accessible.
if (auto ptr = e_ptr->getSharedStoragePtr<MyStorage>()) {
MOFEM_LOG("SYNC", Sev::verbose)
<< data.getIndices() << " : " << ptr->globalIndices;
// Check constancy of data. Stored data are indices, and expected stored
// that should be indices, thus difference between should be zero.
auto diff = data.getIndices() - ptr->globalIndices;
auto dot = inner_prod(diff, diff);
if (dot > 0)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Test failed");
} else {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Pointer to stored data is expected to be set");
}
} else {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Pointer to entity should be set by finite element");
}
}
};
struct OpVolumeAssemble : public VolOp {
OpVolumeAssemble(const std::string &field_name, Vec v)
: VolOp(field_name, OPROW), V(v) {}
MoFEMErrorCode doWork(int side, EntityType type,
auto nb_dofs = data.getIndices().size();
if (nb_dofs) {
std::vector<double> v(nb_dofs, 1);
CHKERR VecSetValues<EssentialBcStorage>(V, data, &v[0], ADD_VALUES);
}
}
private:
Vec V;
};
struct VolRule {
int operator()(int, int, int) const { return 1; }
};
int main(int argc, char *argv[]) {
// initialize petsc
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
// Create MoAB database
moab::Core moab_core;
moab::Interface &moab = moab_core;
// Create MoFEM database and link it to MoAB
MoFEM::Core mofem_core(moab);
MoFEM::Interface &m_field = mofem_core;
// Register DM Manager
DMType dm_name = "DMMOFEM";
// Simple interface
Simple *simple_interface;
CHKERR m_field.getInterface(simple_interface);
{
// get options from command line
CHKERR simple_interface->getOptions();
// load mesh file
CHKERR simple_interface->loadFile();
// add fields
CHKERR simple_interface->addDomainField("FIELD", H1,
// set fields order
CHKERR simple_interface->setFieldOrder("FIELD", 4);
// setup problem
CHKERR simple_interface->setUp();
// get dm
auto dm = simple_interface->getDM();
// create elements
auto domain_fe = boost::make_shared<VolEle>(m_field);
// set integration rule
domain_fe->getRuleHook = VolRule();
auto get_skin = [&]() {
Range ents;
CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, ents, true);
Skinner skin(&m_field.get_moab());
Range skin_ents;
CHKERR skin.find_skin(0, ents, false, skin_ents);
// filter not owned entities, those are not on boundary
ParallelComm *pcomm =
ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
Range proc_skin;
CHKERR pcomm->filter_pstatus(skin_ents,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &proc_skin);
Range adj;
for (auto d : {0, 1, 2})
CHKERR m_field.get_moab().get_adjacencies(proc_skin, d, false, adj,
moab::Interface::UNION);
// proc_skin.merge(adj);
return proc_skin;
};
auto get_mark_skin_dofs = [&](Range &&skin) {
auto problem_manager = m_field.getInterface<ProblemsManager>();
auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
problem_manager->markDofs(simple_interface->getProblemName(), ROW,
*marker_ptr);
return marker_ptr;
};
auto mark_dofs = get_mark_skin_dofs(get_skin());
// set operator to the volume elements
domain_fe->getOpPtrVector().push_back(new OpVolumeSet("FIELD"));
domain_fe->getOpPtrVector().push_back(new OpVolumeTest("FIELD"));
domain_fe->getOpPtrVector().push_back(
new OpSetBc("FIELD", true, mark_dofs));
domain_fe->getOpPtrVector().push_back(new OpVolumeAssemble("FIELD", v));
domain_fe->getOpPtrVector().push_back(new OpUnSetBc("FIELD"));
// make integration in volume (here real calculations starts)
domain_fe);
CHKERR VecAssemblyBegin(v);
CHKERR VecAssemblyEnd(v);
CHKERR VecGhostUpdateBegin(v, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(v, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(simple_interface->getDM(), v,
INSERT_VALUES, SCATTER_REVERSE);
// CHKERR VecView(v,PETSC_VIEWER_STDOUT_WORLD);
const MoFEM::Problem *problem_ptr;
CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
double *array;
CHKERR VecGetArray(v, &array);
for (auto dof : *(problem_ptr->getNumeredRowDofsPtr())) {
auto loc_idx = dof->getPetscLocalDofIdx();
if (dof->getFieldData() != array[loc_idx])
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Data inconsistency");
if ((*mark_dofs)[loc_idx]) {
if (array[loc_idx] != 0)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Dof assembled, but it should be not");
} else {
if (array[loc_idx] == 0)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Dof not assembled, but it should be");
}
}
CHKERR VecRestoreArray(v, &array);
}
}
// finish work cleaning memory, getting statistics, etc.
return 0;
}
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
MoFEM::EntitiesFieldData::EntData::getFieldEntities
const VectorFieldEntities & getFieldEntities() const
get field entities
Definition: EntitiesFieldData.hpp:1264
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM::ProblemsManager::OR
@ OR
Definition: ProblemsManager.hpp:365
OpVolumeSet
Operator set cache stored data, in this example, global indices, but it can be any structure.
Definition: test_cache_on_entities.cpp:37
MoFEM.hpp
VolRule
Set integration rule.
Definition: simple_interface.cpp:88
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
MyStorage
Definition: test_cache_on_entities.cpp:27
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
OpVolumeTest
Test if cached data can be accessed, and check consistency of data.
Definition: test_cache_on_entities.cpp:101
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1167
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:86
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::EntityStorage
Definition: FieldEntsMultiIndices.hpp:16
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
convert.type
type
Definition: convert.py:64
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:105
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:264
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::Problem::getNumeredRowDofsPtr
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Definition: ProblemsMultiIndices.hpp:82
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:341
MoFEM::ProblemsManager::markDofs
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const enum MarkOP op, const Range ents, std::vector< unsigned char > &marker) const
Create vector with marked indices.
Definition: ProblemsManager.cpp:3200
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
help
static char help[]
Definition: test_cache_on_entities.cpp:20
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
main
int main(int argc, char *argv[])
Definition: test_cache_on_entities.cpp:167
MoFEM::OpUnSetBc
Definition: FormsIntegrators.hpp:49
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
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
OpVolumeAssemble
Definition: test_cache_on_entities.cpp:145
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
MoFEM::Types::VectorInt
UBlasVector< int > VectorInt
Definition: Types.hpp:67
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::OpSetBc
Set indices on entities on finite element.
Definition: FormsIntegrators.hpp:38
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:426
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
OpVolumeSet::entsIndices
static std::vector< SharedVecInt > entsIndices
This is global static storage.
Definition: test_cache_on_entities.cpp:92
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::SmartPetscObj< Vec >
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:362
MoFEM::DMoFEMLoopFiniteElements
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:586
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MoFEM::VecSetValues< EssentialBcStorage >
MoFEMErrorCode VecSetValues< EssentialBcStorage >(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Set values to vector in operator.
Definition: FormsIntegrators.cpp:76