static char help[] =
"...\n\n";
constexpr bool debug =
false;
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "SPLIT"));
LogManager::setLog("SPLIT");
try {
PetscBool flg_file = PETSC_FALSE;
PetscBool squash_bit_levels = PETSC_TRUE;
PetscBool flg_list_of_sidesets = PETSC_FALSE;
PetscBool add_prisms = PETSC_FALSE;
PetscBool split_corner_edges = PETSC_FALSE;
int nb_sidesets = 10;
int sidesets[nb_sidesets];
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Split sides options",
"none");
CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
CHKERR PetscOptionsBool(
"-squash_bit_levels",
"squash bit levels",
"",
squash_bit_levels, &squash_bit_levels, NULL);
CHKERR PetscOptionsIntArray(
"-side_sets",
"get list of sidesets",
"",
sidesets, &nb_sidesets, &flg_list_of_sidesets);
CHKERR PetscOptionsBool(
"-output_vtk",
"if true outout vtk file",
"",
CHKERR PetscOptionsBool(
"-add_prisms",
"if true outout vtk file",
"",
add_prisms, &add_prisms, PETSC_NULL);
CHKERR PetscOptionsBool(
"-split_corner_edges",
"if true outout vtk file",
"", split_corner_edges, &split_corner_edges,
PETSC_NULL);
CHKERR PetscOptionsBool(
"-output_vtk",
"if true outout vtk file",
"",
ierr = PetscOptionsEnd();
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
const char *option;
option = "";
if (flg_file != PETSC_TRUE) {
"*** ERROR -my_file (mesh file needed)");
}
if (flg_list_of_sidesets != PETSC_TRUE) {
"List of sidesets not given -my_side_sets ...");
}
auto &meshsets_index = m_mng->getMeshsetsMultindex();
auto &m_by_id_and_type =
for (
auto mit = m_by_type.lower_bound(
SIDESET);
mit != m_by_type.upper_bound(
SIDESET); mit++) {
<< "Sideset on the mesh id = " << mit->getMeshsetId();
}
std::vector<BitRefLevel> bit_levels;
auto update_meshsets = [&]() {
CHKERR bit_mng->updateMeshsetByEntitiesChildren(
cubit_meshset, bit_levels.back(), cubit_meshset, MBMAXTYPE, true);
}
};
if (split_corner_edges) {
auto meshset_of_edges_to_refine_ptr = get_temp_meshset_ptr(moab);
for (int mm = 0; mm != nb_sidesets; mm++) {
auto mit =
m_by_id_and_type.find(boost::make_tuple(sidesets[mm],
SIDESET));
if (mit == m_by_id_and_type.end())
"No sideset in database id = %d", sidesets[mm]);
CHKERR moab.get_entities_by_type(mit->getMeshset(), MBTRI, faces,
true);
CHKERR moab.get_adjacencies(faces, 1,
true, faces_edges,
moab::Interface::UNION);
CHKERR skin.find_skin(0, faces,
false, skin_edges);
CHKERR moab.get_connectivity(skin_edges, skin_verts,
true);
CHKERR moab.get_adjacencies(skin_verts, 1,
true, vertex_edges,
moab::Interface::UNION);
Range vertex_edges_verts;
CHKERR moab.get_connectivity(vertex_edges, vertex_edges_verts,
true);
vertex_edges_verts = subtract(vertex_edges_verts, skin_verts);
Range vertex_edges_verts_edges;
CHKERR moab.get_adjacencies(vertex_edges_verts, 1,
true,
vertex_edges_verts_edges,
moab::Interface::UNION);
vertex_edges = subtract(vertex_edges, vertex_edges_verts_edges);
vertex_edges = subtract(vertex_edges, skin_edges);
vertex_edges = intersect(vertex_edges, faces_edges);
CHKERR moab.add_entities(*meshset_of_edges_to_refine_ptr, vertex_edges);
}
int nb_tris;
CHKERR moab.get_number_entities_by_type(*meshset_of_edges_to_refine_ptr,
MBEDGE, nb_tris, true);
MOFEM_LOG(
"SPLIT", Sev::inform) <<
"Refine corner edges " << nb_tris;
CHKERR moab.write_file(
"out_edges_to_refine.vtk",
"VTK",
"",
meshset_of_edges_to_refine_ptr->get_ptr(), 1);
CHKERR refine_mng->addVerticesInTheMiddleOfEdges(
*meshset_of_edges_to_refine_ptr, bit_levels.back());
CHKERR refine_mng->refineTets(0, bit_levels.back());
}
for (int mm = 0; mm != nb_sidesets; mm++) {
auto mit =
m_by_id_and_type.find(boost::make_tuple(sidesets[mm],
SIDESET));
if (mit == m_by_id_and_type.end())
"No sideset in database id = %d", sidesets[mm]);
<< "Split sideset " << mit->getMeshsetId();
const auto cubit_meshset = mit->getMeshset();
{
auto ref_level_meshset_ptr = get_temp_meshset_ptr(moab);
CHKERR bit_mng->getEntitiesByTypeAndRefLevel(bit_levels.back(),
*ref_level_meshset_ptr);
CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
*ref_level_meshset_ptr);
CHKERR moab.get_entities_by_handle(*ref_level_meshset_ptr,
ref_level_tets, true);
CHKERR interface_ptr->getSides(cubit_meshset, bit_levels.back(),
true,
0);
<< "Add bit level " << bit_levels.size();
bit_levels.push_back(
BitRefLevel().set(bit_levels.size()));
CHKERR interface_ptr->splitSides(*ref_level_meshset_ptr,
bit_levels.back(), cubit_meshset,
add_prisms, true, 0);
}
}
if (squash_bit_levels == PETSC_TRUE) {
for (unsigned int ll = 0; ll != bit_levels.size() - 1; ll++) {
true);
}
CHKERR bit_mng->shiftRightBitRef(bit_levels.size() - 1);
}
auto meshset_ptr = get_temp_meshset_ptr(moab);
if (squash_bit_levels)
else
MBTET, *meshset_ptr);
CHKERR moab.write_file(
"out.vtk",
"VTK",
"", meshset_ptr->get_ptr(), 1);
}
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
if (pcomm == NULL)
"Communicator should be allocated");
CHKERR pcomm->assign_global_ids(0, 3, 1,
false);
CHKERR moab.write_file(
"out.h5m");
}
return 0;
}
static PetscErrorCode ierr
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
#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_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
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
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode delete_ents_by_bit_ref(const BitRefLevel bit, const BitRefLevel mask, const bool remove_parent=false, int verb=DEFAULT_VERBOSITY, MoFEMTypes mf=MF_ZERO)=0
delete entities form mofem and moab database
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.
Mesh refinement interface.
Interface for managing meshsets containing materials and boundary conditions.
Create interface from given surface and insert flat prisms in-between.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.