5 namespace bio = boost::iostreams;
11 static char help[] =
"...\n\n";
13 static 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}};
22 int main(
int argc,
char *argv[]) {
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_NULL);
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");
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;
301 int getRule(
int order) {
return -1; };
314 gaussPts.resize(4, 4 + 6);
316 for (; ff < 4; ff++) {
318 for (;
dd < 3;
dd++) {
326 for (; ee < 6; ee++) {
328 for (;
dd < 3;
dd++) {
329 gaussPts(
dd, 4 + ee) =
332 gaussPts(3, 4 + ee) = 1;
337 }
catch (std::exception &ex) {
338 std::ostringstream ss;
339 ss <<
"thorw in method: " << ex.what() <<
" at line " << __LINE__
340 <<
" in file " << __FILE__;
348 struct OpFacesSkinFluxes
359 mField(m_field), tH1(th1), tH2(th2), mySplit(my_split) {}
367 EntityHandle face = getNumeredEntFiniteElementPtr()->getEnt();
371 (
const void **)&t_ptr);
375 (
const void **)&tn_ptr);
381 int nb_dofs = data.
getN().size2();
382 for (;
dd < nb_dofs;
dd++) {
384 *tn_ptr += -data.
getN()(0,
dd) * val;
387 const double eps = 1e-8;
388 if (fabs(*tn_ptr) >
eps) {
390 "H1 continuity failed %6.4e", *tn_ptr);
393 mySplit.precision(5);
394 mySplit << face <<
" " << fabs(*tn_ptr) << std::endl;
411 mField(m_field), tH1(_th1), tH2(_th2), mySplit(my_split) {}
420 EntityHandle face = getNumeredEntFiniteElementPtr()->getEnt();
424 (
const void **)&t_ptr);
428 (
const void **)&tn_ptr);
433 const double eps = 1e-8;
434 if (fabs(*tn_ptr) >
eps) {
436 "H1 continuity failed %6.4e", *tn_ptr);
439 mySplit.precision(5);
441 mySplit << face <<
" " << fabs(*tn_ptr) << std::endl;
451 int getRule(
int order) {
return -1; };
456 gaussPts.resize(3, 1);
476 mField(m_field), tH1(_th1), tH2(_th2), mySplit(my_split) {}
485 EntityHandle edge = getNumeredEntFiniteElementPtr()->getEnt();
489 (
const void **)&t_ptr);
493 (
const void **)&tn_ptr);
499 int nb_dofs = data.
getN().size2();
501 for (;
dd < nb_dofs;
dd++) {
503 tn += data.
getN()(0,
dd) * val;
523 mySplit.precision(5);
525 mySplit << edge <<
" " << fabs(*tn_ptr) << std::endl;
535 int getRule(
int order) {
return -1; };
540 gaussPts.resize(2, 1);
541 gaussPts(0, 0) = 0.5;
548 MyTetFE tet_fe(m_field);
549 MyTriFE tri_fe(m_field);
550 MyTriFE skin_fe(m_field);
551 MyEdgeFE edge_fe(m_field);
554 double def_val[] = {0, 0, 0};
555 CHKERR moab.tag_get_handle(
"T", 3, MB_TYPE_DOUBLE, th1,
556 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
558 tet_fe.getOpPtrVector().push_back(
560 tet_fe.getOpPtrVector().push_back(
new OpTetFluxes(m_field, th1));
563 CHKERR moab.tag_get_handle(
"TN", 1, MB_TYPE_DOUBLE, th2,
564 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
566 tri_fe.getOpPtrVector().push_back(
568 tri_fe.getOpPtrVector().push_back(
570 tri_fe.getOpPtrVector().push_back(
572 skin_fe.getOpPtrVector().push_back(
573 new OpFacesSkinFluxes(m_field, th1, th2, my_split));
575 edge_fe.getOpPtrVector().push_back(
577 edge_fe.getOpPtrVector().push_back(
579 edge_fe.getOpPtrVector().push_back(
580 new OpEdgesFluxes(m_field, th1, th2, my_split));
582 for (Range::iterator fit = faces.begin(); fit != faces.end(); fit++) {
583 CHKERR moab.tag_set_data(th1, &*fit, 1, &def_val);
585 CHKERR moab.tag_set_data(th2, &*fit, 1, &def_val);
590 my_split <<
"internal\n";
592 my_split <<
"skin\n";
594 my_split <<
"edges\n";
599 CHKERR moab.create_meshset(MESHSET_SET, meshset);
603 CHKERR moab.write_file(
"out.vtk",
"VTK",
"", &meshset, 1);