Testing integration on volume side on contact elementChecking continuity of H1 ( later hdiv and hcurl spaces on faces), and testing methods for integration on volume side on contact element.
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
try {
int rank;
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
PetscBool flg = PETSC_TRUE;
#if PETSC_VERSION_GE(3, 6, 4)
255, &flg);
#else
#endif
if (flg != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
const char *option;
option = "";
enum side_contact { MASTERSIDE, SLAVESIDE };
enum bases { AINSWORTH, DEMKOWICZ, LASTBASEOP };
const char *list_bases[] = {"ainsworth", "demkowicz"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
LASTBASEOP, &choice_base_value, &flg);
if (flg == PETSC_TRUE) {
if (choice_base_value == AINSWORTH)
else if (choice_base_value == DEMKOWICZ)
return base;
}
};
const char *list_spaces[] = {"h1", "hdiv", "hcurl"};
PetscBool flg;
PetscInt choice_space_value =
HDIV;
LASTSPACEOP, &choice_space_value, &flg);
if (flg == PETSC_TRUE) {
if (choice_space_value ==
HDIV)
else if (choice_space_value ==
HCURL)
else if (choice_space_value ==
H1)
return space;
}
};
auto add_prism_interface = [&](
Range &tets,
Range &prisms,
std::vector<BitRefLevel> &bit_levels) {
int ll = 1;
if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Insert %s (id: %d)\n",
cit->getName().c_str(), cit->getMeshsetId());
CHKERR moab.get_entities_by_type(cubit_meshset, MBTRI, tris,
true);
master_tris.merge(tris);
{
CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
->getEntitiesByTypeAndRefLevel(bit_levels.back(),
ref_level_meshset);
->getEntitiesByTypeAndRefLevel(bit_levels.back(),
ref_level_meshset);
CHKERR moab.get_entities_by_handle(ref_level_meshset,
ref_level_tets, true);
0);
cubit_meshset, true, true, 0);
CHKERR moab.delete_entities(&ref_level_meshset, 1);
}
->updateMeshsetByEntitiesChildren(
cubit_meshset, bit_levels.back(), cubit_meshset, MBMAXTYPE,
true);
}
}
}
for (unsigned int ll = 0; ll != bit_levels.size() - 1; ll++) {
true);
}
bit_levels.size() - 1);
CHKERR moab.create_meshset(MESHSET_SET, meshset_tets);
CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
->getEntitiesByTypeAndRefLevel(bit_levels[0],
BitRefLevel().set(),
MBTET, meshset_tets);
CHKERR moab.get_entities_by_handle(meshset_tets, tets,
true);
->getEntitiesByTypeAndRefLevel(bit_levels[0],
BitRefLevel().set(),
MBPRISM, meshset_prisms);
CHKERR moab.get_entities_by_handle(meshset_prisms, prisms);
CHKERR moab.get_adjacencies(prisms, 2,
false, tris,
moab::Interface::UNION);
tris = tris.subset_by_type(MBTRI);
slave_tris = subtract(tris, master_tris);
};
Range all_tets, contact_prisms, master_tris, slave_tris;
std::vector<BitRefLevel> bit_levels;
0, 3, bit_levels[0]);
CHKERR add_prism_interface(all_tets, contact_prisms, master_tris,
slave_tris, meshset_tets, meshset_prisms,
bit_levels);
}
};
struct OnContactSideMaster
elemData(elem_data) {}
if (
type == MBTRI && side == getFaceSideNumber()) {
getCoordsAtGaussPts() - getMasterCoordsAtGaussPts();
constexpr
double eps = 1e-12;
if (norm_inf(diff) >
eps)
"coordinates at integration pts are different");
const size_t nb_dofs = data.
getN().size2() / 3;
const size_t nb_integration_pts = data.
getN().size1();
ptr_elem_data = &elemData.shapeFunH1VolSide;
elem_data.resize(nb_integration_pts, nb_dofs, false);
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
for (size_t bb = 0; bb != nb_dofs; ++bb) {
elem_data(gg, bb) = t_base;
++t_base;
}
}
}
}
}
};
FACEMASTER),
volSideFe(m_field), elemData(elem_data) {
new OnContactSideMaster::OpVolSide(elemData));
}
if (
type == MBTRI && side == 0) {
const size_t nb_dofs = data.
getN().size2() / 3;
const size_t nb_integration_pts = data.
getN().size1();
elemData.shapeFunH1Values.resize(nb_integration_pts, nb_dofs,
false);
elemData.shapeFunH1VolSide.resize(nb_integration_pts, nb_dofs,
false);
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
for (size_t bb = 0; bb != nb_dofs; ++bb) {
elemData.shapeFunH1Values(gg, bb) = t_base;
++t_base;
}
}
}
std::string side_fe_name = "V1";
CHKERR loopSideVolumes(side_fe_name, volSideFe, 3, tri_master);
auto check_continuity_of_base = [&](auto &vol_dot_data) {
if (vol_dot_data.size1() != elemData.dotNormalFace.size1())
"Inconsistent number of integration points");
if (vol_dot_data.size2() != elemData.dotNormalFace.size2())
"Inconsistent number of base functions");
constexpr
double eps = 1e-12;
for (size_t gg = 0; gg != vol_dot_data.size1(); ++gg)
for (size_t bb = 0; bb != vol_dot_data.size2(); ++bb) {
const double error = std::abs(vol_dot_data(gg, bb) -
elemData.dotNormalFace(gg, bb));
"Inconsistency %3.4e != %3.4e", vol_dot_data(gg, bb),
elemData.dotNormalFace(gg, bb));
}
};
auto check_continuity_of_h1_base = [&](auto &vol_data) {
if (vol_data.size1() != elemData.shapeFunH1Values.size1())
"Inconsistent number of integration points");
if (vol_data.size2() != elemData.shapeFunH1Values.size2())
"Inconsistent number of base functions");
constexpr
double eps = 1e-12;
for (size_t gg = 0; gg != vol_data.size1(); ++gg)
for (size_t bb = 0; bb != vol_data.size2(); ++bb) {
const double error = std::abs(
vol_data(gg, bb) - elemData.shapeFunH1Values(gg, bb));
"Inconsistency %3.4e != %3.4e", vol_data(gg, bb),
elemData.shapeFunH1Values(gg, bb));
}
};
if (elemData.dotNormalEleLeft.size2() != 0)
CHKERR check_continuity_of_base(elemData.dotNormalEleLeft);
else if (elemData.dotNormalEleRight.size2() != 0)
CHKERR check_continuity_of_base(elemData.dotNormalEleRight);
else if (elemData.shapeFunH1VolSide.size2() != 0)
CHKERR check_continuity_of_h1_base(elemData.shapeFunH1VolSide);
}
}
};
struct OnContactSideSlave
elemData(elem_data) {}
if (
type == MBTRI && side == getFaceSideNumber()) {
getCoordsAtGaussPts() - getMasterCoordsAtGaussPts();
constexpr
double eps = 1e-12;
if (norm_inf(diff) >
eps)
"coordinates at integration pts are different");
const size_t nb_dofs = data.
getN().size2() / 3;
const size_t nb_integration_pts = data.
getN().size1();
ptr_elem_data = &elemData.shapeFunH1VolSide;
elem_data.resize(nb_integration_pts, nb_dofs, false);
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
for (size_t bb = 0; bb != nb_dofs; ++bb) {
elem_data(gg, bb) = t_base;
++t_base;
}
}
}
}
}
};
FACESLAVE),
volSideFe(m_field), elemData(elem_data) {
new OnContactSideSlave::OpVolSide(elemData));
}
if (
type == MBTRI && side == 0) {
const size_t nb_dofs = data.
getN().size2() / 3;
const size_t nb_integration_pts = data.
getN().size1();
elemData.shapeFunH1Values.resize(nb_integration_pts, nb_dofs,
false);
elemData.shapeFunH1VolSide.resize(nb_integration_pts, nb_dofs,
false);
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
for (size_t bb = 0; bb != nb_dofs; ++bb) {
elemData.shapeFunH1Values(gg, bb) = t_base;
++t_base;
}
}
}
std::string side_fe_name = "V1";
CHKERR loopSideVolumes(side_fe_name, volSideFe, 3, tri_slave);
auto check_continuity_of_base = [&](auto &vol_dot_data) {
if (vol_dot_data.size1() != elemData.dotNormalFace.size1())
"Inconsistent number of integration points");
if (vol_dot_data.size2() != elemData.dotNormalFace.size2())
"Inconsistent number of base functions");
constexpr
double eps = 1e-12;
for (size_t gg = 0; gg != vol_dot_data.size1(); ++gg)
for (size_t bb = 0; bb != vol_dot_data.size2(); ++bb) {
const double error = std::abs(vol_dot_data(gg, bb) -
elemData.dotNormalFace(gg, bb));
"Inconsistency %3.4e != %3.4e", vol_dot_data(gg, bb),
elemData.dotNormalFace(gg, bb));
}
};
auto check_continuity_of_h1_base = [&](auto &vol_data) {
if (vol_data.size1() != elemData.shapeFunH1Values.size1())
"Inconsistent number of integration points");
if (vol_data.size2() != elemData.shapeFunH1Values.size2())
"Inconsistent number of base functions");
constexpr
double eps = 1e-12;
for (size_t gg = 0; gg != vol_data.size1(); ++gg)
for (size_t bb = 0; bb != vol_data.size2(); ++bb) {
const double error = std::abs(
vol_data(gg, bb) - elemData.shapeFunH1Values(gg, bb));
"Inconsistency %3.4e != %3.4e", vol_data(gg, bb),
elemData.shapeFunH1Values(gg, bb));
}
};
if (elemData.dotNormalEleLeft.size2() != 0)
CHKERR check_continuity_of_base(elemData.dotNormalEleLeft);
else if (elemData.dotNormalEleRight.size2() != 0)
CHKERR check_continuity_of_base(elemData.dotNormalEleRight);
else if (elemData.shapeFunH1VolSide.size2() != 0)
CHKERR check_continuity_of_h1_base(elemData.shapeFunH1VolSide);
}
}
};
new OnContactSideMaster(m_field, elem_data));
new OnContactSideSlave(m_field, elem_data));
}
return 0;
}