v0.13.1
forces_and_sources_testing_users_base.cpp

Primarily this is used for testing if the code can handle user base. It is also, an example of how to build and use user approximation base. This is a test, so we used RT base by Demkowicz recipe.

Note that triple defines approximation element; element entity type, approximation space and approximation base. Entity type determines the integration method; approximation space determines the adjacency of the matrix and approximation base determines together with space the regularity of approximation.

/** \file forces_and_sources_testing_users_base.cpp
* \example forces_and_sources_testing_users_base.cpp
*
* Primarily this is used for testing if the code can handle user base. It is
* also, an example of how to build and use user approximation base. This is a
* test, so we used RT base by Demkowicz recipe.
*
* Note that triple defines approximation element; element entity type,
* approximation space and approximation base. Entity type determines the
* integration method; approximation space determines the adjacency of the
* matrix and approximation base determines together with space the regularity
* of approximation.
*
*/
#include <MoFEM.hpp>
#include <Hdiv.hpp>
namespace bio = boost::iostreams;
using bio::stream;
using bio::tee_device;
using namespace MoFEM;
static char help[] = "...\n\n";
/**
* @brief Class used to calculate base functions at integration points
*
*/
/**
* @brief Return interface to this class when one ask for for tetrahedron,
* otherisw return interface class for generic class.
*
* @param iface interface class
* @return MoFEMErrorCode
*/
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
UnknownInterface **iface) const {
*iface = const_cast<SomeUserPolynomialBase *>(this);
}
/**
* @brief Calculate base functions at intergeneration points
*
* @param pts
* @param ctx_ptr
* @return MoFEMErrorCode
*/
boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
int nb_gauss_pts = pts.size2();
if (!nb_gauss_pts) {
}
if (pts.size1() < 3) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong dimension of pts, should be at least 3 rows with "
"coordinates");
}
switch (cTx->sPace) {
case HDIV:
break;
default:
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
}
}
private:
const FieldApproximationBase base = cTx->bAse;
// This should be used only in case USER_BASE is selected
if (cTx->bAse != USER_BASE) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong base, should be USER_BASE");
}
// This is example, simply use Demkowicz HDiv base to generate base
// functions
EntitiesFieldData &data = cTx->dAta;
int nb_gauss_pts = pts.size2();
// calculate shape functions, i.e. barycentric coordinates
shapeFun.resize(nb_gauss_pts, 4, false);
CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
&pts(2, 0), nb_gauss_pts);
// derivatives of shape functions
double diff_shape_fun[12];
CHKERR ShapeDiffMBTET(diff_shape_fun);
int volume_order = data.dataOnEntities[MBTET][0].getOrder();
int p_f[4];
double *phi_f[4];
double *diff_phi_f[4];
// Calculate base function on tet faces
for (int ff = 0; ff != 4; ff++) {
int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
int order = volume_order > face_order ? volume_order : face_order;
data.dataOnEntities[MBTRI][ff].getN(base).resize(
nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
p_f[ff] = order;
phi_f[ff] = &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
diff_phi_f[ff] =
&*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
continue;
&data.facesNodes(ff, 0), order, &*shapeFun.data().begin(),
diff_shape_fun, phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
}
// Calculate base functions in tet interior
if (NBVOLUMETET_DEMKOWICZ_HDIV(volume_order) > 0) {
data.dataOnEntities[MBTET][0].getN(base).resize(
nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
data.dataOnEntities[MBTET][0].getDiffN(base).resize(
nb_gauss_pts, 9 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
double *phi_v = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
double *diff_phi_v =
&*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
volume_order, &*shapeFun.data().begin(), diff_shape_fun, p_f, phi_f,
diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
}
// Set size of face base correctly
for (int ff = 0; ff != 4; ff++) {
int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
data.dataOnEntities[MBTRI][ff].getN(base).resize(
nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
}
}
};
int main(int argc, char *argv[]) {
// Initialise MoFEM, MPI and petsc
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
// create moab
moab::Core mb_instance;
// get interface to moab databse
moab::Interface &moab = mb_instance;
// get file
PetscBool flg = PETSC_TRUE;
char mesh_file_name[255];
#if PETSC_VERSION_GE(3, 6, 4)
CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
255, &flg);
#else
CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
mesh_file_name, 255, &flg);
#endif
if (flg != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"*** ERROR -my_file (MESH FILE NEEDED)");
}
// create MoFEM database
MoFEM::Core core(moab);
// get interface to moab database
MoFEM::Interface &m_field = core;
// load mesh file
const char *option;
option = "";
CHKERR moab.load_file(mesh_file_name, 0, option);
// set bit refinement level
CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
0, 3, BitRefLevel().set(0));
// Create fields, field "FIELD_CGG" has user base, it means that recipe how
// to construct approximation is provided by user. Is set that user provided
// base is in h-div space.
CHKERR m_field.add_field("FILED_CGG", HDIV, USER_BASE, 1);
CHKERR m_field.add_field("FILED_RT", HDIV, DEMKOWICZ_JACOBI_BASE, 1);
// get access to "FIELD_CGG" data structure
auto field_ptr = m_field.get_field_structure("FILED_CGG");
// get table associating number of dofs to entities depending on
// approximation order set on those entities.
auto field_order_table =
const_cast<Field *>(field_ptr)->getFieldOrderTable();
// function set zero number of dofs
auto get_cgg_bubble_order_zero = [](int p) { return 0; };
// function set non-zero number of dofs on tetrahedrons
auto get_cgg_bubble_order_face = [](int p) {
};
auto get_cgg_bubble_order_tet = [](int p) {
};
field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
field_order_table[MBTRI] = get_cgg_bubble_order_face;
field_order_table[MBTET] = get_cgg_bubble_order_tet;
const_cast<Field *>(field_ptr)->rebuildDofsOrderMap();
auto &dof_order_map = field_ptr->getDofOrderMap(MBTET);
for(auto d = 0; d!=10; ++d) {
MOFEM_LOG("WORLD", Sev::noisy) << "dof " << dof_order_map[d];
}
CHKERR m_field.add_finite_element("FE");
// define rows/cols and element data
CHKERR m_field.modify_finite_element_add_field_row("FE", "FILED_CGG");
CHKERR m_field.modify_finite_element_add_field_col("FE", "FILED_CGG");
CHKERR m_field.modify_finite_element_add_field_data("FE", "FILED_CGG");
CHKERR m_field.modify_finite_element_add_field_row("FE", "FILED_RT");
CHKERR m_field.modify_finite_element_add_field_col("FE", "FILED_RT");
CHKERR m_field.modify_finite_element_add_field_data("FE", "FILED_RT");
// add problem
CHKERR m_field.add_problem("PROBLEM");
// set finite elements for problem
CHKERR m_field.modify_problem_add_finite_element("PROBLEM", "FE");
// set refinement level for problem
BitRefLevel().set(0));
// meshset consisting all entities in mesh
EntityHandle root_set = moab.get_root_set();
// add entities to field
CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "FILED_CGG");
CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "FILED_RT");
// add entities to finite element
CHKERR m_field.add_ents_to_finite_element_by_type(root_set, MBTET, "FE");
// set app. order
int order = 3;
CHKERR m_field.set_field_order(root_set, MBTRI, "FILED_CGG", order);
CHKERR m_field.set_field_order(root_set, MBTET, "FILED_CGG", order);
CHKERR m_field.set_field_order(root_set, MBTRI, "FILED_RT", order);
CHKERR m_field.set_field_order(root_set, MBTET, "FILED_RT", order);
/****/
// build database
// build field
CHKERR m_field.build_fields();
// build finite elemnts
// build adjacencies
// build problem
CHKERR m_field.getInterface<ProblemsManager>()->buildProblem("PROBLEM",
true);
// dofs partitioning
CHKERR m_field.getInterface<ProblemsManager>()->partitionSimpleProblem(
"PROBLEM");
CHKERR m_field.getInterface<ProblemsManager>()->partitionFiniteElements(
"PROBLEM");
// what are ghost nodes, see Petsc Manual
CHKERR m_field.getInterface<ProblemsManager>()->partitionGhostDofs(
"PROBLEM");
typedef tee_device<std::ostream, std::ofstream> TeeDevice;
typedef stream<TeeDevice> TeeStream;
std::ofstream ofs("forces_and_sources_testing_users_base.txt");
TeeDevice my_tee(std::cout, ofs);
TeeStream my_split(my_tee);
/**
* Simple user data operator which main purpose is to print values
* of base functions at intergation points.
*
*/
TeeStream &my_split;
MyOp1(const std::string &row_filed, const std::string &col_field,
TeeStream &_my_split, char type)
row_filed, col_field, type),
my_split(_my_split) {
sYmm = false;
}
MoFEMErrorCode doWork(int side, EntityType type,
if (data.getIndices().empty()) {
}
my_split << rowFieldName << endl;
my_split << "side: " << side << " type: " << type << std::endl;
my_split << data << endl;
my_split << data.getN() << endl;
my_split << endl;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
if (row_data.getIndices().empty())
if (col_data.getIndices().empty())
my_split << rowFieldName << " : " << colFieldName << endl;
my_split << "row side: " << row_side << " row_type: " << row_type
<< std::endl;
my_split << "col side: " << col_side << " col_type: " << col_type
<< std::endl;
my_split << row_data.getIndices().size() << " : "
<< col_data.getIndices().size() << endl;
my_split << endl;
}
};
// create finite element instance
// set class needed to cinstruct user approximation base
fe1.getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new SomeUserPolynomialBase());
// push user data oprators
fe1.getOpPtrVector().push_back(
new MyOp1("FILED_CGG", "FILED_CGG", my_split,
ForcesAndSourcesCore::UserDataOperator::OPROW));
fe1.getOpPtrVector().push_back(
new MyOp1("FILED_CGG", "FILED_RT", my_split,
ForcesAndSourcesCore::UserDataOperator::OPROWCOL));
// iterate over finite elements, and execute user data operators on each
// of them
CHKERR m_field.loop_finite_elements("PROBLEM", "FE", fe1);
}
return 0;
}
static Index< 'p', 3 > p
Implementation of H-curl base function.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#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
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:306
tee_device< std::ostream, std::ofstream > TeeDevice
int main(int argc, char *argv[])
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual Field * get_field_structure(const std::string &name)=0
get field structure
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
#define NBFACETRI_DEMKOWICZ_HDIV(P)
char mesh_file_name[255]
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:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition: Hdiv.cpp:617
MoFEMErrorCode Hdiv_Demkowicz_Interior_MBTET(int p, double *N, double *diffN, int p_face[], double *phi_f[4], double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim)
Definition: Hdiv.cpp:762
CoreTmp< 0 > Core
Definition: Core.hpp:1086
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
Base class if inherited used to calculate base functions.
Managing BitRefLevels.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(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_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
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
Deprecated interface functions.
Class used to pass element data to calculate base functions on tet,triangle,edge.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
data structure for finite element entity
MatrixInt facesNodes
nodes on finite element faces
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Provide data structure for (tensor) field approximation.
const std::array< ApproximationOrder, MAX_DOFS_ON_ENTITY > & getDofOrderMap(const EntityType type) const
get hash-map relating dof index on entity with its order
Problem manager is used to build and partition problems.
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Class used to calculate base functions at integration points.
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts)
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Calculate base functions at intergeneration points.
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Return interface to this class when one ask for for tetrahedron, otherisw return interface class for ...
~SomeUserPolynomialBase()=default
SomeUserPolynomialBase()=default