13static const double face_coords[4][9] = {{0, 0, 0, 1, 0, 0, 0, 0, 1},
14 {1, 0, 0, 0, 1, 0, 0, 0, 1},
15 {0, 0, 0, 0, 1, 0, 0, 0, 1},
16 {0, 0, 0, 1, 0, 0, 0, 1, 0}};
19 {0, 0, 0, 1, 0, 0}, {1, 0, 0, 0, 1, 0}, {0, 1, 0, 0, 0, 0},
20 {0, 0, 0, 0, 0, 1}, {1, 0, 0, 0, 0, 1}, {0, 1, 0, 0, 0, 1}};
22int main(
int argc,
char *argv[]) {
28 moab::Core mb_instance;
29 moab::Interface &moab = mb_instance;
31 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
33 PetscBool flg = PETSC_TRUE;
35#if PETSC_VERSION_GE(3, 6, 4)
44 if (flg != PETSC_TRUE) {
45 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
59 CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
84 "MESH_NODE_POSITIONS");
91 "MESH_NODE_POSITIONS");
98 "MESH_NODE_POSITIONS");
105 "MESH_NODE_POSITIONS");
138 CHKERR skin.find_skin(0, tets,
false, skin_faces);
146 faces = subtract(faces, skin_faces);
150 CHKERR moab.get_adjacencies(faces, 1,
false, edges, moab::Interface::UNION);
154 CHKERR moab.get_adjacencies(skin_faces, 1,
false, skin_edges,
155 moab::Interface::UNION);
207 CHKERR VecSetRandom(
v, PETSC_NULLPTR);
210 "TEST_PROBLEM",
ROW,
v, INSERT_VALUES, SCATTER_REVERSE);
215 typedef tee_device<std::ostream, std::ofstream>
TeeDevice;
217 std::ofstream ofs(
"forces_and_sources_hdiv_continuity_check.txt");
229 "H1", UserDataOperator::OPROW),
230 mField(m_field), tH(
th) {}
241 boost::shared_ptr<const NumeredEntFiniteElement> mofem_fe =
242 getNumeredEntFiniteElementPtr();
245 .find(boost::make_tuple(type, side))
248 int sense = side_table.get<1>()
249 .find(boost::make_tuple(type, side))
254 int nb_dofs = data.
getN().size2();
255 for (
int dd = 0; dd < nb_dofs; dd++) {
261 (
const void **)&t_ptr);
266 if (type == MBEDGE) {
268 boost::shared_ptr<const NumeredEntFiniteElement> mofem_fe =
269 getNumeredEntFiniteElementPtr();
272 .find(boost::make_tuple(type, side))
276 CHKERR mField.
get_moab().get_adjacencies(&edge, 1, 3,
false, adj_tets,
277 moab::Interface::UNION);
278 const int nb_adj_tets = adj_tets.size();
281 int nb_dofs = data.
getN().size2();
282 for (
int dd = 0; dd < nb_dofs; dd++) {
288 (
const void **)&t_ptr);
290 *t_ptr +=
t / nb_adj_tets;
310 gaussPts.resize(4, 4 + 6);
312 for (; ff < 4; ff++) {
314 for (; dd < 3; dd++) {
316 cblas_ddot(3, &N_tri(0, 0), 1, &
face_coords[ff][dd], 3);
322 for (; ee < 6; ee++) {
324 for (; dd < 3; dd++) {
335 struct OpFacesSkinFluxes
345 "H1", UserDataOperator::OPROW),
346 mField(m_field), tH1(th1), tH2(th2), mySplit(my_split) {}
354 EntityHandle face = getNumeredEntFiniteElementPtr()->getEnt();
358 (
const void **)&t_ptr);
362 (
const void **)&tn_ptr);
368 int nb_dofs = data.
getN().size2();
369 for (; dd < nb_dofs; dd++) {
371 *tn_ptr += -data.
getN()(0, dd) * val;
374 const double eps = 1e-8;
375 if (fabs(*tn_ptr) >
eps) {
377 "H1 continuity failed %6.4e", *tn_ptr);
380 mySplit.precision(5);
381 mySplit << face <<
" " << fabs(*tn_ptr) << std::endl;
397 "H1", UserDataOperator::OPROW),
398 mField(m_field), tH1(_th1), tH2(_th2), mySplit(my_split) {}
407 EntityHandle face = getNumeredEntFiniteElementPtr()->getEnt();
411 (
const void **)&t_ptr);
415 (
const void **)&tn_ptr);
420 const double eps = 1e-8;
421 if (fabs(*tn_ptr) >
eps) {
423 "H1 continuity failed %6.4e", *tn_ptr);
426 mySplit.precision(5);
428 mySplit << face <<
" " << fabs(*tn_ptr) << std::endl;
443 gaussPts.resize(3, 1);
462 "H1", UserDataOperator::OPROW),
463 mField(m_field), tH1(_th1), tH2(_th2), mySplit(my_split) {}
472 EntityHandle edge = getNumeredEntFiniteElementPtr()->getEnt();
476 (
const void **)&t_ptr);
480 (
const void **)&tn_ptr);
486 int nb_dofs = data.
getN().size2();
488 for (; dd < nb_dofs; dd++) {
490 tn += data.
getN()(0, dd) * val;
495 mySplit.precision(5);
497 mySplit << edge <<
" " << fabs(*tn_ptr) << std::endl;
512 gaussPts.resize(2, 1);
520 MyTetFE tet_fe(m_field);
521 MyTriFE tri_fe(m_field);
522 MyTriFE skin_fe(m_field);
523 MyEdgeFE edge_fe(m_field);
526 double def_val[] = {0, 0, 0};
527 CHKERR moab.tag_get_handle(
"T", 3, MB_TYPE_DOUBLE, th1,
528 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
530 tet_fe.getOpPtrVector().push_back(
532 tet_fe.getOpPtrVector().push_back(
new OpTetFluxes(m_field, th1));
535 CHKERR moab.tag_get_handle(
"TN", 1, MB_TYPE_DOUBLE, th2,
536 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
538 tri_fe.getOpPtrVector().push_back(
540 tri_fe.getOpPtrVector().push_back(
542 tri_fe.getOpPtrVector().push_back(
544 skin_fe.getOpPtrVector().push_back(
545 new OpFacesSkinFluxes(m_field, th1, th2, my_split));
547 edge_fe.getOpPtrVector().push_back(
549 edge_fe.getOpPtrVector().push_back(
551 edge_fe.getOpPtrVector().push_back(
552 new OpEdgesFluxes(m_field, th1, th2, my_split));
554 for (Range::iterator fit = faces.begin(); fit != faces.end(); fit++) {
555 CHKERR moab.tag_set_data(th1, &*fit, 1, &def_val);
557 CHKERR moab.tag_set_data(th2, &*fit, 1, &def_val);
562 my_split <<
"internal\n";
564 my_split <<
"skin\n";
566 my_split <<
"edges\n";
571 CHKERR moab.create_meshset(MESHSET_SET, meshset);
575 CHKERR moab.write_file(
"out.vtk",
"VTK",
"", &meshset, 1);