v0.10.0
test_cache_on_entities.cpp

Tetsing entities cache. Entities acache is data abstraction enabling user to pass data between operators, finite elements, and problems, in convinient 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
*
* Tetsing entities cache. Entities acache is data abstraction enabling user to
* pass data between operators, finite elements, and problems, in convinient
* 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.
*
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
struct MyStorage : EntityStorage {
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<bool>>();
problem_manager->markDofs(simple_interface->getProblemName(), ROW, skin,
*marker_ptr);
return marker_ptr;
};
SmartPetscObj<Vec> v;
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)
CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
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;
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:349
MoFEM::FaceElementForcesAndSourcesCore FaceEle
@ ROW
Definition: definitions.h:192
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
@ H1
continuous field
Definition: definitions.h:177
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:292
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:132
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
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:445
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:348
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1026
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
MoFEMErrorCode VecSetValues< EssentialBcStorage >(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Set values to vector in operator.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:312
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:343
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:288
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
CoreTmp< 0 > Core
Definition: Core.hpp:1129
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
DataForcesAndSourcesCore::EntData EntData
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition: Core.hpp:77
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
Deprecated interface functions.
keeps basic data about problem
boost::shared_ptr< NumeredDofEntity_multiIndex > & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
Operator set cache stored data, in this example, global indices, but it can be any structure.
static std::vector< SharedVecInt > entsIndices
This is global static storage.
Test if cached data can be accessed, and check consistency of data.
Set integration rule.
int main(int argc, char *argv[])
static char help[]
FaceEle::UserDataOperator FaceOp