17namespace bio = boost::iostreams;
23static char help[] =
"...\n\n";
25static const double face_coords[4][9] = {{0, 0, 0, 1, 0, 0, 0, 0, 1},
26 {1, 0, 0, 0, 1, 0, 0, 0, 1},
27 {0, 0, 0, 0, 1, 0, 0, 0, 1},
28 {0, 0, 0, 1, 0, 0, 0, 1, 0}};
31 {0, 0, 0, 1, 0, 0}, {1, 0, 0, 0, 1, 0}, {0, 1, 0, 0, 0, 0},
32 {0, 0, 0, 0, 0, 1}, {1, 0, 0, 0, 0, 1}, {0, 1, 0, 0, 0, 1}};
34int main(
int argc,
char *argv[]) {
43 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
45 PetscBool flg = PETSC_TRUE;
47#if PETSC_VERSION_GE(3, 6, 4)
56 if (flg != PETSC_TRUE) {
57 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
70 EntityHandle meshset_level0;
71 CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
96 "MESH_NODE_POSITIONS");
103 "MESH_NODE_POSITIONS");
110 "MESH_NODE_POSITIONS");
117 "MESH_NODE_POSITIONS");
134 EntityHandle root_set = moab.get_root_set();
150 CHKERR skin.find_skin(0, tets,
false, skin_faces);
158 faces = subtract(faces, skin_faces);
162 CHKERR moab.get_adjacencies(faces, 1,
false, edges, moab::Interface::UNION);
166 CHKERR moab.get_adjacencies(skin_faces, 1,
false, skin_edges,
167 moab::Interface::UNION);
219 CHKERR VecSetRandom(
v, PETSC_NULL);
222 "TEST_PROBLEM",
ROW,
v, INSERT_VALUES, SCATTER_REVERSE);
227 typedef tee_device<std::ostream, std::ofstream>
TeeDevice;
229 std::ofstream ofs(
"forces_and_sources_hdiv_continuity_check.txt");
241 "H1", UserDataOperator::OPROW),
242 mField(m_field), tH(
th) {}
253 boost::shared_ptr<const NumeredEntFiniteElement> mofem_fe =
254 getNumeredEntFiniteElementPtr();
256 EntityHandle face = side_table.get<1>()
257 .find(boost::make_tuple(
type, side))
260 int sense = side_table.get<1>()
261 .find(boost::make_tuple(
type, side))
266 int nb_dofs = data.
getN().size2();
267 for (
int dd = 0;
dd < nb_dofs;
dd++) {
273 (
const void **)&t_ptr);
278 if (
type == MBEDGE) {
280 boost::shared_ptr<const NumeredEntFiniteElement> mofem_fe =
281 getNumeredEntFiniteElementPtr();
283 EntityHandle edge = side_table.get<1>()
284 .find(boost::make_tuple(
type, side))
288 CHKERR mField.
get_moab().get_adjacencies(&edge, 1, 3,
false, adj_tets,
289 moab::Interface::UNION);
290 const int nb_adj_tets = adj_tets.size();
293 int nb_dofs = data.
getN().size2();
294 for (
int dd = 0;
dd < nb_dofs;
dd++) {
300 (
const void **)&t_ptr);
302 *t_ptr +=
t / nb_adj_tets;
326 gaussPts.resize(4, 4 + 6);
328 for (; ff < 4; ff++) {
330 for (;
dd < 3;
dd++) {
338 for (; ee < 6; ee++) {
340 for (;
dd < 3;
dd++) {
349 }
catch (std::exception &ex) {
350 std::ostringstream ss;
351 ss <<
"thorw in method: " << ex.what() <<
" at line " << __LINE__
352 <<
" in file " << __FILE__;
360 struct OpFacesSkinFluxes
370 "H1", UserDataOperator::OPROW),
371 mField(m_field), tH1(th1), tH2(th2), mySplit(my_split) {}
379 EntityHandle face = getNumeredEntFiniteElementPtr()->getEnt();
383 (
const void **)&t_ptr);
387 (
const void **)&tn_ptr);
393 int nb_dofs = data.
getN().size2();
394 for (;
dd < nb_dofs;
dd++) {
396 *tn_ptr += -data.
getN()(0,
dd) * val;
399 const double eps = 1e-8;
400 if (fabs(*tn_ptr) >
eps) {
402 "H1 continuity failed %6.4e", *tn_ptr);
405 mySplit.precision(5);
406 mySplit << face <<
" " << fabs(*tn_ptr) << std::endl;
422 "H1", UserDataOperator::OPROW),
423 mField(m_field), tH1(_th1), tH2(_th2), mySplit(my_split) {}
432 EntityHandle face = getNumeredEntFiniteElementPtr()->getEnt();
436 (
const void **)&t_ptr);
440 (
const void **)&tn_ptr);
445 const double eps = 1e-8;
446 if (fabs(*tn_ptr) >
eps) {
448 "H1 continuity failed %6.4e", *tn_ptr);
451 mySplit.precision(5);
453 mySplit << face <<
" " << fabs(*tn_ptr) << std::endl;
468 gaussPts.resize(3, 1);
487 "H1", UserDataOperator::OPROW),
488 mField(m_field), tH1(_th1), tH2(_th2), mySplit(my_split) {}
497 EntityHandle edge = getNumeredEntFiniteElementPtr()->getEnt();
501 (
const void **)&t_ptr);
505 (
const void **)&tn_ptr);
511 int nb_dofs = data.
getN().size2();
513 for (;
dd < nb_dofs;
dd++) {
515 tn += data.
getN()(0,
dd) * val;
535 mySplit.precision(5);
537 mySplit << edge <<
" " << fabs(*tn_ptr) << std::endl;
552 gaussPts.resize(2, 1);
560 MyTetFE tet_fe(m_field);
561 MyTriFE tri_fe(m_field);
562 MyTriFE skin_fe(m_field);
563 MyEdgeFE edge_fe(m_field);
566 double def_val[] = {0, 0, 0};
567 CHKERR moab.tag_get_handle(
"T", 3, MB_TYPE_DOUBLE, th1,
568 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
570 tet_fe.getOpPtrVector().push_back(
572 tet_fe.getOpPtrVector().push_back(
new OpTetFluxes(m_field, th1));
575 CHKERR moab.tag_get_handle(
"TN", 1, MB_TYPE_DOUBLE, th2,
576 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
578 tri_fe.getOpPtrVector().push_back(
580 tri_fe.getOpPtrVector().push_back(
582 tri_fe.getOpPtrVector().push_back(
584 skin_fe.getOpPtrVector().push_back(
585 new OpFacesSkinFluxes(m_field, th1, th2, my_split));
587 edge_fe.getOpPtrVector().push_back(
589 edge_fe.getOpPtrVector().push_back(
591 edge_fe.getOpPtrVector().push_back(
592 new OpEdgesFluxes(m_field, th1, th2, my_split));
594 for (Range::iterator fit = faces.begin(); fit != faces.end(); fit++) {
595 CHKERR moab.tag_set_data(th1, &*fit, 1, &def_val);
597 CHKERR moab.tag_set_data(th2, &*fit, 1, &def_val);
602 my_split <<
"internal\n";
604 my_split <<
"skin\n";
606 my_split <<
"edges\n";
610 EntityHandle meshset;
611 CHKERR moab.create_meshset(MESHSET_SET, meshset);
615 CHKERR moab.write_file(
"out.vtk",
"VTK",
"", &meshset, 1);
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_STD_EXCEPTION_THROW
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
int main(int argc, char *argv[])
static const double face_coords[4][9]
static const double edge_coords[6][6]
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
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 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.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
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.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition 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
double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
constexpr double t
plate stiffness
virtual moab::Interface & get_moab()=0
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.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
default operator for EDGE element
EdgeElementForcesAndSourcesCore(Interface &m_field)
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 VectorDouble & getFieldData() const
get dofs values
default operator for TRI element
FaceElementForcesAndSourcesCore(Interface &m_field)
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
MatrixDouble gaussPts
Matrix of integration points.
Calculate HO coordinates at gauss points.
Calculate normals at Gauss points of triangle element.
Calculate tangent vector on edge form HO geometry approximation.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Volume finite element base.
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpFacesFluxes(double &div)